Move the `eigen_qr` function behind the `EigenQR` trait.

This simplifies generic programming.

src/lib.rs

@@ -129,6 +129,7 @@ pub use traits::{
     Diag,
     Dim,
     Dot,
+    EigenQR,
     Eye,
     FloatPnt,
     FloatVec,
@@ -178,7 +179,6 @@ pub use structs::{
 
 pub use linalg::{
     qr,
-    eigen_qr,
     householder_matrix
 };
 
@@ -875,6 +875,15 @@ pub fn mean<N, M: Mean<N>>(observations: &M) -> N {
     Mean::mean(observations)
 }
 
+/*
+ * EigenQR<N, V>
+ */
+/// Computes the eigenvalues and eigenvectors of a square matrix usin the QR algorithm.
+#[inline(always)]
+pub fn eigen_qr<N, V, M: EigenQR<N, V>>(m: &M, eps: &N, niter: uint) -> (M, V) {
+    EigenQR::eigen_qr(m, eps, niter)
+}
+
 //
 //
 // Structure

src/linalg/decompositions.rs

@@ -76,12 +76,9 @@ pub fn qr<N, V, M>(m: &M) -> (M, M)
 pub fn eigen_qr<N, V, VS, M>(m: &M, eps: &N, niter: uint) -> (M, V)
     where N:  Float,
           VS: Indexable<uint, N> + Norm<N>,
-          M:  Indexable<(uint, uint), N> + SquareMat<N, V> + ColSlice<VS> + ApproxEq<N> + Clone {
+          M:  Indexable<(uint, uint), N> + SquareMat<N, V> + Add<M, M> + Sub<M, M> + ColSlice<VS> +
-    let (rows, cols) = m.shape();
+              ApproxEq<N> + Clone {
-
+    let mut eigenvectors: M = ::one::<M>();
-    assert!(rows == cols, "The matrix being decomposed must be square.");
-
-    let mut eigenvectors: M = Eye::new_identity(rows);
     let mut eigenvalues = m.clone();
     // let mut shifter: M = Eye::new_identity(rows);
 
@@ -89,7 +86,7 @@ pub fn eigen_qr<N, V, VS, M>(m: &M, eps: &N, niter: uint) -> (M, V)
     for _ in range(0, niter) {
         let mut stop = true;
 
-        for j in range(0, cols) {
+        for j in range(0, ::dim::<M>()) {
             for i in range(0, j) {
                 if unsafe { eigenvalues.unsafe_at((i, j)) }.abs() >= *eps {
                     stop = false;
@@ -97,7 +94,7 @@ pub fn eigen_qr<N, V, VS, M>(m: &M, eps: &N, niter: uint) -> (M, V)
                 }
             }
 
-            for i in range(j + 1, rows) {
+            for i in range(j + 1, ::dim::<M>()) {
                 if unsafe { eigenvalues.unsafe_at((i, j)) }.abs() >= *eps {
                     stop = false;
                     break;

src/structs/mat.rs

@@ -13,8 +13,9 @@ use structs::dvec::{DVec1, DVec2, DVec3, DVec4, DVec5, DVec6};
 
 use traits::structure::{Cast, Row, Col, Iterable, IterableMut, Dim, Indexable,
                         Eye, ColSlice, RowSlice, Diag, Shape};
-use traits::operations::{Absolute, Transpose, Inv, Outer};
+use traits::operations::{Absolute, Transpose, Inv, Outer, EigenQR};
 use traits::geometry::{ToHomogeneous, FromHomogeneous, Orig};
+use linalg;
 
 
 /// Special identity matrix. All its operation are no-ops.
@@ -128,6 +129,7 @@ diag_impl!(Mat1, Vec1, 1)
 to_homogeneous_impl!(Mat1, Mat2, 1, 2)
 from_homogeneous_impl!(Mat1, Mat2, 1, 2)
 outer_impl!(Vec1, Mat1)
+eigen_qr_impl!(Mat1, Vec1)
 
 /// Square matrix of dimension 2.
 #[deriving(Eq, PartialEq, Encodable, Decodable, Clone, Hash, Rand, Zero, Show)]
@@ -231,6 +233,7 @@ diag_impl!(Mat2, Vec2, 2)
 to_homogeneous_impl!(Mat2, Mat3, 2, 3)
 from_homogeneous_impl!(Mat2, Mat3, 2, 3)
 outer_impl!(Vec2, Mat2)
+eigen_qr_impl!(Mat2, Vec2)
 
 /// Square matrix of dimension 3.
 #[deriving(Eq, PartialEq, Encodable, Decodable, Clone, Hash, Rand, Zero, Show)]
@@ -348,6 +351,7 @@ diag_impl!(Mat3, Vec3, 3)
 to_homogeneous_impl!(Mat3, Mat4, 3, 4)
 from_homogeneous_impl!(Mat3, Mat4, 3, 4)
 outer_impl!(Vec3, Mat3)
+eigen_qr_impl!(Mat3, Vec3)
 
 /// Square matrix of dimension 4.
 #[deriving(Eq, PartialEq, Encodable, Decodable, Clone, Hash, Rand, Zero, Show)]
@@ -519,6 +523,7 @@ diag_impl!(Mat4, Vec4, 4)
 to_homogeneous_impl!(Mat4, Mat5, 4, 5)
 from_homogeneous_impl!(Mat4, Mat5, 4, 5)
 outer_impl!(Vec4, Mat4)
+eigen_qr_impl!(Mat4, Vec4)
 
 /// Square matrix of dimension 5.
 #[deriving(Eq, PartialEq, Encodable, Decodable, Clone, Hash, Rand, Zero, Show)]
@@ -706,6 +711,7 @@ diag_impl!(Mat5, Vec5, 5)
 to_homogeneous_impl!(Mat5, Mat6, 5, 6)
 from_homogeneous_impl!(Mat5, Mat6, 5, 6)
 outer_impl!(Vec5, Mat5)
+eigen_qr_impl!(Mat5, Vec5)
 
 /// Square matrix of dimension 6.
 #[deriving(Eq, PartialEq, Encodable, Decodable, Clone, Hash, Rand, Zero, Show)]
@@ -943,3 +949,4 @@ col_slice_impl!(Mat6, Vec6, DVec6, 6)
 row_slice_impl!(Mat6, Vec6, DVec6, 6)
 diag_impl!(Mat6, Vec6, 6)
 outer_impl!(Vec6, Mat6)
+eigen_qr_impl!(Mat6, Vec6)

src/structs/mat_macros.rs

@@ -638,3 +638,14 @@ macro_rules! outer_impl(
         }
     )
 )
+
+macro_rules! eigen_qr_impl(
+    ($t: ident, $v: ident) => (
+        impl<N> EigenQR<N, $v<N>> for $t<N>
+            where N: Float + ApproxEq<N> + Clone {
+            fn eigen_qr(m: &$t<N>, eps: &N, niter: uint) -> ($t<N>, $v<N>) {
+                linalg::eigen_qr(m, eps, niter)
+            }
+        }
+    )
+)

src/traits/mod.rs

@@ -9,7 +9,7 @@ pub use traits::structure::{FloatVec, FloatPnt, Basis, Cast, Col, Dim, Indexable
                             ColSlice, RowSlice, Diag, Eye, Shape};
 
 pub use traits::operations::{Absolute, ApproxEq, Axpy, Cov, Det, Inv, LMul, Mean, Outer, POrd,
-                             RMul, ScalarAdd, ScalarSub, ScalarMul, ScalarDiv, Transpose};
+                             RMul, ScalarAdd, ScalarSub, ScalarMul, ScalarDiv, Transpose, EigenQR};
 pub use traits::operations::{POrdering, PartialLess, PartialEqual, PartialGreater, NotComparable};
 
 pub mod geometry;

src/traits/operations.rs

@@ -1,6 +1,6 @@
 //! Low level operations on vectors and matrices.
 
-
+use traits::structure::SquareMat;
 
 /// Result of a partial ordering.
 #[deriving(Eq, PartialEq, Encodable, Decodable, Clone, Show)]
@@ -250,6 +250,12 @@ pub trait Mean<N> {
     fn mean(v: &Self) -> N;
 }
 
+/// Trait for computing the eigenvector and eigenvalues of a square matrix usin the QR algorithm.
+pub trait EigenQR<N, V>: SquareMat<N, V> {
+    /// Computes the eigenvectors and eigenvalues of this matrix.
+    fn eigen_qr(m: &Self, eps: &N, niter: uint) -> (Self, V);
+}
+
 
 // /// Cholesky decomposition.
 // pub trait Chol {

src/traits/structure.rs

@@ -1,8 +1,8 @@
 //! Traits giving structural informations on linear algebra objects or the space they live in.
 
-use std::num::Zero;
+use std::num::{Zero, One};
 use std::slice::{Items, MutItems};
-use traits::operations::{RMul, LMul, Axpy, Transpose};
+use traits::operations::{RMul, LMul, Axpy, Transpose, Inv};
 use traits::geometry::{Dot, Norm, Orig};
 
 /// Traits of objects which can be created from an object of type `T`.
@@ -21,12 +21,12 @@ impl<N, M, R, C> Mat<N, R, C> for M
 }
 
 /// Trait implemented by square matrices.
-pub trait SquareMat<N, V>: Mat<N, V, V> + Mul<Self, Self> + Eye + Transpose + Add<Self, Self> +
+pub trait SquareMat<N, V>: Mat<N, V, V> + Mul<Self, Self> + Eye + Transpose + Diag<V> + Inv + Dim +
-                           Sub<Self, Self> + Diag<V> {
+                           One {
 }
 
 impl<N, V, M> SquareMat<N, V> for M
-    where M: Mat<N, V, V> + Mul<M, M> + Eye + Transpose + Add<M, M> + Sub<M, M> + Diag<V> {
+    where M: Mat<N, V, V> + Mul<M, M> + Eye + Transpose + Diag<V> + Inv + Dim + One {
 }
 
 /// Trait for constructing the identity matrix