Use layout in QR factorization

Author: termoshtt

src/matrix.rs

@@ -96,19 +96,16 @@ impl<A> Matrix for Array<A, Ix2>
         let ua = Array::from_vec(u).into_shape((n, n)).unwrap();
         let va = Array::from_vec(vt).into_shape((m, m)).unwrap();
         match layout {
-            Layout::RowMajor => Ok((ua, sv, va)), 
+            Layout::RowMajor => Ok((ua, sv, va)),
             Layout::ColumnMajor => Ok((ua.reversed_axes(), sv, va.reversed_axes())),
         }
     }
     fn qr(self) -> Result<(Self, Self), LinalgError> {
         let (n, m) = self.size();
         let strides = self.strides();
         let k = min(n, m);
-        let (q, r) = if strides[0] < strides[1] {
-            try!(ImplQR::qr(m, n, self.clone().into_raw_vec()))
-        } else {
-            try!(ImplQR::lq(n, m, self.clone().into_raw_vec()))
-        };
+        let layout = self.layout()?;
+        let (q, r) = try!(ImplQR::qr(layout, m, n, self.clone().into_raw_vec()));
         let (qa, ra) = if strides[0] < strides[1] {
             (Array::from_vec(q).into_shape((m, n)).unwrap().reversed_axes(),
              Array::from_vec(r).into_shape((m, n)).unwrap().reversed_axes())

src/qr.rs

@@ -1,111 +1,41 @@
 //! Implement QR decomposition
 
 use std::cmp::min;
-use lapack::fortran::*;
+use lapack::c::*;
 use num_traits::Zero;
 
 use error::LapackError;
 
 pub trait ImplQR: Sized {
-    fn qr(n: usize, m: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError>;
-    fn lq(n: usize, m: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError>;
+    fn qr(layout: Layout, n: usize, m: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError>;
 }
 
 macro_rules! impl_qr {
-    ($geqrf:path, $orgqr:path, $gelqf:path, $orglq:path) => {
-// XXX These codes are most same, but the argument of $orgqr and $orglq are different!
-fn qr(n: usize, m: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
-    let n = n as i32;
-    let m = m as i32;
-    let mut info = 0;
-    let k = min(m, n);
-    let lda = m;
-    let lw_default = 1000;
-    let mut tau = vec![Self::zero(); k as usize];
-    let mut work = vec![Self::zero(); lw_default];
-// estimate lwork
-    $geqrf(m, n, &mut a, lda, &mut tau, &mut work, -1, &mut info);
-    let lwork_r = work[0] as i32;
-    if lwork_r > lw_default as i32 {
-        work = vec![Self::zero(); lwork_r as usize];
-    }
-// calc R
-    $geqrf(m, n, &mut a, lda, &mut tau, &mut work, lwork_r, &mut info);
-    if info != 0 {
-        return Err(From::from(info));
-    }
-    let r = a.clone();
-// re-estimate lwork
-    $orgqr(m, k, k, &mut a, lda, &mut tau, &mut work, -1, &mut info);
-    let lwork_q = work[0] as i32;
-    if lwork_q > lwork_r {
-        work = vec![Self::zero(); lwork_q as usize];
-    }
-// calc Q
-    $orgqr(m,
-           k,
-           k,
-           &mut a,
-           lda,
-           &mut tau,
-           &mut work,
-           lwork_q,
-           &mut info);
-    if info == 0 {
-        Ok((a, r))
-    } else {
-        Err(From::from(info))
-    }
-}
-fn lq(n: usize, m: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
-    let n = n as i32;
-    let m = m as i32;
-    let mut info = 0;
-    let k = min(m, n);
-    let lda = m;
-    let lw_default = 1000;
-    let mut tau = vec![Self::zero(); k as usize];
-    let mut work = vec![Self::zero(); lw_default];
-// estimate lwork
-    $gelqf(m, n, &mut a, lda, &mut tau, &mut work, -1, &mut info);
-    let lwork_r = work[0] as i32;
-    if lwork_r > lw_default as i32 {
-        work = vec![Self::zero(); lwork_r as usize];
-    }
-// calc R
-    $gelqf(m, n, &mut a, lda, &mut tau, &mut work, lwork_r, &mut info);
-    if info != 0 {
-        return Err(From::from(info));
-    }
-    let r = a.clone();
-// re-estimate lwork
-    $orglq(k, n, k, &mut a, lda, &mut tau, &mut work, -1, &mut info);
-    let lwork_q = work[0] as i32;
-    if lwork_q > lwork_r {
-        work = vec![Self::zero(); lwork_q as usize];
-    }
-// calc Q
-    $orglq(k,
-           n,
-           k,
-           &mut a,
-           lda,
-           &mut tau,
-           &mut work,
-           lwork_q,
-           &mut info);
-    if info == 0 {
-        Ok((a, r))
-    } else {
-        Err(From::from(info))
+    ($scalar:ty, $geqrf:path, $orgqr:path) => {
+impl ImplQR for $scalar {
+    fn qr(layout: Layout, n: usize, m: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
+        let n = n as i32;
+        let m = m as i32;
+        let k = min(m, n);
+        let lda = match layout {
+            Layout::ColumnMajor => m,
+            Layout::RowMajor => n,
+        };
+        let mut tau = vec![Self::zero(); k as usize];
+        let info = $geqrf(layout, m, n, &mut a, lda, &mut tau);
+        if info != 0 {
+            return Err(From::from(info));
+        }
+        let r = a.clone();
+        let info = $orgqr(layout, m, k, k, &mut a, lda, &mut tau);
+        if info == 0 {
+            Ok((a, r))
+        } else {
+            Err(From::from(info))
+        }
     }
 }
 }} // endmacro
 
-impl ImplQR for f64 {
-    impl_qr!(dgeqrf, dorgqr, dgelqf, dorglq);
-}
-
-impl ImplQR for f32 {
-    impl_qr!(sgeqrf, sorgqr, sgelqf, sorglq);
-}
+impl_qr!(f64, dgeqrf, dorgqr);
+impl_qr!(f32, sgeqrf, sorgqr);