From 167655d01d1653727f795fdacb1178cf610e806d Mon Sep 17 00:00:00 2001
From: Stefan Kroboth <stefan.kroboth@gmail.com>
Date: Sun, 6 Oct 2019 08:54:28 +0200
Subject: [PATCH 1/3] Update ndarray to version 0.13
---
Cargo.toml | 2 +-
src/convert.rs | 4 ++--
src/eigh.rs | 7 ++++---
src/krylov/arnoldi.rs | 4 ++--
src/krylov/householder.rs | 7 ++++---
src/krylov/mgs.rs | 4 ++--
src/layout.rs | 8 ++++----
src/operator.rs | 2 +-
src/qr.rs | 4 ++--
src/solveh.rs | 2 +-
src/svd.rs | 2 +-
src/svddc.rs | 22 +++++-----------------
tests/det.rs | 8 ++++----
tests/triangular.rs | 2 +-
14 files changed, 34 insertions(+), 44 deletions(-)
diff --git a/Cargo.toml b/Cargo.toml
index 6686b0ae..3559dfec 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -30,7 +30,7 @@ num-complex = "0.2.1"
rand = "0.5"
[dependencies.ndarray]
-version = "0.12"
+version = "0.13"
features = ["blas"]
default-features = false
diff --git a/src/convert.rs b/src/convert.rs
index 1562650f..92b0e99e 100644
--- a/src/convert.rs
+++ b/src/convert.rs
@@ -113,14 +113,14 @@ where
assert!(a.is_square());
match uplo {
UPLO::Upper => {
- for row in 0..a.rows() {
+ for row in 0..a.nrows() {
for col in 0..row {
a[(row, col)] = a[(col, row)].conj();
}
}
}
UPLO::Lower => {
- for col in 0..a.cols() {
+ for col in 0..a.ncols() {
for row in 0..col {
a[(row, col)] = a[(col, row)].conj();
}
diff --git a/src/eigh.rs b/src/eigh.rs
index 967d7bf0..5e2162b3 100644
--- a/src/eigh.rs
+++ b/src/eigh.rs
@@ -8,6 +8,7 @@ use crate::layout::*;
use crate::operator::LinearOperator;
use crate::types::*;
use crate::UPLO;
+use std::iter::FromIterator;
/// Eigenvalue decomposition of Hermite matrix reference
pub trait Eigh {
@@ -70,7 +71,7 @@ where
MatrixLayout::F(_) => {}
}
let s = unsafe { A::eigh(true, self.square_layout()?, uplo, self.as_allocated_mut()?)? };
- Ok((ArrayBase::from_vec(s), self))
+ Ok((ArrayBase::from(s), self))
}
}
@@ -126,7 +127,7 @@ where
fn eigvalsh_inplace(&mut self, uplo: UPLO) -> Result<Self::EigVal> {
let s = unsafe { A::eigh(true, self.square_layout()?, uplo, self.as_allocated_mut()?)? };
- Ok(ArrayBase::from_vec(s))
+ Ok(ArrayBase::from(s))
}
}
@@ -164,7 +165,7 @@ where
fn ssqrt_into(self, uplo: UPLO) -> Result<Self::Output> {
let (e, v) = self.eigh_into(uplo)?;
- let e_sqrt = Array1::from_iter(e.iter().map(|r| Scalar::from_real(r.sqrt())));
+ let e_sqrt = Array::from_iter(e.iter().map(|r| Scalar::from_real(r.sqrt())));
let ev = e_sqrt.into_diagonal().apply2(&v.t());
Ok(v.apply2(&ev))
}
diff --git a/src/krylov/arnoldi.rs b/src/krylov/arnoldi.rs
index 5246a59b..297cafbf 100644
--- a/src/krylov/arnoldi.rs
+++ b/src/krylov/arnoldi.rs
@@ -38,7 +38,7 @@ where
assert!(ortho.tolerance() < One::one());
// normalize before append because |v| may be smaller than ortho.tolerance()
let norm = v.norm_l2();
- azip!(mut v(&mut v) in { *v = v.div_real(norm) });
+ azip!((v in &mut v) *v = v.div_real(norm));
ortho.append(v.view());
Arnoldi {
a,
@@ -82,7 +82,7 @@ where
self.a.apply_mut(&mut self.v);
let result = self.ortho.div_append(&mut self.v);
let norm = self.v.norm_l2();
- azip!(mut v(&mut self.v) in { *v = v.div_real(norm) });
+ azip!((v in &mut self.v) *v = v.div_real(norm));
match result {
AppendResult::Added(coef) => {
self.h.push(coef.clone());
diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs
index 17b7d905..90ac58c1 100644
--- a/src/krylov/householder.rs
+++ b/src/krylov/householder.rs
@@ -18,7 +18,8 @@ where
let alpha = -x[0].mul_real(norm / x[0].abs());
x[0] -= alpha;
let inv_rev_norm = A::Real::one() / x.norm_l2();
- azip!(mut a(x) in { *a = a.mul_real(inv_rev_norm)});
+ // azip!(mut a(x) in { *a = a.mul_real(inv_rev_norm)});
+ azip!((a in x) *a = a.mul_real(inv_rev_norm));
}
/// Take a reflection `P = I - 2ww^T`
@@ -107,7 +108,7 @@ impl<A: Scalar + Lapack> Householder<A> {
let k = self.len();
let res = a.slice(s![k..]).norm_l2();
let mut c = Array1::zeros(k + 1);
- azip!(mut c(c.slice_mut(s![..k])), a(a.slice(s![..k])) in { *c = a });
+ azip!((c in c.slice_mut(s![..k]), &a in a.slice(s![..k])) *c = a);
if k < a.len() {
let ak = a[k];
c[k] = -ak.mul_real(res / ak.abs());
@@ -123,7 +124,7 @@ impl<A: Scalar + Lapack> Householder<A> {
S: DataMut<Elem = A>,
{
let k = self.len();
- azip!(mut a( a.slice_mut(s![..k])) in { *a = A::zero() });
+ azip!((a in a.slice_mut(s![..k])) *a = A::zero());
self.backward_reflection(a);
}
}
diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs
index b6a8d2bf..2e2e1fef 100644
--- a/src/krylov/mgs.rs
+++ b/src/krylov/mgs.rs
@@ -51,7 +51,7 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
for i in 0..self.len() {
let q = &self.q[i];
let c = q.inner(&a);
- azip!(mut a (&mut *a), q (q) in { *a = *a - c * q } );
+ azip!((a in &mut *a, &q in q) *a = *a - c * q);
coef[i] = c;
}
let nrm = a.norm_l2();
@@ -88,7 +88,7 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
// Linearly dependent
return AppendResult::Dependent(coef);
}
- azip!(mut a(&mut *a) in { *a = *a / A::from_real(nrm) });
+ azip!((a in &mut *a) *a = *a / A::from_real(nrm));
self.q.push(a.to_owned());
AppendResult::Added(coef)
}
diff --git a/src/layout.rs b/src/layout.rs
index 730756fe..51099e6f 100644
--- a/src/layout.rs
+++ b/src/layout.rs
@@ -96,10 +96,10 @@ where
let shape = self.shape();
let strides = self.strides();
if shape[0] == strides[1] as usize {
- return Ok(MatrixLayout::F((self.cols() as i32, self.rows() as i32)));
+ return Ok(MatrixLayout::F((self.ncols() as i32, self.nrows() as i32)));
}
if shape[1] == strides[0] as usize {
- return Ok(MatrixLayout::C((self.rows() as i32, self.cols() as i32)));
+ return Ok(MatrixLayout::C((self.nrows() as i32, self.ncols() as i32)));
}
Err(LinalgError::InvalidStride {
s0: strides[0],
@@ -122,8 +122,8 @@ where
Ok(())
} else {
Err(LinalgError::NotSquare {
- rows: self.rows() as i32,
- cols: self.cols() as i32,
+ rows: self.nrows() as i32,
+ cols: self.ncols() as i32,
})
}
}
diff --git a/src/operator.rs b/src/operator.rs
index a39c43d1..368650bf 100644
--- a/src/operator.rs
+++ b/src/operator.rs
@@ -25,7 +25,7 @@ pub trait LinearOperator {
S: DataMut<Elem = Self::Elem>,
{
let b = self.apply(a);
- azip!(mut a(a), b in { *a = b });
+ azip!((a in a, &b in &b) *a = b);
}
/// Apply operator with move
diff --git a/src/qr.rs b/src/qr.rs
index 65443ddc..be2de0c2 100644
--- a/src/qr.rs
+++ b/src/qr.rs
@@ -103,8 +103,8 @@ where
type R = Array2<A>;
fn qr_into(mut self) -> Result<(Self::Q, Self::R)> {
- let n = self.rows();
- let m = self.cols();
+ let n = self.nrows();
+ let m = self.ncols();
let k = ::std::cmp::min(n, m);
let l = self.layout()?;
let r = unsafe { A::qr(l, self.as_allocated_mut()?)? };
diff --git a/src/solveh.rs b/src/solveh.rs
index a54caa54..582eeb26 100644
--- a/src/solveh.rs
+++ b/src/solveh.rs
@@ -309,7 +309,7 @@ where
let mut ln_det = A::Real::zero();
let mut ipiv_enum = ipiv_iter.enumerate();
while let Some((k, ipiv_k)) = ipiv_enum.next() {
- debug_assert!(k < a.rows() && k < a.cols());
+ debug_assert!(k < a.nrows() && k < a.ncols());
if ipiv_k > 0 {
// 1x1 block at k, must be real.
let elem = unsafe { a.uget((k, k)) }.re();
diff --git a/src/svd.rs b/src/svd.rs
index 3acce0c2..7ad2dc8f 100644
--- a/src/svd.rs
+++ b/src/svd.rs
@@ -81,7 +81,7 @@ where
let vt = svd_res
.vt
.map(|vt| into_matrix(l.resized(m, m), vt).expect("Size of VT mismatches"));
- let s = ArrayBase::from_vec(svd_res.s);
+ let s = ArrayBase::from(svd_res.s);
Ok((u, s, vt))
}
}
diff --git a/src/svddc.rs b/src/svddc.rs
index 10a49668..82f7fe5d 100644
--- a/src/svddc.rs
+++ b/src/svddc.rs
@@ -28,10 +28,7 @@ pub trait SVDDCInto {
type U;
type VT;
type Sigma;
- fn svddc_into(
- self,
- uvt_flag: UVTFlag,
- ) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)>;
+ fn svddc_into(self, uvt_flag: UVTFlag) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)>;
}
/// Singular-value decomposition of matrix reference by divide-and-conquer
@@ -39,10 +36,7 @@ pub trait SVDDCInplace {
type U;
type VT;
type Sigma;
- fn svddc_inplace(
- &mut self,
- uvt_flag: UVTFlag,
- ) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)>;
+ fn svddc_inplace(&mut self, uvt_flag: UVTFlag) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)>;
}
impl<A, S> SVDDC for ArrayBase<S, Ix2>
@@ -68,10 +62,7 @@ where
type VT = Array2<A>;
type Sigma = Array1<A::Real>;
- fn svddc_into(
- mut self,
- uvt_flag: UVTFlag,
- ) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)> {
+ fn svddc_into(mut self, uvt_flag: UVTFlag) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)> {
self.svddc_inplace(uvt_flag)
}
}
@@ -85,10 +76,7 @@ where
type VT = Array2<A>;
type Sigma = Array1<A::Real>;
- fn svddc_inplace(
- &mut self,
- uvt_flag: UVTFlag,
- ) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)> {
+ fn svddc_inplace(&mut self, uvt_flag: UVTFlag) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)> {
let l = self.layout()?;
let svd_res = unsafe { A::svddc(l, uvt_flag, self.as_allocated_mut()?)? };
let (m, n) = l.size();
@@ -104,7 +92,7 @@ where
let vt = svd_res
.vt
.map(|vt| into_matrix(l.resized(ldvt, tdvt), vt).expect("Size of VT mismatches"));
- let s = ArrayBase::from_vec(svd_res.s);
+ let s = ArrayBase::from(svd_res.s);
Ok((u, s, vt))
}
}
diff --git a/tests/det.rs b/tests/det.rs
index 065685bc..610b0b59 100644
--- a/tests/det.rs
+++ b/tests/det.rs
@@ -8,9 +8,9 @@ where
A: Scalar,
S: Data<Elem = A>,
{
- let mut select_rows = (0..a.rows()).collect::<Vec<_>>();
+ let mut select_rows = (0..a.nrows()).collect::<Vec<_>>();
select_rows.remove(row);
- let mut select_cols = (0..a.cols()).collect::<Vec<_>>();
+ let mut select_cols = (0..a.ncols()).collect::<Vec<_>>();
select_cols.remove(col);
a.select(Axis(0), &select_rows).select(Axis(1), &select_cols)
}
@@ -24,8 +24,8 @@ where
A: Scalar,
S: Data<Elem = A>,
{
- assert_eq!(a.rows(), a.cols());
- match a.cols() {
+ assert_eq!(a.nrows(), a.ncols());
+ match a.ncols() {
0 => A::one(),
1 => a[(0, 0)],
cols => (0..cols)
diff --git a/tests/triangular.rs b/tests/triangular.rs
index 21f12c71..ca609c5c 100644
--- a/tests/triangular.rs
+++ b/tests/triangular.rs
@@ -20,7 +20,7 @@ fn test2d<A, Sa, Sb>(uplo: UPLO, a: &ArrayBase<Sa, Ix2>, b: &ArrayBase<Sb, Ix2>,
where
A: Scalar + Lapack,
Sa: Data<Elem = A>,
- Sb: DataMut<Elem = A> + DataOwned + DataClone,
+ Sb: DataMut<Elem = A> + DataOwned + Data + RawDataClone,
{
println!("a = {:?}", a);
println!("b = {:?}", b);
From d8d635a9b80e4cdfb9c06b08441e29ae1b161a37 Mon Sep 17 00:00:00 2001
From: Stefan Kroboth <stefan.kroboth@gmail.com>
Date: Wed, 9 Oct 2019 07:50:58 +0200
Subject: [PATCH 2/3] Apply suggestions from code review
Co-Authored-By: Jim Turner <github@turner.link>
---
src/krylov/householder.rs | 1 -
src/krylov/mgs.rs | 4 ++--
2 files changed, 2 insertions(+), 3 deletions(-)
diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs
index 90ac58c1..9515e7b1 100644
--- a/src/krylov/householder.rs
+++ b/src/krylov/householder.rs
@@ -18,7 +18,6 @@ where
let alpha = -x[0].mul_real(norm / x[0].abs());
x[0] -= alpha;
let inv_rev_norm = A::Real::one() / x.norm_l2();
- // azip!(mut a(x) in { *a = a.mul_real(inv_rev_norm)});
azip!((a in x) *a = a.mul_real(inv_rev_norm));
}
diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs
index 2e2e1fef..589df726 100644
--- a/src/krylov/mgs.rs
+++ b/src/krylov/mgs.rs
@@ -51,7 +51,7 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
for i in 0..self.len() {
let q = &self.q[i];
let c = q.inner(&a);
- azip!((a in &mut *a, &q in q) *a = *a - c * q);
+ azip!((a in a, &q in q) *a = *a - c * q);
coef[i] = c;
}
let nrm = a.norm_l2();
@@ -88,7 +88,7 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
// Linearly dependent
return AppendResult::Dependent(coef);
}
- azip!((a in &mut *a) *a = *a / A::from_real(nrm));
+ azip!((a in a) *a = *a / A::from_real(nrm));
self.q.push(a.to_owned());
AppendResult::Added(coef)
}
From a482367925a84388242608072ee393f61c072ebe Mon Sep 17 00:00:00 2001
From: Stefan Kroboth <stefan.kroboth@gmail.com>
Date: Wed, 9 Oct 2019 08:22:26 +0200
Subject: [PATCH 3/3] revert back to &mut *a in azip! macros of krylov/mgs.rs
---
src/krylov/mgs.rs | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs
index 589df726..2e2e1fef 100644
--- a/src/krylov/mgs.rs
+++ b/src/krylov/mgs.rs
@@ -51,7 +51,7 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
for i in 0..self.len() {
let q = &self.q[i];
let c = q.inner(&a);
- azip!((a in a, &q in q) *a = *a - c * q);
+ azip!((a in &mut *a, &q in q) *a = *a - c * q);
coef[i] = c;
}
let nrm = a.norm_l2();
@@ -88,7 +88,7 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
// Linearly dependent
return AppendResult::Dependent(coef);
}
- azip!((a in a) *a = *a / A::from_real(nrm));
+ azip!((a in &mut *a) *a = *a / A::from_real(nrm));
self.q.push(a.to_owned());
AppendResult::Added(coef)
}