Implement solve_lower

Author: termoshtt

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>
@@ -29,4 +29,12 @@ impl<A> TriangularMatrix for Array<A, Ix2>
         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

@@ -52,3 +52,39 @@ fn solve_upper_t() {
     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);
+}