Merge branch 'trianglar'

Author: termoshtt

src/lib.rs

@@ -35,6 +35,7 @@ pub mod vector;
 pub mod matrix;
 pub mod square;
 pub mod hermite;
+pub mod triangular;
 
 pub mod qr;
 pub mod svd;

src/prelude.rs

@@ -3,3 +3,4 @@ pub use vector::Vector;
 pub use matrix::Matrix;
 pub use square::SquareMatrix;
 pub use hermite::HermiteMatrix;
+pub use triangular::TriangularMatrix;

src/solve.rs

@@ -6,28 +6,29 @@ use std::cmp::min;
 use error::LapackError;
 
 pub trait ImplSolve: Sized {
-    fn inv(layout: Layout, size: usize, a: Vec<Self>) -> Result<Vec<Self>, LapackError>;
+    /// execute LU decomposition
     fn lu(layout: Layout, m: usize, n: usize, a: Vec<Self>) -> Result<(Vec<i32>, Vec<Self>), LapackError>;
+    /// calc inverse matrix with LU factorized matrix
+    fn inv(layout: Layout, size: usize, a: Vec<Self>, ipiv: &Vec<i32>) -> Result<Vec<Self>, LapackError>;
+    /// solve linear problem with LU factorized matrix
+    fn solve(layout: Layout,
+             size: usize,
+             a: &Vec<Self>,
+             ipiv: &Vec<i32>,
+             b: Vec<Self>)
+             -> Result<Vec<Self>, LapackError>;
+    /// solve triangular linear problem
+    fn solve_triangle(layout: Layout,
+                      uplo: u8,
+                      size: usize,
+                      a: &[Self],
+                      b: Vec<Self>)
+                      -> Result<Vec<Self>, LapackError>;
 }
 
 macro_rules! impl_solve {
-    ($scalar:ty, $getrf:path, $getri:path, $laswp:path) => {
+    ($scalar:ty, $getrf:path, $getri:path, $getrs:path, $trtrs:path) => {
 impl ImplSolve for $scalar {
-    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 info = $getrf(layout, n, n, &mut a, lda, &mut ipiv);
-        if info != 0 {
-            return Err(From::from(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;
@@ -44,8 +45,42 @@ impl ImplSolve for $scalar {
             Err(From::from(info))
         }
     }
+    fn inv(layout: Layout, size: usize, mut a: Vec<Self>, ipiv: &Vec<i32>) -> Result<Vec<Self>, LapackError> {
+        let n = size as i32;
+        let lda = n;
+        let info = $getri(layout, n, &mut a, lda, &ipiv);
+        if info == 0 {
+            Ok(a)
+        } else {
+            Err(From::from(info))
+        }
+    }
+    fn solve(layout: Layout, size: usize, a: &Vec<Self>, ipiv: &Vec<i32>, mut b: Vec<Self>) -> Result<Vec<Self>, LapackError> {
+        let n = size as i32;
+        let lda = n;
+        let info = $getrs(layout, 'N' as u8, n, 1, a, lda, ipiv, &mut b, n);
+        if info == 0 {
+            Ok(b)
+        } else {
+            Err(From::from(info))
+        }
+    }
+    fn solve_triangle(layout: Layout, uplo: u8, size: usize, a: &[Self], mut b: Vec<Self>) -> Result<Vec<Self>, LapackError> {
+        let n = size as i32;
+        let lda = n;
+        let ldb = match layout {
+            Layout::ColumnMajor => n,
+            Layout::RowMajor => 1,
+        };
+        let info = $trtrs(layout, uplo, 'N' as u8, 'N' as u8, n, 1, a, lda, &mut b, ldb);
+        if info == 0 {
+            Ok(b)
+        } else {
+            Err(From::from(info))
+        }
+    }
 }
 }} // end macro_rules
 
-impl_solve!(f64, dgetrf, dgetri, dlaswp);
-impl_solve!(f32, sgetrf, sgetri, slaswp);
+impl_solve!(f64, dgetrf, dgetri, dgetrs, dtrtrs);
+impl_solve!(f32, sgetrf, sgetri, sgetrs, strtrs);

src/square.rs

@@ -44,7 +44,8 @@ impl<A> SquareMatrix for Array<A, Ix2>
         self.check_square()?;
         let (n, _) = self.size();
         let layout = self.layout()?;
-        let a = ImplSolve::inv(layout, n, self.into_raw_vec())?;
+        let (ipiv, a) = ImplSolve::lu(layout, n, n, self.into_raw_vec())?;
+        let a = ImplSolve::inv(layout, n, a, &ipiv)?;
         let m = Array::from_vec(a).into_shape((n, n)).unwrap();
         match layout {
             Layout::RowMajor => Ok(m),

src/triangular.rs

@@ -0,0 +1,40 @@
+
+use ndarray::{Ix2, Array, LinalgScalar};
+use std::fmt::Debug;
+use num_traits::float::Float;
+
+use matrix::Matrix;
+use square::SquareMatrix;
+use error::LinalgError;
+use qr::ImplQR;
+use svd::ImplSVD;
+use norm::ImplNorm;
+use solve::ImplSolve;
+
+pub trait TriangularMatrix: Matrix + SquareMatrix {
+    /// solve a triangular system with upper triangular matrix
+    fn solve_upper(&self, Self::Vector) -> Result<Self::Vector, LinalgError>;
+    /// solve a triangular system with lower triangular matrix
+    fn solve_lower(&self, Self::Vector) -> Result<Self::Vector, LinalgError>;
+}
+
+impl<A> TriangularMatrix for Array<A, Ix2>
+    where A: ImplQR + ImplNorm + ImplSVD + ImplSolve + LinalgScalar + Float + Debug
+{
+    fn solve_upper(&self, b: Self::Vector) -> Result<Self::Vector, LinalgError> {
+        self.check_square()?;
+        let (n, _) = self.size();
+        let layout = self.layout()?;
+        let a = self.as_slice_memory_order().unwrap();
+        let x = ImplSolve::solve_triangle(layout, 'U' as u8, n, a, b.into_raw_vec())?;
+        Ok(Array::from_vec(x))
+    }
+    fn solve_lower(&self, b: Self::Vector) -> Result<Self::Vector, LinalgError> {
+        self.check_square()?;
+        let (n, _) = self.size();
+        let layout = self.layout()?;
+        let a = self.as_slice_memory_order().unwrap();
+        let x = ImplSolve::solve_triangle(layout, 'L' as u8, n, a, b.into_raw_vec())?;
+        Ok(Array::from_vec(x))
+    }
+}

tests/triangular.rs

@@ -0,0 +1,90 @@
+
+extern crate rand;
+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;
+
+fn all_close<D: Dimension>(a: Array<f64, D>, b: Array<f64, D>) {
+    if !a.all_close(&b, 1.0e-7) {
+        panic!("\nTwo matrices are not equal:\na = \n{:?}\nb = \n{:?}\n",
+               a,
+               b);
+    }
+}
+
+#[test]
+fn solve_upper() {
+    let r_dist = Range::new(0., 1.);
+    let mut a = Array::<f64, _>::random((3, 3), r_dist);
+    for ((i, j), val) in a.indexed_iter_mut() {
+        if i > j {
+            *val = 0.0;
+        }
+    }
+    println!("a = \n{:?}", &a);
+    let b = Array::<f64, _>::random(3, r_dist);
+    println!("b = \n{:?}", &b);
+    let x = a.solve_upper(b.clone()).unwrap();
+    println!("x = \n{:?}", &x);
+    println!("Ax = \n{:?}", a.dot(&x));
+    all_close(a.dot(&x), b);
+}
+
+#[test]
+fn solve_upper_t() {
+    let r_dist = Range::new(0., 1.);
+    let mut a = Array::<f64, _>::random((3, 3), r_dist).reversed_axes();
+    for ((i, j), val) in a.indexed_iter_mut() {
+        if i > j {
+            *val = 0.0;
+        }
+    }
+    println!("a = \n{:?}", &a);
+    let b = Array::<f64, _>::random(3, r_dist);
+    println!("b = \n{:?}", &b);
+    let x = a.solve_upper(b.clone()).unwrap();
+    println!("x = \n{:?}", &x);
+    println!("Ax = \n{:?}", a.dot(&x));
+    all_close(a.dot(&x), b);
+}
+
+#[test]
+fn solve_lower() {
+    let r_dist = Range::new(0., 1.);
+    let mut a = Array::<f64, _>::random((3, 3), r_dist);
+    for ((i, j), val) in a.indexed_iter_mut() {
+        if i < j {
+            *val = 0.0;
+        }
+    }
+    println!("a = \n{:?}", &a);
+    let b = Array::<f64, _>::random(3, r_dist);
+    println!("b = \n{:?}", &b);
+    let x = a.solve_lower(b.clone()).unwrap();
+    println!("x = \n{:?}", &x);
+    println!("Ax = \n{:?}", a.dot(&x));
+    all_close(a.dot(&x), b);
+}
+
+#[test]
+fn solve_lower_t() {
+    let r_dist = Range::new(0., 1.);
+    let mut a = Array::<f64, _>::random((3, 3), r_dist).reversed_axes();
+    for ((i, j), val) in a.indexed_iter_mut() {
+        if i < j {
+            *val = 0.0;
+        }
+    }
+    println!("a = \n{:?}", &a);
+    let b = Array::<f64, _>::random(3, r_dist);
+    println!("b = \n{:?}", &b);
+    let x = a.solve_lower(b.clone()).unwrap();
+    println!("x = \n{:?}", &x);
+    println!("Ax = \n{:?}", a.dot(&x));
+    all_close(a.dot(&x), b);
+}