Merge pull request #81 from termoshtt/simplify_solve_mod

Refactoring Factorized{H}

Author: termoshtt

src/solve.rs

@@ -70,7 +70,7 @@ pub use lapack_traits::{Pivot, Transpose};
 ///
 /// If you plan to solve many equations with the same `A` matrix but different
 /// `b` vectors, it's faster to factor the `A` matrix once using the
-/// `Factorize` trait, and then solve using the `Factorized` struct.
+/// `Factorize` trait, and then solve using the `LUFactorized` struct.
 pub trait Solve<A: Scalar> {
     /// Solves a system of linear equations `A * x = b` where `A` is `self`, `b`
     /// is the argument, and `x` is the successful result.
@@ -125,15 +125,15 @@ pub trait Solve<A: Scalar> {
 }
 
 /// Represents the LU factorization of a matrix `A` as `A = P*L*U`.
-pub struct Factorized<S: Data> {
+pub struct LUFactorized<S: Data> {
     /// The factors `L` and `U`; the unit diagonal elements of `L` are not
     /// stored.
     pub a: ArrayBase<S, Ix2>,
     /// The pivot indices that define the permutation matrix `P`.
     pub ipiv: Pivot,
 }
 
-impl<A, S> Solve<A> for Factorized<S>
+impl<A, S> Solve<A> for LUFactorized<S>
 where
     A: Scalar,
     S: Data<Elem = A>,
@@ -213,46 +213,29 @@ where
     }
 }
 
