Add test for solve_upper

Author: termoshtt

src/solve.rs

@@ -21,7 +21,7 @@ pub trait ImplSolve: Sized {
     fn solve_triangle(layout: Layout,
                       uplo: u8,
                       size: usize,
-                      a: &Vec<Self>,
+                      a: &[Self],
                       b: Vec<Self>)
                       -> Result<Vec<Self>, LapackError>;
 }
@@ -65,10 +65,14 @@ impl ImplSolve for $scalar {
             Err(From::from(info))
         }
     }
-    fn solve_triangle(layout: Layout, uplo: u8, size: usize, a: &Vec<Self>, mut b: Vec<Self>) -> Result<Vec<Self>, LapackError> {
+    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 info = $trtrs(layout, uplo, 'N' as u8, 'N' as u8, n, 1, a, lda, &mut b, 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 {

src/triangular.rs

@@ -14,8 +14,8 @@ 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>;
+    // /// 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>
@@ -25,12 +25,8 @@ impl<A> TriangularMatrix for Array<A, Ix2>
         self.check_square()?;
         let (n, _) = self.size();
         let layout = self.layout()?;
-        let x = ImplSolve::solve_triangle(layout, 'U' as u8, n, self.as_slice().unwrap(), b)?;
-    }
-    fn solve_lower(&self, b: Self::Vector) -> Result<Self::Vector, LinalgError> {
-        self.check_square()?;
-        let (n, _) = self.size();
-        let layout = self.layout()?;
-        let x = ImplSolve::solve_triangle(layout, 'U' as u8, n, self.as_slice().unwrap(), b)?;
+        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))
     }
 }

tests/triangular.rs

@@ -0,0 +1,54 @@
+
+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);
+}