Skip to content

Commit 9bd1409

Browse files
doraneko94doraneko94
committed
add eig-test, fix calc method (eig)
1 parent 8e1191d commit 9bd1409

File tree

5 files changed

+155
-36
lines changed

5 files changed

+155
-36
lines changed

‎examples/eig.rs

Lines changed: 4 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,12 @@
1-
extern crate ndarray;
2-
extern crate ndarray_linalg;
3-
41
use ndarray::*;
52
use ndarray_linalg::*;
63

74
fn main() {
8-
let a = arr2(&[[0.0, -1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 3.0, 2.0]]);
5+
let a = arr2(&[[2.0, 1.0, 2.0], [-2.0, 2.0, 1.0], [1.0, 2.0, -2.0]]);
96
let (e, vecs) = a.clone().eig().unwrap();
107
println!("eigenvalues = \n{:?}", e);
11-
println!("V = \n{:?}", vecs.t());
12-
let av = a.dot(&vecs.t());
8+
println!("V = \n{:?}", vecs);
9+
let a_c: Array2<c64> = a.map(|f| c64::new(*f, 0.0));
10+
let av = a_c.dot(&vecs);
1311
println!("AV = \n{:?}", av);
1412
}

‎src/eig.rs

Lines changed: 2 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -18,19 +18,13 @@ where
1818
S: Data<Elem = A>,
1919
{
2020
type EigVal = Array1<A::Complex>;
21-
type EigVec = Array2<A>;
21+
type EigVec = Array2<A::Complex>;
2222

2323
fn eig(&self) -> Result<(Self::EigVal, Self::EigVec)> {
2424
let mut a = self.to_owned();
2525
let layout = a.square_layout()?;
26-
// XXX Force layout to be Fortran (see #146)
27-
match layout {
28-
MatrixLayout::C(_) => a.swap_axes(0, 1),
29-
MatrixLayout::F(_) => {}
30-
}
31-
let (s, t) = unsafe { A::eig(true, a.square_layout()?, a.as_allocated_mut()?)? };
26+
let (s, t) = unsafe { A::eig(true, layout, a.as_allocated_mut()?)? };
3227
let (n, _) = layout.size();
33-
println!("{:?}", t);
3428
Ok((ArrayBase::from(s), ArrayBase::from(t).into_shape((n as usize, n as usize)).unwrap()))
3529
}
3630
}
@@ -50,12 +44,6 @@ where
5044

5145
fn eigvals(&self) -> Result<Self::EigVal> {
5246
let mut a = self.to_owned();
53-
let layout = a.square_layout()?;
54-
// XXX Force layout to be Fortran (see #146)
55-
match layout {
56-
MatrixLayout::C(_) => a.swap_axes(0, 1),
57-
MatrixLayout::F(_) => {}
58-
}
5947
let (s, _) = unsafe { A::eig(true, a.square_layout()?, a.as_allocated_mut()?)? };
6048
Ok(ArrayBase::from(s))
6149
}

‎src/lapack/eig.rs

Lines changed: 36 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -11,44 +11,64 @@ use super::into_result;
1111

1212
/// Wraps `*geev` for real/complex
1313
pub trait Eig_: Scalar {
14-
unsafe fn eig(calc_v: bool, l: MatrixLayout, a: &mut [Self]) -> Result<(Vec<Self::Complex>, Vec<Self>)>;
14+
unsafe fn eig(calc_v: bool, l: MatrixLayout, a: &mut [Self]) -> Result<(Vec<Self::Complex>, Vec<Self::Complex>)>;
1515
}
1616

1717
macro_rules! impl_eig_complex {
1818
($scalar:ty, $ev:path) => {
1919
impl Eig_ for $scalar {
2020
unsafe fn eig(calc_v: bool, l: MatrixLayout, mut a: &mut [Self]) -> Result<(Vec<Self::Complex>, Vec<Self::Complex>)> {
2121
let (n, _) = l.size();
22-
let jobvl = if calc_v { b'V' } else { b'N' };
22+
let jobvr = if calc_v { b'V' } else { b'N' };
2323
let mut w = vec![Self::Complex::zero(); n as usize];
24-
let mut vl = vec![Self::Complex::zero(); (n * n) as usize];
25-
let mut vr = Vec::new();
26-
let info = $ev(l.lapacke_layout(), jobvl, b'N', n, &mut a, n, &mut w, &mut vl, n, &mut vr, n);
27-
into_result(info, (w, vl))
24+
let mut vl = Vec::new();
25+
let mut vr = vec![Self::Complex::zero(); (n * n) as usize];
26+
let info = $ev(l.lapacke_layout(), b'N', jobvr, n, &mut a, n, &mut w, &mut vl, n, &mut vr, n);
27+
into_result(info, (w, vr))
2828
}
2929
}
3030
};
3131
}
3232

