Merge pull request #18 from termoshtt/lu

LU decomposition

Author: termoshtt

src/hermite.rs

@@ -1,5 +1,6 @@
 //! Define trait for Hermite matrices
 
+use std::fmt::Debug;
 use ndarray::prelude::*;
 use ndarray::LinalgScalar;
 use num_traits::float::Float;
@@ -22,7 +23,7 @@ pub trait HermiteMatrix: SquareMatrix + Matrix {
 }
 
 impl<A> HermiteMatrix for Array<A, (Ix, Ix)>
-    where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + ImplEigh + LinalgScalar + Float
+    where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + ImplEigh + LinalgScalar + Float + Debug
 {
     fn eigh(self) -> Result<(Self::Vector, Self), LinalgError> {
         try!(self.check_square());

src/lib.rs

@@ -25,8 +25,9 @@
 //! - [WIP] Cholesky factorization
 
 extern crate lapack;
-extern crate ndarray;
 extern crate num_traits;
+#[macro_use(s)]
+extern crate ndarray;
 
 pub mod prelude;
 pub mod error;

src/matrix.rs

@@ -1,20 +1,26 @@
 //! Define trait for general matrix
 
 use std::cmp::min;
+use std::fmt::Debug;
 use ndarray::prelude::*;
 use ndarray::LinalgScalar;
+use lapack::c::Layout;
 
 use error::LapackError;
 use qr::ImplQR;
 use svd::ImplSVD;
 use norm::ImplNorm;
+use solve::ImplSolve;
 
 /// Methods for general matrices
 pub trait Matrix: Sized {
     type Scalar;
     type Vector;
+    type Permutator;
     /// 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;
     /// Operator norm for L-inf norm
@@ -25,16 +31,35 @@ pub trait Matrix: Sized {
     fn svd(self) -> Result<(Self, Self::Vector, Self), LapackError>;
     /// QR decomposition
     fn qr(self) -> Result<(Self, Self), LapackError>;
+    /// LU decomposition
+    fn lu(self) -> Result<(Self::Permutator, Self, Self), LapackError>;
+    /// permutate matrix (inplace)
+    fn permutate(&mut self, p: &Self::Permutator);
+    /// permutate matrix (outplace)
+    fn permutated(mut self, p: &Self::Permutator) -> Self {
+        self.permutate(p);
+        self
+    }
 }
 
 impl<A> Matrix for Array<A, (Ix, Ix)>
-    where A: ImplQR + ImplSVD + ImplNorm + LinalgScalar
+    where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + LinalgScalar + Debug
 {
     type Scalar = A;
     type Vector = Array<A, Ix>;
+    type Permutator = Vec<i32>;
+
     fn size(&self) -> (usize, usize) {
         (self.rows(), self.cols())
     }
+    fn layout(&self) -> Layout {
+        let strides = self.strides();
+        if strides[0] < strides[1] {
+            Layout::ColumnMajor
+        } else {
+            Layout::RowMajor
+        }
+    }
     fn norm_1(&self) -> Self::Scalar {
         let (m, n) = self.size();
         let strides = self.strides();
@@ -112,4 +137,54 @@ impl<A> Matrix for Array<A, (Ix, Ix)>
         }
         Ok((qm, rm))
     }
+    fn lu(self) -> Result<(Self::Permutator, Self, Self), LapackError> {
+        let (n, m) = self.size();
+        println!("n={}, m={}", n, m);
+        let k = min(n, m);
+        let (p, mut a) = match self.layout() {
+            Layout::ColumnMajor => {
+                println!("ColumnMajor");
+                let (p, l) = ImplSolve::lu(self.layout(), n, m, self.clone().into_raw_vec())?;
+                (p, Array::from_vec(l).into_shape((m, n)).unwrap().reversed_axes())
+            }
+            Layout::RowMajor => {
+                println!("RowMajor");
+                let (p, l) = ImplSolve::lu(self.layout(), n, m, self.clone().into_raw_vec())?;
+                (p, Array::from_vec(l).into_shape((n, m)).unwrap())
+            }
+        };
+        println!("a (after LU) = \n{:?}", &a);
+        let mut lm = Array::zeros((n, k));
+        for ((i, j), val) in lm.indexed_iter_mut() {
+            if i > j {
+                *val = a[(i, j)];
+            } else if i == j {
+                *val = A::one();
+            }
+        }
+        for ((i, j), val) in a.indexed_iter_mut() {
+            if i > j {
+                *val = A::zero();
+            }
+        }
+        let am = if n > k {
+            a.slice(s![0..k as isize, ..]).to_owned()
+        } else {
+            a
+        };
+        println!("am = \n{:?}", am);
+        Ok((p, lm, am))
+    }
+    fn permutate(&mut self, ipiv: &Self::Permutator) {
+        let (_, m) = self.size();
+        for (i, j_) in ipiv.iter().enumerate().rev() {
+            let j = (j_ - 1) as usize;
+            if i == j {
+                continue;
+            }
+            for k in 0..m {
+                self.swap((i, k), (j, k));
+            }
+        }
+    }
 }

src/solve.rs

@@ -1,37 +1,55 @@
 //! Implement linear solver and inverse matrix
 
-use lapack::fortran::*;
-use num_traits::Zero;
+use lapack::c::*;
+use std::cmp::min;
 
 use error::LapackError;
 
 pub trait ImplSolve: Sized {
-    fn inv(size: usize, mut a: Vec<Self>) -> Result<Vec<Self>, LapackError>;
+    fn inv(layout: Layout, size: usize, a: Vec<Self>) -> Result<Vec<Self>, LapackError>;
+    fn lu(layout: Layout,
+          m: usize,
+          n: usize,
+          a: Vec<Self>)
+          -> Result<(Vec<i32>, Vec<Self>), LapackError>;
 }
 
 macro_rules! impl_solve {
-    ($scalar:ty, $getrf:path, $getri:path) => {
+    ($scalar:ty, $getrf:path, $getri:path, $laswp:path) => {
 impl ImplSolve for $scalar {
-    fn inv(size: usize, mut a: Vec<Self>) -> Result<Vec<Self>, LapackError> {
+    fn inv(layout: Layout, size: usize, mut a: Vec<Self>) -> Result<Vec<Self>, LapackError> {
         let n = size as i32;
         let lda = n;
         let mut ipiv = vec![0; size];
-        let mut info = 0;
-        $getrf(n, n, &mut a, lda, &mut ipiv, &mut info);
+        let info = $getrf(layout, n, n, &mut a, lda, &mut ipiv);
         if info != 0 {
             return Err(From::from(info));
         }
-        let lwork = n;
-        let mut work = vec![Self::zero(); size];
-        $getri(n, &mut a, lda, &mut ipiv, &mut work, lwork, &mut info);
+        let info = $getri(layout, n, &mut a, lda, &mut ipiv);
         if info == 0 {
             Ok(a)
         } else {
             Err(From::from(info))
         }
     }
+    fn lu(layout: Layout, m: usize, n: usize, mut a: Vec<Self>) -> Result<(Vec<i32>, Vec<Self>), LapackError> {
+        let m = m as i32;
+        let n = n as i32;
+        let k = min(m, n);
+        let lda = match layout {
+            Layout::ColumnMajor => m,
+            Layout::RowMajor => n,
+        };
+        let mut ipiv = vec![0; k as usize];
+        let info = $getrf(layout, m, n, &mut a, lda, &mut ipiv);
+        if info == 0 {
+            Ok((ipiv, a))
+        } else {
+            Err(From::from(info))
+        }
+    }
 }
 }} // end macro_rules
 
-impl_solve!(f64, dgetrf, dgetri);
-impl_solve!(f32, sgetrf, sgetri);
+impl_solve!(f64, dgetrf, dgetri, dlaswp);
+impl_solve!(f32, sgetrf, sgetri, slaswp);

src/square.rs

@@ -1,5 +1,6 @@
 //! Define trait for Hermite matrices
 
+use std::fmt::Debug;
 use ndarray::prelude::*;
 use ndarray::LinalgScalar;
 use num_traits::float::Float;
@@ -17,7 +18,6 @@ use solve::ImplSolve;
 /// but does not assure that the matrix is square.
 /// If not square, `NotSquareError` will be thrown.
 pub trait SquareMatrix: Matrix {
-    // fn lu(self) -> (Self, Self);
     // fn eig(self) -> (Self::Vector, Self);
     /// inverse matrix
     fn inv(self) -> Result<Self, LinalgError>;
@@ -38,13 +38,13 @@ pub trait SquareMatrix: Matrix {
 }
 
 impl<A> SquareMatrix for Array<A, (Ix, Ix)>
-    where A: ImplQR + ImplNorm + ImplSVD + ImplSolve + LinalgScalar + Float
+    where A: ImplQR + ImplNorm + ImplSVD + ImplSolve + LinalgScalar + Float + Debug
 {
     fn inv(self) -> Result<Self, LinalgError> {
         try!(self.check_square());
         let (n, _) = self.size();
         let is_fortran_align = self.strides()[0] > self.strides()[1];
-        let a = try!(ImplSolve::inv(n, self.into_raw_vec()));
+        let a = try!(ImplSolve::inv(self.layout(), n, self.into_raw_vec()));
         let m = Array::from_vec(a).into_shape((n, n)).unwrap();
         if is_fortran_align {
             Ok(m)