Triangular matrix

src/generate.rs

@@ -19,8 +19,18 @@ pub fn conjugate<A, Si, So>(a: &ArrayBase<Si, Ix2>) -> ArrayBase<So, Ix2>
     a
 }
 
+/// Random vector
+pub fn random_vector<A, S>(n: usize) -> ArrayBase<S, Ix1>
+    where A: RandNormal,
+          S: DataOwned<Elem = A>
+{
+    let mut rng = thread_rng();
+    let v: Vec<A> = (0..n).map(|_| A::randn(&mut rng)).collect();
+    ArrayBase::from_vec(v)
+}
+
 /// Random matrix
-pub fn random<A, S>(n: usize, m: usize) -> ArrayBase<S, Ix2>
+pub fn random_matrix<A, S>(n: usize, m: usize) -> ArrayBase<S, Ix2>
     where A: RandNormal,
           S: DataOwned<Elem = A>
 {
@@ -34,7 +44,7 @@ pub fn random_square<A, S>(n: usize) -> ArrayBase<S, Ix2>
     where A: RandNormal,
           S: DataOwned<Elem = A>
 {
-    random(n, n)
+    random_matrix(n, n)
 }
 
 /// Random Hermite matrix

src/impl2/mod.rs

@@ -5,13 +5,15 @@ pub mod svd;
 pub mod solve;
 pub mod cholesky;
 pub mod eigh;
+pub mod triangular;
 
 pub use self::opnorm::*;
 pub use self::qr::*;
 pub use self::svd::*;
 pub use self::solve::*;
 pub use self::cholesky::*;
 pub use self::eigh::*;
+pub use self::triangular::*;
 
 use super::error::*;
 
@@ -20,7 +22,8 @@ trait_alias!(LapackScalar: OperatorNorm_,
              SVD_,
              Solve_,
              Cholesky_,
-             Eigh_);
+             Eigh_,
+             Triangular_);
 
 pub fn into_result<T>(info: i32, val: T) -> Result<T> {
     if info == 0 {
@@ -36,3 +39,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,54 @@
+//! Implement linear solver and inverse matrix
+
+use lapack::c;
+
+use error::*;
+use types::*;
+use layout::Layout;
+use super::{UPLO, Transpose, into_result};
+
+#[derive(Debug, Clone, Copy)]
+#[repr(u8)]
+pub enum Diag {
+    Unit = b'U',
+    NonUnit = b'N',
+}
+
+pub trait Triangular_: Sized {
+    fn inv_triangular(l: Layout, UPLO, Diag, a: &mut [Self]) -> Result<()>;
+    fn solve_triangular(al: Layout, bl: Layout, UPLO, Diag, a: &[Self], b: &mut [Self]) -> Result<()>;
+}
+
+macro_rules! impl_triangular {
+    ($scalar:ty, $trtri:path, $trtrs:path) => {
+
+impl Triangular_ for $scalar {
+    fn inv_triangular(l: Layout, uplo: UPLO, diag: Diag, a: &mut [Self]) -> Result<()> {
+        let (n, _) = l.size();
+        let lda = l.lda();
+        let info = $trtri(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: Diag, a: &[Self], mut b: &mut [Self]) -> Result<()> {
+        let (n, _) = al.size();
+        let lda = al.lda();
+        let (_, nrhs) = bl.size();
+        let ldb = bl.lda();
+        println!("al = {:?}", al);
+        println!("bl = {:?}", bl);
+        println!("n = {}", n);
+        println!("lda = {}", lda);
+        println!("nrhs = {}", nrhs);
+        println!("ldb = {}", ldb);
+        let info = $trtrs(al.lapacke_layout(), uplo as u8, Transpose::No as u8, diag as u8, n, nrhs, a, lda, &mut b, ldb);
+        into_result(info, ())
+    }
+}
+
+}} // impl_triangular!
+
+impl_triangular!(f64, c::dtrtri, c::dtrtrs);
+impl_triangular!(f32, c::strtri, c::strtrs);
+impl_triangular!(c64, c::ztrtri, c::ztrtrs);
+impl_triangular!(c32, c::ctrtri, c::ctrtrs);