Add solveh module (without test)

Author: termoshtt

src/lapack_traits/mod.rs

@@ -24,7 +24,8 @@ use super::types::*;
 pub type Pivot = Vec<i32>;
 
 pub trait LapackScalar
-    : OperatorNorm_ + QR_ + SVD_ + Solve_ + Cholesky_ + Eigh_ + Triangular_ {
+    : OperatorNorm_ + QR_ + SVD_ + Solve_ + Solveh_ + Cholesky_ + Eigh_ + Triangular_
+    {
 }
 
 impl LapackScalar for f32 {}

src/lapack_traits/solveh.rs

@@ -14,9 +14,9 @@ pub trait Solveh_: Sized {
     /// Bunch-Kaufman: wrapper of `*sytrf` and `*hetrf`
     unsafe fn bk(MatrixLayout, UPLO, a: &mut [Self]) -> Result<Pivot>;
     /// Wrapper of `*sytri` and `*hetri`
-    unsafe fn inv(MatrixLayout, UPLO, a: &mut [Self], &Pivot) -> Result<()>;
+    unsafe fn invh(MatrixLayout, UPLO, a: &mut [Self], &Pivot) -> Result<()>;
     /// Wrapper of `*sytrs` and `*hetrs`
-    unsafe fn solve(MatrixLayout, UPLO, a: &[Self], &Pivot, b: &mut [Self]) -> Result<()>;
+    unsafe fn solveh(MatrixLayout, UPLO, a: &[Self], &Pivot, b: &mut [Self]) -> Result<()>;
 }
 
 macro_rules! impl_solveh {
@@ -30,13 +30,13 @@ impl Solveh_ for $scalar {
         into_result(info, ipiv)
     }
 
-    unsafe fn inv(l: MatrixLayout, uplo: UPLO, a: &mut [Self], ipiv: &Pivot) -> Result<()> {
+    unsafe fn invh(l: MatrixLayout, uplo: UPLO, a: &mut [Self], ipiv: &Pivot) -> Result<()> {
         let (n, _) = l.size();
         let info = $tri(l.lapacke_layout(), uplo as u8, n, a, l.lda(), ipiv);
         into_result(info, ())
     }
 
-    unsafe fn solve(l: MatrixLayout, uplo: UPLO, a: &[Self], ipiv: &Pivot, b: &mut [Self]) -> Result<()> {
+    unsafe fn solveh(l: MatrixLayout, uplo: UPLO, a: &[Self], ipiv: &Pivot, b: &mut [Self]) -> Result<()> {
         let (n, _) = l.size();
         let nrhs = 1;
         let ldb = 1;

src/lib.rs

@@ -42,6 +42,7 @@ pub mod operator;
 pub mod opnorm;
 pub mod qr;
 pub mod solve;
+pub mod solveh;
 pub mod svd;
 pub mod trace;
 pub mod triangular;
@@ -59,6 +60,7 @@ pub use operator::*;
 pub use opnorm::*;
 pub use qr::*;
 pub use solve::*;
+pub use solveh::*;
 pub use svd::*;
 pub use trace::*;
 pub use triangular::*;

src/solveh.rs

@@ -0,0 +1,125 @@
+//! Solve Hermite/Symmetric linear problems
+
+use ndarray::*;
+
+use super::convert::*;
+use super::error::*;
+use super::layout::*;
+use super::types::*;
+
+pub use lapack_traits::{Pivot, UPLO};
+
+pub struct FactorizedH<S: Data> {
+    pub a: ArrayBase<S, Ix2>,
+    pub ipiv: Pivot,
+}
+
+impl<A, S> FactorizedH<S>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    pub fn solveh<Sb>(&self, mut rhs: ArrayBase<Sb, Ix1>) -> Result<ArrayBase<Sb, Ix1>>
+    where
+        Sb: DataMut<Elem = A>,
+    {
+        unsafe {
+            A::solveh(
+                self.a.square_layout()?,
+                UPLO::Upper,
+                self.a.as_allocated()?,
+                &self.ipiv,
+                rhs.as_slice_mut().unwrap(),
+            )?
+        };
+        Ok(rhs)
+    }
+}
+
+impl<A, S> FactorizedH<S>
+where
+    A: Scalar,
+    S: DataMut<Elem = A>,
+{
+    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)
+    }
+}
+
+pub trait FactorizeH<S: Data> {
+    fn factorizeh(&self) -> Result<FactorizedH<S>>;
+}
+
+pub trait FactorizeHInto<S: Data> {
+    fn factorizeh_into(self) -> Result<FactorizedH<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>> {
+        let ipiv = unsafe { A::bk(self.layout()?, UPLO::Upper, self.as_allocated_mut()?)? };
+        Ok(FactorizedH {
+            a: self,
+            ipiv: ipiv,
+        })
+    }
+}
+
+impl<A, Si> FactorizeH<OwnedRepr<A>> for ArrayBase<Si, Ix2>
+where
+    A: Scalar,
+    Si: Data<Elem = A>,
+{
+    fn factorizeh(&self) -> Result<FactorizedH<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 })
+    }
+}
+
+pub trait InverseH {
+    type Output;
+    fn invh(&self) -> Result<Self::Output>;
+}
+
+pub trait InverseHInto {
+    type Output;
+    fn invh_into(self) -> Result<Self::Output>;
+}
+
+impl<A, S> InverseHInto for ArrayBase<S, Ix2>
+where
+    A: Scalar,
+    S: DataMut<Elem = A>,
+{
+    type Output = Self;
+
+    fn invh_into(self) -> Result<Self::Output> {
+        let f = self.factorizeh_into()?;
+        f.into_inverseh()
+    }
+}
+
+impl<A, Si> InverseH for ArrayBase<Si, Ix2>
+where
+    A: Scalar,
+    Si: Data<Elem = A>,
+{
+    type Output = Array2<A>;
+
+    fn invh(&self) -> Result<Self::Output> {
+        let f = self.factorizeh()?;
+        f.into_inverseh()
+    }
+}