diff --git a/src/dmat.rs b/src/dmat.rs
index 26bab478..366d154f 100644
--- a/src/dmat.rs
+++ b/src/dmat.rs
@@ -1,5 +1,7 @@
 //! Matrix with dimensions unknown at compile-time.
 
+#[doc(hidden)]; // we hide doc to not have to document the $trhs double dispatch trait.
+
 use std::rand::Rand;
 use std::rand;
 use std::num::{One, Zero};
@@ -7,28 +9,28 @@ use std::vec;
 use std::cmp::ApproxEq;
 use std::util;
 use dvec::{DVec, DVecMulRhs};
-use traits::operations::{Inv, Transpose};
+use traits::operations::{Inv, Transpose, Mean, Cov};
+
+#[doc(hidden)]
+mod metal;
 
 /// Matrix with dimensions unknown at compile-time.
-#[deriving(Eq, ToStr, Clone)]
+#[deriving(Eq, Clone)]
 pub struct DMat<N> {
     priv nrows: uint,
     priv ncols: uint,
     priv mij: ~[N]
 }
 
-/// Trait of object `o` which can be multiplied by a `DMat` `d`: `d * o`.
-pub trait DMatMulRhs<N, Res> {
-    /// Multiplies a `DMat` by `Self`.
-    fn binop(left: &DMat<N>, right: &Self) -> Res;
-}
+double_dispatch_binop_decl_trait!(DMat, DMatMulRhs)
+double_dispatch_binop_decl_trait!(DMat, DMatDivRhs)
+double_dispatch_binop_decl_trait!(DMat, DMatAddRhs)
+double_dispatch_binop_decl_trait!(DMat, DMatSubRhs)
 
-impl<N, Rhs: DMatMulRhs<N, Res>, Res> Mul<Rhs, Res> for DMat<N> {
-    #[inline(always)]
-    fn mul(&self, other: &Rhs) -> Res {
-        DMatMulRhs::binop(self, other)
-    }
-}
+mul_redispatch_impl!(DMat, DMatMulRhs)
+div_redispatch_impl!(DMat, DMatDivRhs)
+add_redispatch_impl!(DMat, DMatAddRhs)
+sub_redispatch_impl!(DMat, DMatSubRhs)
 
 impl<N> DMat<N> {
     /// Creates an uninitialized matrix.
@@ -89,6 +91,20 @@ impl<N: Clone> DMat<N> {
             mij:   vec::from_elem(nrows * ncols, val)
         }
     }
