QR factorization for fixed size matrices.

The ColSlice implementation for fixed size matrices returns a DVec,
while this is probably not optimal performance-wise, the dynamic nature
of the result makes this necessary. Using a data type presenting the
ImmutableVector trait would solve this, but it looks like a non-trivial
change.

src/structs/mat.rs

@@ -8,9 +8,10 @@ use traits::operations::ApproxEq;
 use std::slice::{Items, MutItems};
 use structs::vec::{Vec1, Vec2, Vec3, Vec4, Vec5, Vec6, Vec1MulRhs, Vec4MulRhs,
                    Vec5MulRhs, Vec6MulRhs};
+use structs::dvec::DVec;
 
 use traits::structure::{Cast, Row, Col, Iterable, IterableMut, Dim, Indexable,
-                        Eye};
+                        Eye, ColSlice, RowSlice};
 use traits::operations::{Absolute, Transpose, Inv, Outer};
 use traits::geometry::{ToHomogeneous, FromHomogeneous};
 
@@ -119,6 +120,8 @@ transpose_impl!(Mat1, 1)
 approx_eq_impl!(Mat1)
 row_impl!(Mat1, Vec1, 1)
 col_impl!(Mat1, Vec1, 1)
+col_slice_impl!(Mat1, Vec1, 1)
+row_slice_impl!(Mat1, Vec1, 1)
 to_homogeneous_impl!(Mat1, Mat2, 1, 2)
 from_homogeneous_impl!(Mat1, Mat2, 1, 2)
 outer_impl!(Vec1, Mat1)
@@ -218,6 +221,8 @@ transpose_impl!(Mat2, 2)
 approx_eq_impl!(Mat2)
 row_impl!(Mat2, Vec2, 2)
 col_impl!(Mat2, Vec2, 2)
+col_slice_impl!(Mat2, Vec2, 2)
+row_slice_impl!(Mat2, Vec2, 2)
 to_homogeneous_impl!(Mat2, Mat3, 2, 3)
 from_homogeneous_impl!(Mat2, Mat3, 2, 3)
 outer_impl!(Vec2, Mat2)
@@ -331,6 +336,8 @@ transpose_impl!(Mat3, 3)
 approx_eq_impl!(Mat3)
 // (specialized) row_impl!(Mat3, Vec3, 3)
 // (specialized) col_impl!(Mat3, Vec3, 3)
+col_slice_impl!(Mat3, Vec3, 3)
+row_slice_impl!(Mat3, Vec3, 3)
 to_homogeneous_impl!(Mat3, Mat4, 3, 4)
 from_homogeneous_impl!(Mat3, Mat4, 3, 4)
 outer_impl!(Vec3, Mat3)
@@ -496,6 +503,8 @@ transpose_impl!(Mat4, 4)
 approx_eq_impl!(Mat4)
 row_impl!(Mat4, Vec4, 4)
 col_impl!(Mat4, Vec4, 4)
+col_slice_impl!(Mat4, Vec4, 4)
+row_slice_impl!(Mat4, Vec4, 4)
 to_homogeneous_impl!(Mat4, Mat5, 4, 5)
 from_homogeneous_impl!(Mat4, Mat5, 4, 5)
 outer_impl!(Vec4, Mat4)
@@ -677,6 +686,8 @@ transpose_impl!(Mat5, 5)
 approx_eq_impl!(Mat5)
 row_impl!(Mat5, Vec5, 5)
 col_impl!(Mat5, Vec5, 5)
+col_slice_impl!(Mat5, Vec5, 5)
+row_slice_impl!(Mat5, Vec5, 5)
 to_homogeneous_impl!(Mat5, Mat6, 5, 6)
 from_homogeneous_impl!(Mat5, Mat6, 5, 6)
 outer_impl!(Vec5, Mat5)
@@ -910,4 +921,6 @@ transpose_impl!(Mat6, 6)
 approx_eq_impl!(Mat6)
 row_impl!(Mat6, Vec6, 6)
 col_impl!(Mat6, Vec6, 6)
+col_slice_impl!(Mat6, Vec6, 6)
+row_slice_impl!(Mat6, Vec6, 6)
 outer_impl!(Vec6, Mat6)

src/structs/mat_macros.rs

@@ -225,6 +225,26 @@ macro_rules! indexable_impl(
   )
 )
 
