Add ImplSolve::solve to reuse LU factorized matrix

Author: termoshtt

src/solve.rs

@@ -6,28 +6,22 @@ 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>;
 }
 
 macro_rules! impl_solve {
-    ($scalar:ty, $getrf:path, $getri:path, $laswp:path) => {
+    ($scalar:ty, $getrf:path, $getri:path, $getrs: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 +38,28 @@ 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))
+        }
+    }
 }
 }} // end macro_rules
 
-impl_solve!(f64, dgetrf, dgetri, dlaswp);
-impl_solve!(f32, sgetrf, sgetri, slaswp);
+impl_solve!(f64, dgetrf, dgetri, dgetrs);
+impl_solve!(f32, sgetrf, sgetri, sgetrs);

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),