+
+    /// Builds a matrix filled with the components provided by a vector.
+    ///
+    /// The vector must have at least `nrows * ncols` elements.
+    #[inline]
+    pub fn from_vec(nrows: uint, ncols: uint, vec: &[N]) -> DMat<N> {
+        assert!(nrows * ncols <= vec.len());
+
+        DMat {
+            nrows: nrows,
+            ncols: ncols,
+            mij:   vec.slice_to(nrows * ncols).to_owned()
+        }
+    }
 }
 
 impl<N> DMat<N> {
@@ -116,7 +132,7 @@ impl<N> DMat<N> {
 
     /// Transforms this matrix into an array. This consumes the matrix and is O(1).
     #[inline]
-    pub fn to_array(self) -> ~[N] {
+    pub fn to_vec(self) -> ~[N] {
         self.mij
     }
 }
@@ -350,24 +366,88 @@ Inv for DMat<N> {
 impl<N: Clone> Transpose for DMat<N> {
     #[inline]
     fn transposed(&self) -> DMat<N> {
-        let mut res = self.clone();
+        if self.nrows == self.ncols {
+            let mut res = self.clone();
 
-        res.transpose();
+            res.transpose();
 
-        res
+            res
+        }
+        else {
+            let mut res = unsafe { DMat::new_uninitialized(self.ncols, self.nrows) };
+
+            for i in range(0u, self.nrows) {
+                for j in range(0u, self.ncols) {
+                    unsafe {
+                        res.set_fast(j, i, self.at_fast(i, j))
+                    }
+                }
+            }
+
+            res
+        }
     }
 
+    #[inline]
     fn transpose(&mut self) {
-        for i in range(1u, self.nrows) {
-            for j in range(0u, self.ncols - 1) {
-                let off_i_j = self.offset(i, j);
-                let off_j_i = self.offset(j, i);
+        if self.nrows == self.ncols {
+            for i in range(1u, self.nrows) {
+                for j in range(0u, self.ncols - 1) {
+                    let off_i_j = self.offset(i, j);
+                    let off_j_i = self.offset(j, i);
 
-                self.mij.swap(off_i_j, off_j_i);
+                    self.mij.swap(off_i_j, off_j_i);
+                }
+            }
+
+            util::swap(&mut self.nrows, &mut self.ncols);
+        }
+        else {
+            // FIXME: implement a better algorithm which does that in-place.
+            *self = self.transposed();
+        }
+    }
+}
+
+impl<N: Num + NumCast + Clone> Mean<DVec<N>> for DMat<N> {
+    fn mean(&self) -> DVec<N> {
+        let mut res: DVec<N> = DVec::new_zeros(self.ncols);
+        let normalizer: N    = NumCast::from(1.0f64 / NumCast::from(self.nrows));
+
+        for i in range(0u, self.nrows) {
+            for j in range(0u, self.ncols) {
+                unsafe {
+                    let acc = res.at_fast(j) + self.at_fast(i, j) * normalizer;
+                    res.set_fast(j, acc);
+                }
             }
         }
 
-        util::swap(&mut self.nrows, &mut self.ncols);
+        res
+    }
+}
+
+impl<N: Clone + Num + NumCast + DMatDivRhs<N, DMat<N>> + ToStr > Cov<DMat<N>> for DMat<N> {
+    // FIXME: this could be heavily optimized, removing all temporaries by merging loops.
+    fn cov(&self) -> DMat<N> {
+        assert!(self.nrows > 1);
+
+        let mut centered = unsafe { DMat::new_uninitialized(self.nrows, self.ncols) };
+        let mean = self.mean();
+
+        // FIXME: use the rows iterator when available
+        for i in range(0u, self.nrows) {
+            for j in range(0u, self.ncols) {
+                unsafe {
+                    centered.set_fast(i, j, self.at_fast(i, j) - mean.at_fast(j));
+                }
+            }
+        }
+
+        // FIXME: return a triangular matrix?
+        let normalizer: N = NumCast::from(self.nrows() - 1);
+        // FIXME: this will do 2 allocations for temporaries!
+        (centered.transposed() * centered) / normalizer
     }
 }
 
@@ -398,3 +478,137 @@ impl<N: ApproxEq<N>> ApproxEq<N> for DMat<N> {
         }
     }
 }
+
+macro_rules! scalar_mul_impl (
+    ($n: ident) => (
+        impl DMatMulRhs<$n, DMat<$n>> for $n {
+            #[inline]
+            fn binop(left: &DMat<$n>, right: &$n) -> DMat<$n> {
+                DMat {
+                    nrows: left.nrows,
+                    ncols: left.ncols,
+                    mij:   left.mij.iter().map(|a| a * *right).collect()
+                }
+            }
+        }
+    )
+)
+
+macro_rules! scalar_div_impl (
+    ($n: ident) => (
+        impl DMatDivRhs<$n, DMat<$n>> for $n {
+            #[inline]
+            fn binop(left: &DMat<$n>, right: &$n) -> DMat<$n> {
+                DMat {
+                    nrows: left.nrows,
+                    ncols: left.ncols,
+                    mij:   left.mij.iter().map(|a| a / *right).collect()
+                }
+            }
+        }
+    )
+)
+
+macro_rules! scalar_add_impl (
+    ($n: ident) => (
+        impl DMatAddRhs<$n, DMat<$n>> for $n {
+            #[inline]
+            fn binop(left: &DMat<$n>, right: &$n) -> DMat<$n> {
+                DMat {
+                    nrows: left.nrows,
+                    ncols: left.ncols,
+                    mij:   left.mij.iter().map(|a| a + *right).collect()
+                }
+            }
+        }
+    )
+)
+
+macro_rules! scalar_sub_impl (
+    ($n: ident) => (
+        impl DMatSubRhs<$n, DMat<$n>> for $n {
+            #[inline]
+            fn binop(left: &DMat<$n>, right: &$n) -> DMat<$n> {
+                DMat {
+                    nrows: left.nrows,
+                    ncols: left.ncols,
+                    mij:   left.mij.iter().map(|a| a - *right).collect()
+                }
+            }
+        }
+    )
+)
+
+scalar_mul_impl!(f64)
+scalar_mul_impl!(f32)
+scalar_mul_impl!(u64)
+scalar_mul_impl!(u32)
+scalar_mul_impl!(u16)
+scalar_mul_impl!(u8)
+scalar_mul_impl!(i64)
+scalar_mul_impl!(i32)
+scalar_mul_impl!(i16)
+scalar_mul_impl!(i8)
+scalar_mul_impl!(float)
+scalar_mul_impl!(uint)
+scalar_mul_impl!(int)
+
+scalar_div_impl!(f64)
+scalar_div_impl!(f32)
+scalar_div_impl!(u64)
+scalar_div_impl!(u32)
+scalar_div_impl!(u16)
+scalar_div_impl!(u8)
+scalar_div_impl!(i64)
+scalar_div_impl!(i32)
+scalar_div_impl!(i16)
+scalar_div_impl!(i8)
+scalar_div_impl!(float)
+scalar_div_impl!(uint)
+scalar_div_impl!(int)
+
+scalar_add_impl!(f64)
+scalar_add_impl!(f32)
+scalar_add_impl!(u64)
+scalar_add_impl!(u32)
+scalar_add_impl!(u16)
+scalar_add_impl!(u8)
+scalar_add_impl!(i64)
+scalar_add_impl!(i32)
+scalar_add_impl!(i16)
+scalar_add_impl!(i8)
+scalar_add_impl!(float)
+scalar_add_impl!(uint)
+scalar_add_impl!(int)
+
+scalar_sub_impl!(f64)
+scalar_sub_impl!(f32)
+scalar_sub_impl!(u64)
+scalar_sub_impl!(u32)
+scalar_sub_impl!(u16)
+scalar_sub_impl!(u8)
+scalar_sub_impl!(i64)
+scalar_sub_impl!(i32)
+scalar_sub_impl!(i16)
+scalar_sub_impl!(i8)
+scalar_sub_impl!(float)
+scalar_sub_impl!(uint)
+scalar_sub_impl!(int)
+
+impl<N: ToStr + Clone> ToStr for DMat<N> {
+    fn to_str(&self) -> ~str {
+        let mut res = ~"DMat ";
+        res = res + self.nrows.to_str() + " " + self.ncols.to_str() + " {\n";
+
+        for i in range(0u, self.nrows) {
+            for j in range(0u, self.ncols) {
+                res = res + " " + unsafe { self.at_fast(i, j).to_str() };
+            }
+
+            res = res + "\n";
+        }
+        res = res + "}";
+
+        res
+    }
+}
diff --git a/src/dvec.rs b/src/dvec.rs
index b77274c8..30d33fc8 100644
--- a/src/dvec.rs
+++ b/src/dvec.rs
@@ -12,6 +12,7 @@ use std::iter::FromIterator;
 use traits::geometry::{Dot, Norm, Translation};
 use traits::structure::{Iterable, IterableMut};
 
+#[doc(hidden)]
 mod metal;
 
 /// Vector with a dimension unknown at compile-time.
@@ -85,6 +86,16 @@ impl<N> DVec<N> {
             at: vec
         }
     }
+
+    #[inline]
+    pub unsafe fn set_fast(&mut self, i: uint, val: N) {
+        *self.at.unsafe_mut_ref(i) = val
+    }
+
+    #[inline]
+    pub fn to_vec(self) -> ~[N] {
+        self.at
+    }
 }
 
 impl<N: Clone> DVec<N> {
@@ -93,6 +104,18 @@ impl<N: Clone> DVec<N> {
     pub fn from_elem(dim: uint, elem: N) -> DVec<N> {
         DVec { at: vec::from_elem(dim, elem) }
     }
+
+    /// Builds a vector filled with the components provided by a vector.
+    ///
+    /// The vector must have at least `dim` elements.
+    #[inline]
+    pub fn from_vec(dim: uint, vec: &[N]) -> DVec<N> {
+        assert!(dim <= vec.len());
+
+        DVec {
+            at: vec.slice_to(dim).to_owned()
+        }
+    }
 }
 
 impl<N> DVec<N> {
diff --git a/src/mat.rs b/src/mat.rs
index 01f1bf0b..07e30da2 100644
--- a/src/mat.rs
+++ b/src/mat.rs
@@ -10,7 +10,7 @@ use vec::*;
 
 // traits
 pub use traits::structure::{Mat, Dim, Indexable, Iterable, IterableMut, MatCast, Row, Col};
-pub use traits::operations::{Absolute, ScalarSub, ScalarAdd, Inv, RMul, Transpose};
+pub use traits::operations::{Absolute, ScalarSub, ScalarAdd, Inv, RMul, Transpose, Mean, Cov};
 pub use traits::geometry::{Rotation, RotationMatrix, Rotate, Transformation, Transform,
                            Translation, Translate, ToHomogeneous, FromHomogeneous,
                            RotationWithTranslation, AbsoluteRotate};
diff --git a/src/metal.rs b/src/metal.rs
index 3e179e39..378c9996 100644
--- a/src/metal.rs
+++ b/src/metal.rs
@@ -1,5 +1,7 @@
 #[macro_escape];
 
+#[doc(hidden)]; // we hide doc to not have to document the $trhs double dispatch trait.
+
 // Create the traits needed to do fancy operator oveloading.
 // This is a meta version of
 // http://smallcultfollowing.com/babysteps/blog/2012/10/04/refining-traits-slash-impls/ 
diff --git a/src/tests/mat.rs b/src/tests/mat.rs
index 301d9753..18aca05f 100644
--- a/src/tests/mat.rs
+++ b/src/tests/mat.rs
@@ -87,12 +87,12 @@ fn test_inv_mat6() {
 
 #[test]
 fn test_rotation2() {
-  do 10000.times {
-    let randmat: Rotmat<Mat2<f64>> = One::one();
-    let ang     = &Vec1::new(abs::<f64>(random()) % Real::pi());
+    do 10000.times {
+        let randmat: Rotmat<Mat2<f64>> = One::one();
+        let ang     = &Vec1::new(abs::<f64>(random()) % Real::pi());
 
-    assert!(randmat.rotated(ang).rotation().approx_eq(ang));
-  }
+        assert!(randmat.rotated(ang).rotation().approx_eq(ang));
+    }
 }
 
 #[test]
@@ -104,12 +104,54 @@ fn test_index_mat2() {
 
 #[test]
 fn test_inv_rotation3() {
-  do 10000.times {
-    let randmat: Rotmat<Mat3<f64>> = One::one();
-    let dir: Vec3<f64> = random();
-    let ang            = &(dir.normalized() * (abs::<f64>(random()) % Real::pi()));
-    let rot            = randmat.rotated(ang);
+    do 10000.times {
+        let randmat: Rotmat<Mat3<f64>> = One::one();
+        let dir: Vec3<f64> = random();
+        let ang            = &(dir.normalized() * (abs::<f64>(random()) % Real::pi()));
+        let rot            = randmat.rotated(ang);
 
-    assert!((rot.transposed() * rot).approx_eq(&One::one()));
-  }
+        assert!((rot.transposed() * rot).approx_eq(&One::one()));
+    }
+}
+
+#[test]
+fn test_mean_dmat() {
+    let mat = DMat::from_vec(
+        3,
+        3,
+        [
+            1.0f64, 2.0, 3.0,
+            4.0f64, 5.0, 6.0,
+            7.0f64, 8.0, 9.0,
+        ]
+    );
+
+    assert!(mat.mean().approx_eq(&DVec::from_vec(3, [4.0f64, 5.0, 6.0])));
+}
+
+#[test]
+fn test_cov_dmat() {
+    let mat = DMat::from_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 expected = DMat::from_vec(
+        3,
+        3,
+        [
+            0.025,   0.0075,  0.00175,
+            0.0075,  0.007,   0.00135,
+            0.00175, 0.00135, 0.00043
+        ]
+    );
+
+    assert!(mat.cov().approx_eq(&expected));
 }
diff --git a/src/traits/operations.rs b/src/traits/operations.rs
index b2f59b7a..4575a8fa 100644
--- a/src/traits/operations.rs
+++ b/src/traits/operations.rs
@@ -31,6 +31,22 @@ pub trait Outer<V, M> {
     fn outer(&self, other: &V) -> M;
 }
 
+/// Trait for computing the covariance of a set of data.
+pub trait Cov<M> {
+    /// Computes the covariance of the obsevations stored by `self`:
+    ///     * for matrices, observations are stored in its rows.
+    ///     * for vectors, observations are stored in its components (thus are 1-dimensional).
+    fn cov(&self) -> M;
+}
+
+/// Trait for computing the covariance of a set of data.
+pub trait Mean<N> {
+    /// Computes the mean of the observations stored by `self`.
+    ///     * for matrices, observations are stored in its rows.
+    ///     * for vectors, observations are stored in its components (thus are 1-dimensional).
+    fn mean(&self) -> N;
+}
+
 // XXX: those two traits should not exist since there is generalized operator overloading of Add
 // and Sub.
 // However, using the same trait multiple time as a trait bound (ex: impl<T: Add<N, V> + Add<V, V>)
diff --git a/src/traits/structure.rs b/src/traits/structure.rs
index 64146416..7e105272 100644
--- a/src/traits/structure.rs
+++ b/src/traits/structure.rs
@@ -105,16 +105,24 @@ pub trait Row<R> {
     fn row(&self, i: uint) -> R;
     /// Writes the `i`-th row of `self`.
     fn set_row(&mut self, i: uint, R);
+
+    // FIXME: add iterators on rows: this could be a very good way to generalize _and_ optimize
+    // a lot of operations.
 }
 
 /// Trait to access columns of a matrix or vector.
 pub trait Col<C> {
     /// The number of column of this matrix or vector.
     fn num_cols(&self) -> uint;
+
     /// Reads the `i`-th column of `self`.
     fn col(&self, i: uint) -> C;
+
     /// Writes the `i`-th column of `self`.
     fn set_col(&mut self, i: uint, C);
+
+    // FIXME: add iterators on columns: this could be a very good way to generalize _and_ optimize
+    // a lot of operations.
 }
 
 /// Trait of objects having a spacial dimension known at compile time.