Add impl2/triangular.rs

Author: termoshtt

src/impl2/mod.rs

@@ -5,6 +5,7 @@ pub mod svd;
 pub mod solve;
 pub mod cholesky;
 pub mod eigh;
+pub mod triangular;
 
 pub use self::opnorm::*;
 pub use self::qr::*;
@@ -32,3 +33,11 @@ pub enum UPLO {
     Upper = b'U',
     Lower = b'L',
 }
+
+#[derive(Debug, Clone, Copy)]
+#[repr(u8)]
+pub enum Transpose {
+    No = b'N',
+    Transpose = b'T',
+    Hermite = b'C',
+}

src/impl2/solve.rs

@@ -5,18 +5,10 @@ use types::*;
 use error::*;
 use layout::Layout;
 
-use super::into_result;
+use super::{Transpose, into_result};
 
 pub type Pivot = Vec<i32>;
 
-#[derive(Debug, Clone, Copy)]
-#[repr(u8)]
-pub enum Transpose {
-    No = b'N',
-    Transpose = b'T',
-    Hermite = b'C',
-}
-
 pub trait Solve_: Sized {
     fn lu(Layout, a: &mut [Self]) -> Result<Pivot>;
     fn inv(Layout, a: &mut [Self], &Pivot) -> Result<()>;

src/impl2/triangular.rs

@@ -0,0 +1,52 @@
+//! Implement linear solver and inverse matrix
+
+use lapack::c;
+
+use error::*;
+use layout::Layout;
+use super::{UPLO, Transpose, into_result};
+
+#[derive(Debug, Clone, Copy)]
+#[repr(u8)]
+pub enum TriangularDiag {
+    Unit = b'U',
+    NonUnit = b'N',
+}
+
+pub trait Triangular_: Sized {
+    fn inv_triangular(l: Layout, UPLO, TriangularDiag, a: &mut [Self]) -> Result<()>;
+    fn solve_triangular(al: Layout, bl: Layout, UPLO, TriangularDiag, a: &[Self], b: &mut [Self]) -> Result<()>;
+}
+
+impl Triangular_ for f64 {
+    fn inv_triangular(l: Layout, uplo: UPLO, diag: TriangularDiag, a: &mut [Self]) -> Result<()> {
+        let (n, _) = l.size();
+        let lda = l.lda();
+        let info = c::dtrtri(l.lapacke_layout(), uplo as u8, diag as u8, n, a, lda);
+        into_result(info, ())
+    }
+
+    fn solve_triangular(al: Layout,
+                        bl: Layout,
+                        uplo: UPLO,
+                        diag: TriangularDiag,
+                        a: &[Self],
+                        mut b: &mut [Self])
+                        -> Result<()> {
+        let (n, _) = al.size();
+        let lda = al.lda();
+        let nrhs = bl.len();
+        let ldb = bl.lda();
+        let info = c::dtrtrs(al.lapacke_layout(),
+                             uplo as u8,
+                             Transpose::No as u8,
+                             diag as u8,
+                             n,
+                             nrhs,
+                             a,
+                             lda,
+                             &mut b,
+                             ldb);
+        into_result(info, ())
+    }
+}

src/layout.rs

@@ -5,6 +5,7 @@ use lapack::c;
 use super::error::*;
 
 pub type LDA = i32;
+pub type LEN = i32;
 pub type Col = i32;
 pub type Row = i32;
 
@@ -36,6 +37,13 @@ impl Layout {
         }
     }
 
+    pub fn len(&self) -> LEN {
+        match *self {
+            Layout::C((row, _)) => row,
+            Layout::F((col, _)) => col,
+        }
+    }
+
     pub fn lapacke_layout(&self) -> c::Layout {
         match *self {
             Layout::C(_) => c::Layout::RowMajor,