Merge branch 'master' into index_type

Conflicts:
	src/hermite.rs
	src/matrix.rs
	src/square.rs

Author: termoshtt

.travis.yml

@@ -0,0 +1,28 @@
+language: rust
+sudo: required
+
+rust:
+  - stable
+  - beta
+  - nightly
+
+os:
+  - linux
+
+cache: cargo
+
+matrix:
+  allow-failures:
+    - rust: nightly
+
+script:
+  - cargo test -vv --no-default-features --features=$FEATURE
+
+addons:
+  apt:
+    sources:
+      - kubuntu-backports
+
+    packages:
+      - cmake
+      - gfortran

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "ndarray-linalg"
-version = "0.1.2"
+version = "0.1.3"
 authors = ["Toshiki Teramura <[email protected]>"]
 
 description = "Linear algebra package for rust-ndarray using LAPACK"

README.md

@@ -1,4 +1,4 @@
-ndarray-linalg [![Crate](http://meritbadge.herokuapp.com/ndarray-linalg)](https://crates.io/crates/ndarray-linalg) [![docs.rs](https://docs.rs/ndarray-linalg/badge.svg)](https://docs.rs/ndarray-linalg) [![wercker status](https://app.wercker.com/status/a45df26fa97eab7debf53b32fc576b35/s/master "wercker status")](https://app.wercker.com/project/byKey/a45df26fa97eab7debf53b32fc576b35)
+ndarray-linalg [![Crate](http://meritbadge.herokuapp.com/ndarray-linalg)](https://crates.io/crates/ndarray-linalg) [![docs.rs](https://docs.rs/ndarray-linalg/badge.svg)](https://docs.rs/ndarray-linalg) [![Build Status](https://travis-ci.org/termoshtt/ndarray-linalg.svg?branch=master)](https://travis-ci.org/termoshtt/ndarray-linalg)
 ===============
 Linear algebra package for [rust-ndarray](https://github.com/bluss/rust-ndarray) using LAPACK via [stainless-steel/lapack](https://github.com/stainless-steel/lapack)
 

src/cholesky.rs

@@ -0,0 +1,24 @@
+
+use lapack::c::*;
+use error::LapackError;
+
+pub trait ImplCholesky: Sized {
+    fn cholesky(layout: Layout, n: usize, a: Vec<Self>) -> Result<Vec<Self>, LapackError>;
+}
+
+macro_rules! impl_cholesky {
+    ($scalar:ty, $potrf:path) => {
+impl ImplCholesky for $scalar {
+    fn cholesky(layout: Layout, n: usize, mut a: Vec<Self>) -> Result<Vec<Self>, LapackError> {
+        let info = $potrf(layout, b'U', n as i32, &mut a, n as i32);
+        if info == 0 {
+            Ok(a)
+        } else {
+            Err(From::from(info))
+        }
+    }
+}
+}} // end macro_rules
+
+impl_cholesky!(f64, dpotrf);
+impl_cholesky!(f32, spotrf);

src/hermite.rs

@@ -1,7 +1,9 @@
 //! Define trait for Hermite matrices
 
 use ndarray::{Ix2, Array, LinalgScalar};
+use std::fmt::Debug;
 use num_traits::float::Float;
+use lapack::c::Layout;
 
 use matrix::Matrix;
 use square::SquareMatrix;
@@ -11,17 +13,20 @@ use qr::ImplQR;
 use svd::ImplSVD;
 use norm::ImplNorm;
 use solve::ImplSolve;
+use cholesky::ImplCholesky;
 
 /// Methods for Hermite matrix
 pub trait HermiteMatrix: SquareMatrix + Matrix {
     /// eigenvalue decomposition
     fn eigh(self) -> Result<(Self::Vector, Self), LinalgError>;
     /// symmetric square root of Hermite matrix
     fn ssqrt(self) -> Result<Self, LinalgError>;
+    /// Cholesky factorization
+    fn cholesky(self) -> Result<Self, LinalgError>;
 }
 
 impl<A> HermiteMatrix for Array<A, Ix2>
-    where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + ImplEigh + LinalgScalar + Float
+    where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + ImplEigh + ImplCholesky + LinalgScalar + Float + Debug
 {
     fn eigh(self) -> Result<(Self::Vector, Self), LinalgError> {
         try!(self.check_square());
@@ -42,4 +47,21 @@ impl<A> HermiteMatrix for Array<A, Ix2>
         }
         Ok(v.dot(&res))
     }
+    fn cholesky(self) -> Result<Self, LinalgError> {
+        try!(self.check_square());
+        println!("layout = {:?}", self.layout());
+        let (n, _) = self.size();
+        let layout = self.layout();
+        let a = try!(ImplCholesky::cholesky(layout, n, self.into_raw_vec()));
+        let mut c = match layout {
+          Layout::RowMajor => Array::from_vec(a).into_shape((n, n)).unwrap(),
+          Layout::ColumnMajor => Array::from_vec(a).into_shape((n, n)).unwrap().reversed_axes(),
+        };
+        for ((i, j), val) in c.indexed_iter_mut() {
+            if i > j {
+                *val = A::zero();
+            }
+        }
+        Ok(c)
+    }
 }

src/lib.rs

@@ -10,23 +10,24 @@
 //! - [operator norm for L1 norm](matrix/trait.Matrix.html#tymethod.norm_1)
 //! - [operator norm for L-inf norm](matrix/trait.Matrix.html#tymethod.norm_i)
 //! - [Frobeiuns norm](matrix/trait.Matrix.html#tymethod.norm_f)
+//! - [LU factorization](matrix/trait.Matrix.html#tymethod.lu)
 //!
 //! SquareMatrix
 //! -------------
 //! - [inverse of matrix](square/trait.SquareMatrix.html#tymethod.inv)
 //! - [trace of matrix](square/trait.SquareMatrix.html#tymethod.trace)
 //! - [WIP] eigenvalue
-//! - [WIP] LU factorization
 //!
 //! HermiteMatrix
 //! --------------
 //! - [eigenvalue analysis](hermite/trait.HermiteMatrix.html#tymethod.eigh)
 //! - [symmetric square root](hermite/trait.HermiteMatrix.html#tymethod.ssqrt)
-//! - [WIP] Cholesky factorization
+//! - [Cholesky factorization](hermite/trait.HermiteMatrix.html#tymethod.cholesky)
 
 extern crate lapack;
-extern crate ndarray;
 extern crate num_traits;
+#[macro_use(s)]
+extern crate ndarray;
 
 pub mod prelude;
 pub mod error;
@@ -40,3 +41,4 @@ pub mod svd;
 pub mod eigh;
 pub mod norm;
 pub mod solve;
+pub mod cholesky;

src/matrix.rs

@@ -1,19 +1,26 @@
 //! Define trait for general matrix
 
 use std::cmp::min;
-use ndarray::{Ix1, Ix2, Array, Axis, LinalgScalar};
+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
@@ -24,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, Ix2>
-    where A: ImplQR + ImplSVD + ImplNorm + LinalgScalar
+    where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + LinalgScalar + Debug
 {
     type Scalar = A;
-    type Vector = Array<A, Ix1>;
+    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();
@@ -111,4 +137,54 @@ impl<A> Matrix for Array<A, Ix2>
         }
         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,6 +1,7 @@
 //! Define trait for Hermite matrices
 
 use ndarray::{Ix2, Array, LinalgScalar};
+use std::fmt::Debug;
 use num_traits::float::Float;
 
 use matrix::Matrix;
@@ -16,7 +17,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>;
@@ -37,13 +37,13 @@ pub trait SquareMatrix: Matrix {
 }
 
 impl<A> SquareMatrix for Array<A, Ix2>
-    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)