Eigenvalue analysis for general matrix

Cargo.toml

@@ -12,6 +12,7 @@ license = "MIT"
 [dependencies]
 lapack = "0.11.1"
 num-traits = "0.1.36"
+num-complex = "*"
 
 [dependencies.ndarray]
 version = "0.6.9"

src/eig.rs

@@ -0,0 +1,44 @@
+//! Implement eigenvalue decomposition of general matrix
+
+use lapack::fortran::*;
+use num_traits::Zero;
+
+use error::LapackError;
+
+pub trait ImplEig: Sized {
+    fn eig(n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>, Vec<Self>), LapackError>;
+}
+
+macro_rules! impl_eig {
+    ($scalar:ty, $geev:path) => {
+impl ImplEig for $scalar {
+    fn eig(n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>, Vec<Self>), LapackError> {
+        let lda = n as i32;
+        let ldvr = n as i32;
+        let ldvl = n as i32;
+        let lw_default = 1000;
+        let mut wr = vec![Self::zero(); n];
+        let mut wi = vec![Self::zero(); n];
+        let mut vr = vec![Self::zero(); n * n];
+        let mut vl = Vec::new();
+        let mut work = vec![Self::zero(); lw_default];
+        let mut info = 0;
+        $geev(b'N', b'V', n as i32, &mut a, lda, &mut wr, &mut wi,
+              &mut vl, ldvl, &mut vr, ldvr, &mut work, -1, &mut info);
+        let lwork = work[0] as i32;
+        if lwork > lw_default as i32 {
+            work = vec![Self::zero(); lwork as usize];
+        }
+        $geev(b'N', b'V', n as i32, &mut a, lda, &mut wr, &mut wi,
+              &mut vl, ldvl, &mut vr, ldvr, &mut work, lwork, &mut info);
+        if info == 0 {
+            Ok((wr, wi, vr))
+        } else {
+            Err(From::from(info))
+        }
+    }
+}
+}} // end macro_rules
+
+impl_eig!(f64, dgeev);
+impl_eig!(f32, sgeev);

src/hermite.rs

@@ -9,6 +9,7 @@ use lapack::c::Layout;
 use matrix::Matrix;
 use square::SquareMatrix;
 use error::LinalgError;
+use eig::ImplEig;
 use eigh::ImplEigh;
 use qr::ImplQR;
 use svd::ImplSVD;
@@ -27,7 +28,7 @@ pub trait HermiteMatrix: SquareMatrix + Matrix {
 }
 
 impl<A> HermiteMatrix for Array<A, (Ix, Ix)>