3333
macro_rules! impl_eig_real {
34-
($scalar:ty, $ev:path, $cn:ty) => {
34+
($scalar:ty, $ev:path) => {
3535
impl Eig_ for $scalar {
36-
unsafe fn eig(calc_v: bool, l: MatrixLayout, mut a: &mut [Self]) -> Result<(Vec<Self::Complex>, Vec<Self::Real>)> {
36+
unsafe fn eig(calc_v: bool, l: MatrixLayout, mut a: &mut [Self]) -> Result<(Vec<Self::Complex>, Vec<Self::Complex>)> {
3737
let (n, _) = l.size();
38-
let jobvl = if calc_v { b'V' } else { b'N' };
38+
let jobvr = if calc_v { b'V' } else { b'N' };
3939
let mut wr = vec![Self::Real::zero(); n as usize];
4040
let mut wi = vec![Self::Real::zero(); n as usize];
41-
let mut vl = vec![Self::Real::zero(); (n * n) as usize];
42-
let mut vr = Vec::new();
43-
let info = $ev(l.lapacke_layout(), jobvl, b'N', n, &mut a, n, &mut wr, &mut wi, &mut vl, n, &mut vr, n);
44-
let w: Vec<$cn> = wr.iter().zip(wi.iter()).map(|(&r, &i)| <$cn>::new(r, i)).collect();
45-
into_result(info, (w, vl))
41+
let mut vl = Vec::new();
42+
let mut vr = vec![Self::Real::zero(); (n * n) as usize];
43+
let info = $ev(l.lapacke_layout(), b'N', jobvr, n, &mut a, n, &mut wr, &mut wi, &mut vl, n, &mut vr, n);
44+
let w: Vec<Self::Complex> = wr.iter().zip(wi.iter()).map(|(&r, &i)| Self::Complex::new(r, i)).collect();
45+
let n = n as usize;
46+
let mut conj = false;
47+
let mut v = vec![Self::Complex::zero(); n * n];
48+
for (i, c) in w.iter().enumerate() {
49+
if conj {
50+
for j in 0..n {
51+
v[n*j+i] = Self::Complex::new(vr[n*j+i-1], -vr[n*j+i]);
52+
}
53+
conj = false;
54+
} else if c.im != 0.0 {
55+
conj = true;
56+
for j in 0..n {
57+
v[n*j+i] = Self::Complex::new(vr[n*j+i], vr[n*j+i+1]);
58+
}
59+
} else {
60+
for j in 0..n {
61+
v[n*j+i] = Self::Complex::new(vr[n*j+i], 0.0);
62+
}
63+
}
64+
}
65+
into_result(info, (w, v))
4666
}
4767
}
4868
};
4969
}
5070

51-
impl_eig_real!(f64, lapacke::dgeev, c64);
52-
impl_eig_real!(f32, lapacke::sgeev, c32);
71+
impl_eig_real!(f64, lapacke::dgeev);
72+
impl_eig_real!(f32, lapacke::sgeev);
5373
impl_eig_complex!(c64, lapacke::zgeev);
5474
impl_eig_complex!(c32, lapacke::cgeev);

‎src/lib.rs

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
//! - Decomposition methods:
88
//! - [QR decomposition](qr/index.html)
99
//! - [Cholesky/LU decomposition](cholesky/index.html)
10+
//! - [Eigenvalue decomposition](eig/index.html)
1011
//! - [Eigenvalue decomposition for Hermite matrices](eigh/index.html)
1112
//! - [**S**ingular **V**alue **D**ecomposition](svd/index.html)
1213
//! - Solution of linear systems:

‎tests/eig.rs

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
use ndarray::*;
2+
use ndarray_linalg::*;
3+
4+
#[test]
5+
fn dgeev() {
6+
// https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/dgeev_ex.f.htm
7+
let a: Array2<f64> = arr2(&[[-1.01, 0.86, -4.60, 3.31, -4.81],
8+
[ 3.98, 0.53, -7.04, 5.29, 3.55],
9+
[ 3.30, 8.26, -3.89, 8.20, -1.51],
10+
[ 4.43, 4.96, -7.66, -7.33, 6.18],
11+
[ 7.31, -6.43, -6.16, 2.47, 5.58]]);
12+
let (e, vecs): (Array1<_>, Array2<_>) = (&a).eig().unwrap();
13+
assert_close_l2!(&e,
14+
&arr1(&[c64::new( 2.86, 10.76), c64::new( 2.86,-10.76), c64::new( -0.69, 4.70), c64::new( -0.69, -4.70), c64::new(-10.46, 0.00)]),
15+
1.0e-3);
16+
17+
/*
18+
let answer = &arr2(&[[c64::new( 0.11, 0.17), c64::new( 0.11, -0.17), c64::new( 0.73, 0.00), c64::new( 0.73, 0.00), c64::new( 0.46, 0.00)],
19+
[c64::new( 0.41, -0.26), c64::new( 0.41, 0.26), c64::new( -0.03, -0.02), c64::new( -0.03, 0.02), c64::new( 0.34, 0.00)],
20+
[c64::new( 0.10, -0.51), c64::new( 0.10, 0.51), c64::new( 0.19, -0.29), c64::new( 0.19, 0.29), c64::new( 0.31, 0.00)],
21+
[c64::new( 0.40, -0.09), c64::new( 0.40, 0.09), c64::new( -0.08, -0.08), c64::new( -0.08, 0.08), c64::new( -0.74, 0.00)],
22+
[c64::new( 0.54, 0.00), c64::new( 0.54, 0.00), c64::new( -0.29, -0.49), c64::new( -0.29, 0.49), c64::new( 0.16, 0.00)]]);
23+
*/
24+
25+
let a_c: Array2<c64> = a.map(|f| c64::new(*f, 0.0));
26+
for (i, v) in vecs.axis_iter(Axis(1)).enumerate() {
27+
let av = a_c.dot(&v);
28+
let ev = v.mapv(|f| e[i] * f);
29+
assert_close_l2!(&av, &ev, 1.0e-7);
30+
}
31+
}
32+
33+
#[test]
34+
fn fgeev() {
35+
// https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/dgeev_ex.f.htm
36+
let a: Array2<f32> = arr2(&[[-1.01, 0.86, -4.60, 3.31, -4.81],
37+
[ 3.98, 0.53, -7.04, 5.29, 3.55],
38+
[ 3.30, 8.26, -3.89, 8.20, -1.51],
39+
[ 4.43, 4.96, -7.66, -7.33, 6.18],
40+
[ 7.31, -6.43, -6.16, 2.47, 5.58]]);
41+
let (e, vecs): (Array1<_>, Array2<_>) = (&a).eig().unwrap();
42+
assert_close_l2!(&e,
43+
&arr1(&[c32::new( 2.86, 10.76), c32::new( 2.86,-10.76), c32::new( -0.69, 4.70), c32::new( -0.69, -4.70), c32::new(-10.46, 0.00)]),
44+
1.0e-3);
45+
46+
/*
47+
let answer = &arr2(&[[c32::new( 0.11, 0.17), c32::new( 0.11, -0.17), c32::new( 0.73, 0.00), c32::new( 0.73, 0.00), c32::new( 0.46, 0.00)],
48+
[c32::new( 0.41, -0.26), c32::new( 0.41, 0.26), c32::new( -0.03, -0.02), c32::new( -0.03, 0.02), c32::new( 0.34, 0.00)],
49+
[c32::new( 0.10, -0.51), c32::new( 0.10, 0.51), c32::new( 0.19, -0.29), c32::new( 0.19, 0.29), c32::new( 0.31, 0.00)],
50+
[c32::new( 0.40, -0.09), c32::new( 0.40, 0.09), c32::new( -0.08, -0.08), c32::new( -0.08, 0.08), c32::new( -0.74, 0.00)],
51+
[c32::new( 0.54, 0.00), c32::new( 0.54, 0.00), c32::new( -0.29, -0.49), c32::new( -0.29, 0.49), c32::new( 0.16, 0.00)]]);
52+
*/
53+
54+
let a_c: Array2<c32> = a.map(|f| c32::new(*f, 0.0));
55+
for (i, v) in vecs.axis_iter(Axis(1)).enumerate() {
56+
let av = a_c.dot(&v);
57+
let ev = v.mapv(|f| e[i] * f);
58+
assert_close_l2!(&av, &ev, 1.0e-5);
59+
}
60+
}
61+
62+
#[test]
63+
fn zgeev() {
64+
// https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/zgeev_ex.f.htm
65+
let a: Array2<c64> = arr2(&[[c64::new( -3.84, 2.25), c64::new( -8.94, -4.75), c64::new( 8.95, -6.53), c64::new( -9.87, 4.82)],
66+
[c64::new( -0.66, 0.83), c64::new( -4.40, -3.82), c64::new( -3.50, -4.26), c64::new( -3.15, 7.36)],
67+
[c64::new( -3.99, -4.73), c64::new( -5.88, -6.60), c64::new( -3.36, -0.40), c64::new( -0.75, 5.23)],
68+
[c64::new( 7.74, 4.18), c64::new( 3.66, -7.53), c64::new( 2.58, 3.60), c64::new( 4.59, 5.41)],]);
69+
let (e, vecs): (Array1<_>, Array2<_>) = (&a).eig().unwrap();
70+
assert_close_l2!(&e,
71+
&arr1(&[c64::new( -9.43,-12.98), c64::new( -3.44, 12.69), c64::new( 0.11, -3.40), c64::new( 5.76, 7.13)]),
72+
1.0e-3);
73+
74+
/*
75+
let answer = &arr2(&[[c64::new( 0.43, 0.33), c64::new( 0.83, 0.00), c64::new( 0.60, 0.00), c64::new( -0.31, 0.03)],
76+
[c64::new( 0.51, -0.03), c64::new( 0.08, -0.25), c64::new( -0.40, -0.20), c64::new( 0.04, 0.34)],
77+
[c64::new( 0.62, 0.00), c64::new( -0.25, 0.28), c64::new( -0.09, -0.48), c64::new( 0.36, 0.06)],
78+
[c64::new( -0.23, 0.11), c64::new( -0.10, -0.32), c64::new( -0.43, 0.13), c64::new( 0.81, 0.00)]]);
79+
*/
80+
81+
for (i, v) in vecs.axis_iter(Axis(1)).enumerate() {
82+
let av = a.dot(&v);
83+
let ev = v.mapv(|f| e[i] * f);
84+
assert_close_l2!(&av, &ev, 1.0e-7);
85+
}
86+
}
87+
88+
#[test]
89+
fn cgeev() {
90+
// https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/zgeev_ex.f.htm
91+
let a: Array2<c32> = arr2(&[[c32::new( -3.84, 2.25), c32::new( -8.94, -4.75), c32::new( 8.95, -6.53), c32::new( -9.87, 4.82)],
92+
[c32::new( -0.66, 0.83), c32::new( -4.40, -3.82), c32::new( -3.50, -4.26), c32::new( -3.15, 7.36)],
93+
[c32::new( -3.99, -4.73), c32::new( -5.88, -6.60), c32::new( -3.36, -0.40), c32::new( -0.75, 5.23)],
94+
[c32::new( 7.74, 4.18), c32::new( 3.66, -7.53), c32::new( 2.58, 3.60), c32::new( 4.59, 5.41)],]);
95+
let (e, vecs): (Array1<_>, Array2<_>) = (&a).eig().unwrap();
96+
assert_close_l2!(&e,
97+
&arr1(&[c32::new( -9.43,-12.98), c32::new( -3.44, 12.69), c32::new( 0.11, -3.40), c32::new( 5.76, 7.13)]),
98+
1.0e-3);
99+
100+
/*
101+
let answer = &arr2(&[[c32::new( 0.43, 0.33), c32::new( 0.83, 0.00), c32::new( 0.60, 0.00), c32::new( -0.31, 0.03)],
102+
[c32::new( 0.51, -0.03), c32::new( 0.08, -0.25), c32::new( -0.40, -0.20), c32::new( 0.04, 0.34)],
103+
[c32::new( 0.62, 0.00), c32::new( -0.25, 0.28), c32::new( -0.09, -0.48), c32::new( 0.36, 0.06)],
104+
[c32::new( -0.23, 0.11), c32::new( -0.10, -0.32), c32::new( -0.43, 0.13), c32::new( 0.81, 0.00)]]);
105+
*/
106+
107+
for (i, v) in vecs.axis_iter(Axis(1)).enumerate() {
108+
let av = a.dot(&v);
109+
let ev = v.mapv(|f| e[i] * f);
110+
assert_close_l2!(&av, &ev, 1.0e-5);
111+
}
112+
}

0 commit comments

Comments
 (0)