-impl<A, S> Factorized<S>
-where
-    A: Scalar,
-    S: DataMut<Elem = A>,
-{
-    /// Computes the inverse of the factorized matrix.
-    pub fn into_inverse(mut self) -> Result<ArrayBase<S, Ix2>> {
-        unsafe {
-            A::inv(
-                self.a.square_layout()?,
-                self.a.as_allocated_mut()?,
-                &self.ipiv,
-            )?
-        };
-        Ok(self.a)
-    }
-}
 
 /// An interface for computing LU factorizations of matrix refs.
 pub trait Factorize<S: Data> {
     /// Computes the LU factorization `A = P*L*U`, where `P` is a permutation
     /// matrix.
-    fn factorize(&self) -> Result<Factorized<S>>;
+    fn factorize(&self) -> Result<LUFactorized<S>>;
 }
 
 /// An interface for computing LU factorizations of matrices.
 pub trait FactorizeInto<S: Data> {
     /// Computes the LU factorization `A = P*L*U`, where `P` is a permutation
     /// matrix.
-    fn factorize_into(self) -> Result<Factorized<S>>;
+    fn factorize_into(self) -> Result<LUFactorized<S>>;
 }
 
 impl<A, S> FactorizeInto<S> for ArrayBase<S, Ix2>
 where
     A: Scalar,
     S: DataMut<Elem = A>,
 {
-    fn factorize_into(mut self) -> Result<Factorized<S>> {
+    fn factorize_into(mut self) -> Result<LUFactorized<S>> {
         let ipiv = unsafe { A::lu(self.layout()?, self.as_allocated_mut()?)? };
-        Ok(Factorized {
+        Ok(LUFactorized {
             a: self,
             ipiv: ipiv,
         })
@@ -264,10 +247,10 @@ where
     A: Scalar,
     Si: Data<Elem = A>,
 {
-    fn factorize(&self) -> Result<Factorized<OwnedRepr<A>>> {
+    fn factorize(&self) -> Result<LUFactorized<OwnedRepr<A>>> {
         let mut a: Array2<A> = replicate(self);
         let ipiv = unsafe { A::lu(a.layout()?, a.as_allocated_mut()?)? };
-        Ok(Factorized { a: a, ipiv: ipiv })
+        Ok(LUFactorized { a: a, ipiv: ipiv })
     }
 }
 
@@ -285,6 +268,41 @@ pub trait InverseInto {
     fn inv_into(self) -> Result<Self::Output>;
 }
 
+impl<A, S> InverseInto for LUFactorized<S>
+where
+    A: Scalar,
+    S: DataMut<Elem = A>,
+{
+    type Output = ArrayBase<S, Ix2>;
+
+    fn inv_into(mut self) -> Result<ArrayBase<S, Ix2>> {
+        unsafe {
+            A::inv(
+                self.a.square_layout()?,
+                self.a.as_allocated_mut()?,
+                &self.ipiv,
+            )?
+        };
+        Ok(self.a)
+    }
+}
+
+impl<A, S> Inverse for LUFactorized<S>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    type Output = Array2<A>;
+
+    fn inv(&self) -> Result<Array2<A>> {
+        let f = LUFactorized {
+            a: replicate(&self.a),
+            ipiv: self.ipiv.clone(),
+        };
+        f.inv_into()
+    }
+}
+
 impl<A, S> InverseInto for ArrayBase<S, Ix2>
 where
     A: Scalar,
@@ -294,7 +312,7 @@ where
 
     fn inv_into(self) -> Result<Self::Output> {
         let f = self.factorize_into()?;
-        f.into_inverse()
+        f.inv_into()
     }
 }
 
@@ -307,6 +325,6 @@ where
 
     fn inv(&self) -> Result<Self::Output> {
         let f = self.factorize()?;
-        f.into_inverse()
+        f.inv_into()
     }
 }

src/solveh.rs

@@ -63,7 +63,7 @@ pub use lapack_traits::{Pivot, UPLO};
 /// If you plan to solve many equations with the same Hermitian (or real
 /// symmetric) coefficient matrix `A` but different `b` vectors, it's faster to
 /// factor the `A` matrix once using the `FactorizeH` trait, and then solve
-/// using the `FactorizedH` struct.
+/// using the `BKFactorized` struct.
 pub trait SolveH<A: Scalar> {
     /// Solves a system of linear equations `A * x = b` with Hermitian (or real
     /// symmetric) matrix `A`, where `A` is `self`, `b` is the argument, and
@@ -89,12 +89,12 @@ pub trait SolveH<A: Scalar> {
 
 /// Represents the Bunch–Kaufman factorization of a Hermitian (or real
 /// symmetric) matrix as `A = P * U * D * U^H * P^T`.
-pub struct FactorizedH<S: Data> {
+pub struct BKFactorized<S: Data> {
     pub a: ArrayBase<S, Ix2>,
     pub ipiv: Pivot,
 }
 
-impl<A, S> SolveH<A> for FactorizedH<S>
+impl<A, S> SolveH<A> for BKFactorized<S>
 where
     A: Scalar,
     S: Data<Elem = A>,
@@ -131,54 +131,30 @@ where
 }
 
 
-impl<A, S> FactorizedH<S>
-where
-    A: Scalar,
-    S: DataMut<Elem = A>,
-{
-    /// Computes the inverse of the factorized matrix.
-    ///
-    /// **Warning: The inverse is stored only in the upper triangular portion
-    /// of the result matrix!** If you want the lower triangular portion to be
-    /// correct, you must fill it in according to the results in the upper
-    /// triangular portion.
-    pub fn into_inverseh(mut self) -> Result<ArrayBase<S, Ix2>> {
-        unsafe {
-            A::invh(
-                self.a.square_layout()?,
-                UPLO::Upper,
-                self.a.as_allocated_mut()?,
-                &self.ipiv,
-            )?
-        };
-        Ok(self.a)
-    }
-}
-
 /// An interface for computing the Bunch–Kaufman factorization of Hermitian (or
 /// real symmetric) matrix refs.
 pub trait FactorizeH<S: Data> {
     /// Computes the Bunch–Kaufman factorization of a Hermitian (or real
     /// symmetric) matrix.
-    fn factorizeh(&self) -> Result<FactorizedH<S>>;
+    fn factorizeh(&self) -> Result<BKFactorized<S>>;
 }
 
 /// An interface for computing the Bunch–Kaufman factorization of Hermitian (or
 /// real symmetric) matrices.
 pub trait FactorizeHInto<S: Data> {
     /// Computes the Bunch–Kaufman factorization of a Hermitian (or real
     /// symmetric) matrix.
-    fn factorizeh_into(self) -> Result<FactorizedH<S>>;
+    fn factorizeh_into(self) -> Result<BKFactorized<S>>;
 }
 
 impl<A, S> FactorizeHInto<S> for ArrayBase<S, Ix2>
 where
     A: Scalar,
     S: DataMut<Elem = A>,
 {
-    fn factorizeh_into(mut self) -> Result<FactorizedH<S>> {
+    fn factorizeh_into(mut self) -> Result<BKFactorized<S>> {
         let ipiv = unsafe { A::bk(self.layout()?, UPLO::Upper, self.as_allocated_mut()?)? };
-        Ok(FactorizedH {
+        Ok(BKFactorized {
             a: self,
             ipiv: ipiv,
         })
@@ -190,10 +166,10 @@ where
     A: Scalar,
     Si: Data<Elem = A>,
 {
-    fn factorizeh(&self) -> Result<FactorizedH<OwnedRepr<A>>> {
+    fn factorizeh(&self) -> Result<BKFactorized<OwnedRepr<A>>> {
         let mut a: Array2<A> = replicate(self);
         let ipiv = unsafe { A::bk(a.layout()?, UPLO::Upper, a.as_allocated_mut()?)? };
-        Ok(FactorizedH { a: a, ipiv: ipiv })
+        Ok(BKFactorized { a: a, ipiv: ipiv })
     }
 }
 
@@ -221,6 +197,42 @@ pub trait InverseHInto {
     fn invh_into(self) -> Result<Self::Output>;
 }
 
+impl<A, S> InverseHInto for BKFactorized<S>
+where
+    A: Scalar,
+    S: DataMut<Elem = A>,
+{
+    type Output = ArrayBase<S, Ix2>;
+
+    fn invh_into(mut self) -> Result<ArrayBase<S, Ix2>> {
+        unsafe {
+            A::invh(
+                self.a.square_layout()?,
+                UPLO::Upper,
+                self.a.as_allocated_mut()?,
+                &self.ipiv,
+            )?
+        };
+        Ok(self.a)
+    }
+}
+
+impl<A, S> InverseH for BKFactorized<S>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    type Output = Array2<A>;
+
+    fn invh(&self) -> Result<Self::Output> {
+        let f = BKFactorized {
+            a: replicate(&self.a),
+            ipiv: self.ipiv.clone(),
+        };
+        f.invh_into()
+    }
+}
+
 impl<A, S> InverseHInto for ArrayBase<S, Ix2>
 where
     A: Scalar,
@@ -230,7 +242,7 @@ where
 
     fn invh_into(self) -> Result<Self::Output> {
         let f = self.factorizeh_into()?;
-        f.into_inverseh()
+        f.invh_into()
     }
 }
 
@@ -243,6 +255,6 @@ where
 
     fn invh(&self) -> Result<Self::Output> {
         let f = self.factorizeh()?;
-        f.into_inverseh()
+        f.invh_into()
     }
 }