+macro_rules! col_slice_impl(
+    ($t: ident, $tv: ident, $dim: expr) => (
+        impl<N: Clone + Zero> ColSlice<DVec<N>> for $t<N> {
+            fn col_slice(& self, cid: uint, rstart: uint, rend: uint) -> DVec<N> {
+                assert!(cid < $dim);
+                assert!(rstart < rend);
+                assert!(rend <= $dim);
+                let col = self.col(cid);
+                let res  = DVec::from_vec(
+                    rend - rstart,
+                    unsafe {
+                        cast::transmute::<&$tv<N>, & [N, ..$dim]>
+                        (&col).slice(rstart, rend)
+                    });
+                res
+            }
+        }
+    )
+)
+
 macro_rules! row_impl(
   ($t: ident, $tv: ident, $dim: expr) => (
     impl<N: Clone + Zero> Row<$tv<N>> for $t<N> {
@@ -254,6 +274,26 @@ macro_rules! row_impl(
   )
 )
 
+macro_rules! row_slice_impl(
+    ($t: ident, $tv: ident, $dim: expr) => (
+        impl<N: Clone + Zero> RowSlice<DVec<N>> for $t<N> {
+            fn row_slice(& self, rid: uint, cstart: uint, cend: uint) -> DVec<N> {
+                assert!(rid < $dim);
+                assert!(cstart < cend);
+                assert!(cend <= $dim);
+                let row = self.row(rid);
+                let res  = DVec::from_vec(
+                    cend - cstart,
+                    unsafe {
+                        cast::transmute::<&$tv<N>, & [N, ..$dim]>
+                        (&row).slice(cstart, cend)
+                    });
+                res
+            }
+        }
+    )
+)
+
 macro_rules! col_impl(
   ($t: ident, $tv: ident, $dim: expr) => (
     impl<N: Clone + Zero> Col<$tv<N>> for $t<N> {

src/tests/mat.rs

@@ -3,6 +3,7 @@ use rand::random;
 use na::{Vec1, Vec3, Mat1, Mat2, Mat3, Mat4, Mat5, Mat6, Rot3, DMat, DVec, Indexable};
 use na;
 use na::decomp_qr;
+use std::cmp::{min, max};
 
 macro_rules! test_inv_mat_impl(
   ($t: ty) => (
@@ -24,6 +25,19 @@ macro_rules! test_transpose_mat_impl(
   );
 )
 
+macro_rules! test_decomp_qr_impl(
+  ($t: ty) => (
+    for _ in range(0, 10000) {
+      let randmat : $t = random();
+
+      let (q, r) = decomp_qr(&randmat);
+      let recomp = q * r;
+
+      assert!(na::approx_eq(&randmat,  &recomp));
+    }
+  );
+)
+
 #[test]
 fn test_transpose_mat1() {
     test_transpose_mat_impl!(Mat1<f64>);
@@ -210,20 +224,45 @@ fn test_dmat_from_vec() {
 
 #[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
-        ]
-    );
+    for _ in range(0, 10) {
+        let dim1: uint = random();
+        let dim2: uint = random();
+        let rows = min(40, max(dim1, dim2));
+        let cols = min(40, min(dim1, dim2));
+        let randmat: DMat<f64> = DMat::new_random(rows, cols);
+        let (q, r) = decomp_qr(&randmat);
+        let recomp = q * r;
 
-    let (q, r) = decomp_qr(&mat);
-    let mat_ = q * r;
-
-    assert!(na::approx_eq(&mat_,  &mat));
+        assert!(na::approx_eq(&randmat,  &recomp));
+    }
+}
+
+#[test]
+fn test_decomp_qr_mat1() {
+    test_decomp_qr_impl!(Mat1<f64>);
+}
+
+#[test]
+fn test_decomp_qr_mat2() {
+    test_decomp_qr_impl!(Mat2<f64>);
+}
+
+#[test]
+fn test_decomp_qr_mat3() {
+    test_decomp_qr_impl!(Mat3<f64>);
+}
+
+#[test]
+fn test_decomp_qr_mat4() {
+    test_decomp_qr_impl!(Mat4<f64>);
+}
+
+#[test]
+fn test_decomp_qr_mat5() {
+    test_decomp_qr_impl!(Mat5<f64>);
+}
+
+#[test]
+fn test_decomp_qr_mat6() {
+    test_decomp_qr_impl!(Mat6<f64>);
 }