Add QR_

Author: termoshtt

src/impl2/mod.rs

@@ -1,6 +1,8 @@
 
 pub mod opnorm;
+pub mod qr;
 pub use self::opnorm::*;
+pub use self::qr::*;
 
 pub trait LapackScalar: OperatorNorm_ {}
 impl<A> LapackScalar for A where A: OperatorNorm_ {}

src/impl2/opnorm.rs

@@ -4,7 +4,7 @@ use lapack::c;
 use lapack::c::Layout::ColumnMajor as cm;
 
 use types::*;
-use layout::*;
+use layout::Layout;
 
 #[repr(u8)]
 pub enum NormType {

src/impl2/qr.rs

@@ -0,0 +1,55 @@
+//! Implement QR decomposition
+
+use std::cmp::min;
+use num_traits::Zero;
+use lapack::c;
+
+use types::*;
+use error::*;
+use layout::Layout;
+
+pub trait QR_: Sized {
+    fn householder(Layout, a: &mut [Self]) -> Result<Vec<Self>>;
+    fn q(Layout, a: &mut [Self], tau: &[Self]) -> Result<()>;
+    fn qr(Layout, a: &mut [Self]) -> Result<Vec<Self>>;
+}
+
+macro_rules! impl_qr {
+    ($scalar:ty, $qrf:path, $gqr:path) => {
+impl QR_ for $scalar {
+    fn householder(l: Layout, mut a: &mut [Self]) -> Result<Vec<Self>> {
+        let (row, col) = l.size();
+        let k = min(row, col);
+        let mut tau = vec![Self::zero(); k as usize];
+        let info = $qrf(l.lapacke_layout(), row, col, &mut a, l.lda(), &mut tau);
+        if info == 0 {
+            Ok(tau)
+        } else {
+            Err(LapackError::new(info).into())
+        }
+    }
+
+    fn q(l: Layout, mut a: &mut [Self], tau: &[Self]) -> Result<()> {
+        let (row, col) = l.size();
+        let k = min(row, col);
+        let info = $gqr(l.lapacke_layout(), row, k, k, &mut a, l.lda(), &tau);
+        if info == 0 {
+            Ok(())
+        } else {
+            Err(LapackError::new(info).into())
+        }
+    }
+
+    fn qr(l: Layout, mut a: &mut [Self]) -> Result<Vec<Self>> {
+        let tau = Self::householder(l, a)?;
+        let r = Vec::from(&*a);
+        Self::q(l, a, &tau)?;
+        Ok(r)
+    }
+}
+}} // endmacro
+
+impl_qr!(f64, c::dgeqrf, c::dorgqr);
+impl_qr!(f32, c::sgeqrf, c::sorgqr);
+impl_qr!(c64, c::zgeqrf, c::zungqr);
+impl_qr!(c32, c::cgeqrf, c::cungqr);

src/layout.rs

@@ -1,12 +1,14 @@
 
 use ndarray::*;
+use lapack::c;
 
 use super::error::*;
 
 pub type LDA = i32;
 pub type Col = i32;
 pub type Row = i32;
 
+#[derive(Debug, Clone, Copy)]
 pub enum Layout {
     C((Row, LDA)),
     F((Col, LDA)),
@@ -19,6 +21,28 @@ impl Layout {
             Layout::F((col, lda)) => (lda, col),
         }
     }
+
+    pub fn row(&self) -> Row {
+        self.size().0
+    }
+
+    pub fn col(&self) -> Col {
+        self.size().1
+    }
+
+    pub fn lda(&self) -> LDA {
+        match *self {
+            Layout::C((_, lda)) => lda,
+            Layout::F((_, lda)) => lda,
+        }
+    }
+
+    pub fn lapacke_layout(&self) -> c::Layout {
+        match *self {
+            Layout::C(_) => c::Layout::RowMajor,
+            Layout::F(_) => c::Layout::ColumnMajor,
+        }
+    }
 }
 
 pub trait AllocatedArray {