Add wrapper for *{sy,he}{trf,tri,trs}

Author: termoshtt

src/lapack_traits/mod.rs

@@ -4,6 +4,7 @@ pub mod opnorm;
 pub mod qr;
 pub mod svd;
 pub mod solve;
+pub mod solveh;
 pub mod cholesky;
 pub mod eigh;
 pub mod triangular;
@@ -13,12 +14,15 @@ pub use self::eigh::*;
 pub use self::opnorm::*;
 pub use self::qr::*;
 pub use self::solve::*;
+pub use self::solveh::*;
 pub use self::svd::*;
 pub use self::triangular::*;
 
 use super::error::*;
 use super::types::*;
 
+pub type Pivot = Vec<i32>;
+
 pub trait LapackScalar
     : OperatorNorm_ + QR_ + SVD_ + Solve_ + Cholesky_ + Eigh_ + Triangular_ {
 }

src/lapack_traits/solve.rs

@@ -6,9 +6,7 @@ use error::*;
 use layout::MatrixLayout;
 use types::*;
 
-use super::{Transpose, into_result};
-
-pub type Pivot = Vec<i32>;
+use super::{Pivot, Transpose, into_result};
 
 /// Wraps `*getrf`, `*getri`, and `*getrs`
 pub trait Solve_: Sized {

src/lapack_traits/solveh.rs

@@ -0,0 +1,53 @@
+//! Solve symmetric linear problem using the Bunch-Kaufman diagonal pivoting method.
+//!
+//! See also [the manual of dsytrf](http://www.netlib.org/lapack/lapack-3.1.1/html/dsytrf.f.html)
+
+use lapack::c;
+
+use error::*;
+use layout::MatrixLayout;
+use types::*;
+
+use super::{Pivot, UPLO, into_result};
+
+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<()>;
+    /// Wrapper of `*sytrs` and `*hetrs`
+    unsafe fn solve(MatrixLayout, UPLO, a: &[Self], &Pivot, b: &mut [Self]) -> Result<()>;
+}
+
+macro_rules! impl_solveh {
+    ($scalar:ty, $trf:path, $tri:path, $trs:path) => {
+
+impl Solveh_ for $scalar {
+    unsafe fn bk(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<Pivot> {
+        let (n, _) = l.size();
+        let mut ipiv = vec![0; n as usize];
+        let info = $trf(l.lapacke_layout(), uplo as u8, n, a, l.lda(), &mut ipiv);
+        into_result(info, ipiv)
+    }
+
+    unsafe fn inv(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<()> {
+        let (n, _) = l.size();
+        let nrhs = 1;
+        let ldb = 1;
+        let info = $trs(l.lapacke_layout(), uplo as u8, n, nrhs, a, l.lda(), ipiv, b, ldb);
+        into_result(info, ())
+    }
+}
+
+}} // impl_solveh!
+
+impl_solveh!(f64, c::dsytrf, c::dsytri, c::dsytrs);
+impl_solveh!(f32, c::ssytrf, c::ssytri, c::ssytrs);
+impl_solveh!(c64, c::zhetrf, c::zhetri, c::zhetrs);
+impl_solveh!(c32, c::chetrf, c::chetri, c::chetrs);