implemented QR factorization

this is a first sketch, the algorithm is not yet initialized and relies
on knowledge of DMat internals. A next step would be to implement this
algorithm in a more generic manner.

src/na.rs

@@ -44,6 +44,7 @@ pub use traits::{
 
 pub use structs::{
     Identity,
+    decomp_qr,
     DMat, DVec,
     Iso2, Iso3, Iso4,
     Mat1, Mat2, Mat3, Mat4,

src/structs/dmat.rs

@@ -4,12 +4,15 @@
 
 use rand::Rand;
 use rand;
-use std::num::{One, Zero};
+use std::num::{One, Zero, Float};
 use traits::operations::ApproxEq;
 use std::mem;
 use structs::dvec::{DVec, DVecMulRhs};
 use traits::operations::{Inv, Transpose, Mean, Cov};
 use traits::structure::Cast;
+use traits::geometry::Norm;
+use std::cmp::min;
+use std::fmt::{Show, Formatter, Result};
 
 #[doc(hidden)]
 mod metal;
@@ -494,6 +497,62 @@ impl<N: Clone + Num + Cast<f32> + DMatDivRhs<N, DMat<N>> + ToStr > Cov<DMat<N>>
     }
 }
 
+
+/// QR decomposition using Householder reflections
+/// # Arguments
+/// * `m` matrix to decompose
+pub fn decomp_qr<N: Clone + Num + Float + Show>(m: &DMat<N>) -> (DMat<N>, DMat<N>) {
+  let rows = m.nrows();
+  let cols = m.ncols();
+  assert!(rows >= cols);
+  let mut q : DMat<N> = DMat::new_identity(rows);
+  let mut r = m.clone();
+
+  let subtract_reflection = |vec: DVec<N>| -> DMat<N> {
+      // FIXME: we don't handle the complex case here
+      let mut qk : DMat<N> = DMat::new_identity(rows);
+      let start = rows - vec.at.len();
+      for j in range(start, rows) {
+          for i in range(start, rows) {
+              unsafe {
+                  let vv = vec.at_fast(i-start)*vec.at_fast(j-start);
+                  let qkij = qk.at_fast(i,j);
+                  qk.set_fast(i, j, qkij - vv - vv);
+              }
+          }
+      }
+      qk
+  };
+
+  let iterations = min(rows-1, cols);
+
+  for ite in range(0u, iterations) {
+      // we get the ite-th column truncated from its first ite elements,
+      // this is fast thanks to the matrix being column major
+      let start= m.offset(ite, ite);
+      let stop = m.offset(rows, ite);
+      let mut v = DVec::from_vec(rows - ite, r.mij.slice(start, stop));
+      let alpha =
+          if unsafe { v.at_fast(ite) } >= Zero::zero() {
+              -Norm::norm(&v)
+          }
+          else {
+              Norm::norm(&v)
+          };
+      unsafe {
+          let x = v.at_fast(0);
+          v.set_fast(0, x - alpha);
+      }
+      let _ = v.normalize();
+      let qk = subtract_reflection(v);
+      r = qk * r;
+      q = q * Transpose::transpose_cpy(&qk);
+  }
+
+  (q, r)
+}
+
+
 impl<N: ApproxEq<N>> ApproxEq<N> for DMat<N> {
     #[inline]
     fn approx_epsilon(_: Option<DMat<N>>) -> N {
@@ -515,6 +574,18 @@ impl<N: ApproxEq<N>> ApproxEq<N> for DMat<N> {
     }
 }
 
+impl<N: Show + Clone> Show for DMat<N> {
+    fn fmt(&self, form:&mut Formatter) -> Result {
+        for i in range(0u, self.nrows()) {
+            for j in range(0u, self.ncols()) {
+                write!(form.buf, "{} ", self.at(i, j));
+            }
+            write!(form.buf, "\n");
+        }
+        write!(form.buf, "\n")
+    }
+}
+
 macro_rules! scalar_mul_impl (
     ($n: ident) => (
         impl DMatMulRhs<$n, DMat<$n>> for $n {

src/structs/mod.rs

@@ -1,6 +1,7 @@
 //! Data structures and implementations.
 
 pub use self::dmat::DMat;
+pub use self::dmat::decomp_qr;
 pub use self::dvec::DVec;
 pub use self::vec::{Vec0, Vec1, Vec2, Vec3, Vec4, Vec5, Vec6};
 pub use self::mat::{Identity, Mat1, Mat2, Mat3, Mat4, Mat5, Mat6};

src/tests/mat.rs

@@ -2,6 +2,7 @@ use std::num::{Float, abs};
 use rand::random;
 use na::{Vec1, Vec3, Mat1, Mat2, Mat3, Mat4, Mat5, Mat6, Rot3, DMat, DVec, Indexable};
 use na;
+use na::decomp_qr;
 
 macro_rules! test_inv_mat_impl(
   ($t: ty) => (
@@ -206,3 +207,23 @@ fn test_dmat_from_vec() {
 
     assert!(mat1 == mat2);
 }
+
+#[test]
+fn test_decomp_qr() {
+    let mat = DMat::from_row_vec(
+        5,
+        3,
+        [
+            4.0, 2.0, 0.60,
+            4.2, 2.1, 0.59,
+            3.9, 2.0, 0.58,
+            4.3, 2.1, 0.62,
+            4.1, 2.2, 0.63
+        ]
+    );
+
+    let (q, r) = decomp_qr(&mat);
+    let mat_ = q * r;
+
+    assert!(na::approx_eq(&mat_,  &mat));
+}