-    where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + ImplEigh + ImplCholesky + LinalgScalar + Float + Debug
+    where A: ImplEig + ImplQR + ImplSVD + ImplNorm + ImplSolve + ImplEigh + LinalgScalar + Float + Debug
 {
     fn eigh(self) -> Result<(Self::Vector, Self), LinalgError> {
         try!(self.check_square());

src/lib.rs

@@ -24,10 +24,11 @@
 //! - [symmetric square root](hermite/trait.HermiteMatrix.html#tymethod.ssqrt)
 //! - [WIP] Cholesky factorization
 
-extern crate lapack;
-extern crate num_traits;
 #[macro_use(s)]
 extern crate ndarray;
+extern crate lapack;
+extern crate num_complex;
+extern crate num_traits;
 
 pub mod prelude;
 pub mod error;
@@ -38,6 +39,7 @@ pub mod hermite;
 
 pub mod qr;
 pub mod svd;
+pub mod eig;
 pub mod eigh;
 pub mod norm;
 pub mod solve;

src/matrix.rs

@@ -2,6 +2,7 @@
 
 use std::cmp::min;
 use std::fmt::Debug;
+use num_complex::Complex;
 use ndarray::prelude::*;
 use ndarray::LinalgScalar;
 use lapack::c::Layout;
@@ -14,19 +15,34 @@ use solve::ImplSolve;
 
 /// Methods for general matrices
 pub trait Matrix: Sized {
+    /// Corresponding real type
+    type Real;
+    /// Corresponding complex type
+    type Complex;
+    /// Scalar type
     type Scalar;
+    /// Corresponding one-dimensional vector
     type Vector;
     type Permutator;
+    /// Corresponding real one-dimensional vector
+    type RealVector;
+    /// Corresponding complex one-dimensional vector
+    type ComplexVector;
+    /// Corresponding real matrix
+    type RealMatrix;
+    /// Corresponding complex matrix
+    type ComplexMatrix;
+
     /// number of (rows, columns)
     fn size(&self) -> (usize, usize);
     /// Layout (C/Fortran) of matrix
     fn layout(&self) -> Layout;
     /// Operator norm for L-1 norm
-    fn norm_1(&self) -> Self::Scalar;
+    fn norm_1(&self) -> Self::Real;
     /// Operator norm for L-inf norm
-    fn norm_i(&self) -> Self::Scalar;
+    fn norm_i(&self) -> Self::Real;
     /// Frobenius norm
-    fn norm_f(&self) -> Self::Scalar;
+    fn norm_f(&self) -> Self::Real;
     /// singular-value decomposition (SVD)
     fn svd(self) -> Result<(Self, Self::Vector, Self), LapackError>;
     /// QR decomposition
@@ -45,10 +61,18 @@ pub trait Matrix: Sized {
 impl<A> Matrix for Array<A, (Ix, Ix)>
     where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + LinalgScalar + Debug
 {
+    type Real = A;
+    type Complex = Complex<A>;
     type Scalar = A;
+
+    type RealVector = Array<A, Ix>;
+    type ComplexVector = Array<Self::Complex, Ix>;
     type Vector = Array<A, Ix>;
     type Permutator = Vec<i32>;
 
+    type RealMatrix = Array<A, (Ix, Ix)>;
+    type ComplexMatrix = Array<Self::Complex, (Ix, Ix)>;
+
     fn size(&self) -> (usize, usize) {
         (self.rows(), self.cols())
     }
@@ -60,7 +84,7 @@ impl<A> Matrix for Array<A, (Ix, Ix)>
             Layout::RowMajor
         }
     }
-    fn norm_1(&self) -> Self::Scalar {
+    fn norm_1(&self) -> Self::Real {
         let (m, n) = self.size();
         let strides = self.strides();
         if strides[0] > strides[1] {
@@ -69,7 +93,7 @@ impl<A> Matrix for Array<A, (Ix, Ix)>
             ImplNorm::norm_1(m, n, self.clone().into_raw_vec())
         }
     }
-    fn norm_i(&self) -> Self::Scalar {
+    fn norm_i(&self) -> Self::Real {
         let (m, n) = self.size();
         let strides = self.strides();
         if strides[0] > strides[1] {
@@ -78,7 +102,7 @@ impl<A> Matrix for Array<A, (Ix, Ix)>
             ImplNorm::norm_i(m, n, self.clone().into_raw_vec())
         }
     }
-    fn norm_f(&self) -> Self::Scalar {
+    fn norm_f(&self) -> Self::Real {
         let (m, n) = self.size();
         ImplNorm::norm_f(m, n, self.clone().into_raw_vec())
     }

src/square.rs

@@ -4,10 +4,12 @@ use std::fmt::Debug;
 use ndarray::prelude::*;
 use ndarray::LinalgScalar;
 use num_traits::float::Float;
+use num_complex::Complex;
 
 use matrix::Matrix;
 use error::{LinalgError, NotSquareError};
 use qr::ImplQR;
+use eig::ImplEig;
 use svd::ImplSVD;
 use norm::ImplNorm;
 use solve::ImplSolve;
@@ -18,7 +20,7 @@ use solve::ImplSolve;
 /// but does not assure that the matrix is square.
 /// If not square, `NotSquareError` will be thrown.
 pub trait SquareMatrix: Matrix {
-    // fn eig(self) -> (Self::Vector, Self);
+    fn eig(self) -> Result<(Self::ComplexVector, Self::ComplexMatrix), LinalgError>;
     /// inverse matrix
     fn inv(self) -> Result<Self, LinalgError>;
     /// trace of matrix
@@ -38,8 +40,40 @@ pub trait SquareMatrix: Matrix {
 }
 
 impl<A> SquareMatrix for Array<A, (Ix, Ix)>
-    where A: ImplQR + ImplNorm + ImplSVD + ImplSolve + LinalgScalar + Float + Debug
+    where A: ImplEig + ImplQR + ImplNorm + ImplSVD + ImplSolve + LinalgScalar + Float + Debug
 {
+    fn eig(self) -> Result<(Self::ComplexVector, Self::ComplexMatrix), LinalgError> {
+        try!(self.check_square());
+        let (n, _) = self.size();
+        let (wr, wi, vv) = try!(ImplEig::eig(n, self.into_raw_vec()));
+        println!("wi = {:?}", &wi);
+        let vr = Array::from_vec(vv).into_shape((n, n)).unwrap();
+        let mut v = Array::<Self::Complex, _>::zeros((n, n));
+        let mut i = 0;
+        while i < n {
+            println!("i = {}", &i);
+            println!("wi[i] = {:?}", &wi[i]);
+            if !wi[i].is_normal() {
+                println!("Real eigenvalue");
+                for j in 0..n {
+                    v[(i, j)] = Complex::new(vr[(i, j)], A::zero());
+                }
+                i += 1;
+            } else {
+                println!("Imaginal eigenvalue");
+                for j in 0..n {
+                    v[(i, j)] = Complex::new(vr[(i, j)], vr[(i + 1, j)]);
+                    v[(i + 1, j)] = Complex::new(vr[(i, j)], -vr[(i + 1, j)]);
+                }
+                i += 2;
+            }
+        }
+        let w = wr.into_iter()
+            .zip(wi.into_iter())
+            .map(|(r, i)| Complex::new(r, i))
+            .collect();
+        Ok((w, v.reversed_axes()))
+    }
     fn inv(self) -> Result<Self, LinalgError> {
         try!(self.check_square());
         let (n, _) = self.size();

tests/eig.rs

@@ -0,0 +1,53 @@
+
+extern crate rand;
+extern crate num_complex;
+extern crate ndarray;
+extern crate ndarray_rand;
+extern crate ndarray_linalg;
+
+use ndarray::prelude::*;
+use ndarray_linalg::prelude::*;
+use rand::distributions::*;
+use ndarray_rand::RandomExt;
+use num_complex::Complex;
+
+type c64 = Complex<f64>;
+
+fn dot(a: Array<c64, (Ix, Ix)>, b: Array<f64, (Ix, Ix)>) {
+    let (n, l) = a.size();
+    let (_, m) = b.size();
+    let mut c = Array::<c64, _>::zeros((n, m));
+    for i in 0..n {
+        for j in 0..m {
+            for k in 0..l {
+                c[(i, j)] = c[(i, j)] + a[(i, k)] * b[(k, j)];
+            }
+        }
+    }
+    c
+}
+
+#[test]
+fn eig_random() {
+    let r_dist = Range::new(0., 1.);
+    let a = Array::<f64, _>::random((3, 3), r_dist);
+    println!("a = \n{:?}", &a);
+    let (w, vr) = a.clone().eig().unwrap();
+    println!("w = \n{:?}", &w);
+    println!("vr = \n{:?}", &vr);
+    let mut lm = Array::zeros((3, 3));
+    for i in 0..3 {
+        lm[(i, i)] = w[i];
+    }
+    println!("lm = \n{:?}", &lm);
+    let mut lv = Array::<Complex<f64>, _>::zeros((3, 3));
+    for i in 0..3 {
+        for j in 0..3 {
+            for k in 0..3 {
+                lv[(i, j)] = lv[(i, j)] + lm[(i, k)] * vr[(k, j)];
+            }
+        }
+    }
+    println!("lv = \n{:?}", &lv);
+    panic!("Manual Kill");
+}