diff --git a/Makefile b/Makefile
index 5c00b6af..3e961e87 100644
--- a/Makefile
+++ b/Makefile
@@ -3,13 +3,17 @@ nalgebra_lib_path=lib
 nalgebra_doc_path=doc
 all:
 	mkdir -p $(nalgebra_lib_path)
-	rust build src/lib.rs --out-dir $(nalgebra_lib_path) --opt-level 3
+	rust build src/lib.rs --out-dir $(nalgebra_lib_path) --opt-level 2
 
 test:
 	mkdir -p $(nalgebra_lib_path)
 	rust test src/lib.rs
 	rm libtest~
 
+bench:
+	rustc --test src/lib.rs --opt-level 2 -o bench~ && ./bench~ --bench
+	rm bench~
+
 doc:
 	mkdir -p $(nalgebra_doc_path)
 	rust doc src/lib.rs --output-dir $(nalgebra_doc_path)
@@ -18,6 +22,7 @@ distcheck:
 	rm -rf $(tmp)
 	git clone --recursive . $(tmp)
 	make -C $(tmp)
+	make -C $(tmp) test
 	rm -rf $(tmp)
 
 .PHONY:doc
diff --git a/src/bench/mat.rs b/src/bench/mat.rs
new file mode 100644
index 00000000..e2ff6486
--- /dev/null
+++ b/src/bench/mat.rs
@@ -0,0 +1,164 @@
+use std::rand::random;
+use extra::test::BenchHarness;
+use mat::*;
+use vec::*;
+
+macro_rules! bench_mul_mat(
+    ($bh: expr, $t: ty) => {
+        {
+            let a: $t     = random();
+            let mut b: $t = random();
+
+            do $bh.iter {
+                do 1000.times {
+                    b = a * b;
+                }
+            }
+        }
+    }
+)
+
+#[bench]
+fn bench_mul_mat2(bh: &mut BenchHarness) {
+    bench_mul_mat!(bh, Mat2<f64>)
+}
+
+#[bench]
+fn bench_mul_mat3(bh: &mut BenchHarness) {
+    bench_mul_mat!(bh, Mat3<f64>)
+}
+
+#[bench]
+fn bench_mul_mat4(bh: &mut BenchHarness) {
+    bench_mul_mat!(bh, Mat4<f64>)
+}
+
+#[bench]
+fn bench_mul_mat5(bh: &mut BenchHarness) {
+    bench_mul_mat!(bh, Mat5<f64>)
+}
+
+#[bench]
+fn bench_mul_mat6(bh: &mut BenchHarness) {
+    bench_mul_mat!(bh, Mat6<f64>)
+}
+
+macro_rules! bench_mul_dmat(
+    ($bh: expr, $nrows: expr, $ncols: expr) => {
+        {
+            let a:     DMat<f64> = DMat::new_random($nrows, $ncols);
+            let mut b: DMat<f64> = DMat::new_random($nrows, $ncols);
+
+            do $bh.iter {
+                do 1000.times {
+                    b = a * b;
+                }
+            }
+        }
+    }
+)
+
+#[bench]
+fn bench_mul_dmat2(bh: &mut BenchHarness) {
+    bench_mul_dmat!(bh, 2, 2)
+}
+
+#[bench]
+fn bench_mul_dmat3(bh: &mut BenchHarness) {
+    bench_mul_dmat!(bh, 3, 3)
+}
+
+#[bench]
+fn bench_mul_dmat4(bh: &mut BenchHarness) {
+    bench_mul_dmat!(bh, 4, 4)
+}
+
+#[bench]
+fn bench_mul_dmat5(bh: &mut BenchHarness) {
+    bench_mul_dmat!(bh, 5, 5)
+}
+
+#[bench]
+fn bench_mul_dmat6(bh: &mut BenchHarness) {
+    bench_mul_dmat!(bh, 6, 6)
+}
+
+macro_rules! bench_mul_mat_vec(
+    ($bh: expr, $tm: ty, $tv: ty) => {
+        {
+            let m : $tm     = random();
+            let mut v : $tv = random();
+
+            do $bh.iter {
+                do 1000.times {
+                    v = m.rmul(&v)
+                }
+            }
+        }
+    }
+)
+
+#[bench]
+fn bench_mul_mat_vec2(bh: &mut BenchHarness) {
+    bench_mul_mat_vec!(bh, Mat2<f64>, Vec2<f64>)
+}
+
+#[bench]
+fn bench_mul_mat_vec3(bh: &mut BenchHarness) {
+    bench_mul_mat_vec!(bh, Mat3<f64>, Vec3<f64>)
+}
+
+#[bench]
+fn bench_mul_mat_vec4(bh: &mut BenchHarness) {
+    bench_mul_mat_vec!(bh, Mat4<f64>, Vec4<f64>)
+}
+
+#[bench]
+fn bench_mul_mat_vec5(bh: &mut BenchHarness) {
+    bench_mul_mat_vec!(bh, Mat5<f64>, Vec5<f64>)
+}
+
+#[bench]
+fn bench_mul_mat_vec6(bh: &mut BenchHarness) {
+    bench_mul_mat_vec!(bh, Mat6<f64>, Vec6<f64>)
+}
+
+macro_rules! bench_mul_dmat_dvec(
+    ($bh: expr, $nrows: expr, $ncols: expr) => {
+        {
+            let m : DMat<f64>     = DMat::new_random($nrows, $ncols);
+            let mut v : DVec<f64> = DVec::new_random($ncols);
+
+            do $bh.iter {
+                do 1000.times {
+                    v = m.rmul(&v)
+                }
+            }
+        }
+    }
+)
+
+#[bench]
+fn bench_mul_dmat_dvec2(bh: &mut BenchHarness) {
+    bench_mul_dmat_dvec!(bh, 2, 2)
+}
+
+#[bench]
+fn bench_mul_dmat_dvec3(bh: &mut BenchHarness) {
+    bench_mul_dmat_dvec!(bh, 3, 3)
+}
+
+#[bench]
+fn bench_mul_dmat_dvec4(bh: &mut BenchHarness) {
+    bench_mul_dmat_dvec!(bh, 4, 4)
+}
+
+#[bench]
+fn bench_mul_dmat_dvec5(bh: &mut BenchHarness) {
+    bench_mul_dmat_dvec!(bh, 5, 5)
+}
+
+#[bench]
+fn bench_mul_dmat_dvec6(bh: &mut BenchHarness) {
+    bench_mul_dmat_dvec!(bh, 6, 6)
+}
diff --git a/src/bench/vec.rs b/src/bench/vec.rs
new file mode 100644
index 00000000..59fc6a8f
--- /dev/null
+++ b/src/bench/vec.rs
@@ -0,0 +1,44 @@
+use std::rand::random;
+use extra::test::BenchHarness;
+use vec::*;
+
+macro_rules! bench_dot_vec(
+    ($bh: expr, $t: ty) => {
+        {
+            let a: $t = random();
+            let b: $t = random();
+            let mut d = 0.0;
+
+            do $bh.iter {
+                do 1000.times {
+                    d = d + a.dot(&b);
+                }
+            }
+        }
+    }
+)
+
+#[bench]
+fn bench_dot_vec2(bh: &mut BenchHarness) {
+    bench_dot_vec!(bh, Vec2<f64>)
+}
+
+#[bench]
+fn bench_dot_vec3(bh: &mut BenchHarness) {
+    bench_dot_vec!(bh, Vec3<f64>)
+}
+
+#[bench]
+fn bench_dot_vec4(bh: &mut BenchHarness) {
+    bench_dot_vec!(bh, Vec4<f64>)
+}
+
+#[bench]
+fn bench_dot_vec5(bh: &mut BenchHarness) {
+    bench_dot_vec!(bh, Vec5<f64>)
+}
+
+#[bench]
+fn bench_dot_vec6(bh: &mut BenchHarness) {
+    bench_dot_vec!(bh, Vec6<f64>)
+}
diff --git a/src/lib.rs b/src/lib.rs
index df4f3665..8b36d633 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -76,5 +76,10 @@ pub mod traits {
 mod tests {
     mod mat;
     mod vec;
-    // mod bench;
+}
+
+#[cfg(test)]
+mod bench {
+    mod mat;
+    mod vec;
 }
diff --git a/src/mat.rs b/src/mat.rs
index fd1f233e..6ce49761 100644
--- a/src/mat.rs
+++ b/src/mat.rs
@@ -66,6 +66,7 @@ absolute_impl!(Mat1, m11)
 one_impl!(Mat1, One::one)
 iterable_impl!(Mat1, 1)
 iterable_mut_impl!(Mat1, 1)
+at_fast_impl!(Mat1, 1)
 dim_impl!(Mat1, 1)
 indexable_impl!(Mat1, 1)
 mul_impl!(Mat1, 1)
@@ -113,6 +114,7 @@ iterable_impl!(Mat2, 2)
 iterable_mut_impl!(Mat2, 2)
 dim_impl!(Mat2, 2)
 indexable_impl!(Mat2, 2)
+at_fast_impl!(Mat2, 2)
 mul_impl!(Mat2, 2)
 rmul_impl!(Mat2, Vec2, 2)
 lmul_impl!(Mat2, Vec2, 2)
@@ -169,6 +171,7 @@ iterable_impl!(Mat3, 3)
 iterable_mut_impl!(Mat3, 3)
 dim_impl!(Mat3, 3)
 indexable_impl!(Mat3, 3)
+at_fast_impl!(Mat3, 3)
 mul_impl!(Mat3, 3)
 rmul_impl!(Mat3, Vec3, 3)
 lmul_impl!(Mat3, Vec3, 3)
@@ -254,6 +257,7 @@ iterable_impl!(Mat4, 4)
 iterable_mut_impl!(Mat4, 4)
 dim_impl!(Mat4, 4)
 indexable_impl!(Mat4, 4)
+at_fast_impl!(Mat4, 4)
 mul_impl!(Mat4, 4)
 rmul_impl!(Mat4, Vec4, 4)
 lmul_impl!(Mat4, Vec4, 4)
@@ -352,6 +356,7 @@ iterable_impl!(Mat5, 5)
 iterable_mut_impl!(Mat5, 5)
 dim_impl!(Mat5, 5)
 indexable_impl!(Mat5, 5)
+at_fast_impl!(Mat5, 5)
 mul_impl!(Mat5, 5)
 rmul_impl!(Mat5, Vec5, 5)
 lmul_impl!(Mat5, Vec5, 5)
@@ -461,6 +466,7 @@ iterable_impl!(Mat6, 6)
 iterable_mut_impl!(Mat6, 6)
 dim_impl!(Mat6, 6)
 indexable_impl!(Mat6, 6)
+at_fast_impl!(Mat6, 6)
 mul_impl!(Mat6, 6)
 rmul_impl!(Mat6, Vec6, 6)
 lmul_impl!(Mat6, Vec6, 6)
diff --git a/src/mat_macros.rs b/src/mat_macros.rs
index aaa22fb4..70538667 100644
--- a/src/mat_macros.rs
+++ b/src/mat_macros.rs
@@ -14,6 +14,24 @@ macro_rules! mat_impl(
   )
 )
 
+macro_rules! at_fast_impl(
+    ($t: ident, $dim: expr) => (
+        impl<N: Clone> $t<N> {
+            #[inline]
+            unsafe fn at_fast(&self, (i, j): (uint, uint)) -> N {
+                (*cast::transmute::<&$t<N>, &[N, ..$dim * $dim]>(self)
+                 .unsafe_ref(i * $dim + j)).clone()
+            }
+
+            #[inline]
+            unsafe fn set_fast(&mut self, (i, j): (uint, uint), val: N) {
+                (*cast::transmute::<&mut $t<N>, &mut [N, ..$dim * $dim]>(self)
+                 .unsafe_mut_ref(i * $dim + j)) = val
+            }
+        }
+    )
+)
+
 macro_rules! mat_cast_impl(
   ($t: ident, $comp0: ident $(,$compN: ident)*) => (
     impl<Nin: NumCast + Clone, Nout: NumCast> MatCast<$t<Nout>> for $t<Nin> {
@@ -285,11 +303,13 @@ macro_rules! mul_impl(
                 for j in range(0u, $dim) {
                     let mut acc: N = Zero::zero();
 
-                    for k in range(0u, $dim) {
-                        acc = acc + self.at((i, k)) * other.at((k, j));
-                    }
+                    unsafe {
+                        for k in range(0u, $dim) {
+                            acc = acc + self.at_fast((i, k)) * other.at_fast((k, j));
+                        }
 
-                    res.set((i, j), acc);
+                        res.set_fast((i, j), acc);
+                    }
                 }
             }
 
@@ -308,8 +328,10 @@ macro_rules! rmul_impl(
 
             for i in range(0u, $dim) {
                 for j in range(0u, $dim) {
-                    let val = res.at(i) + other.at(j) * self.at((i, j));
-                    res.set(i, val)
+                    unsafe {
+                        let val = res.at_fast(i) + other.at_fast(j) * self.at_fast((i, j));
+                        res.set_fast(i, val)
+                    }
                 }
             }
 
@@ -328,8 +350,10 @@ macro_rules! lmul_impl(
 
             for i in range(0u, $dim) {
                 for j in range(0u, $dim) {
-                  let val = res.at(i) + other.at(j) * self.at((j, i));
-                  res.set(i, val)
+                    unsafe {
+                        let val = res.at_fast(i) + other.at_fast(j) * self.at_fast((j, i));
+                        res.set_fast(i, val)
+                    }
                 }
             }
 
diff --git a/src/tests/mat.rs b/src/tests/mat.rs
index 4f8fb69d..301d9753 100644
--- a/src/tests/mat.rs
+++ b/src/tests/mat.rs
@@ -1,13 +1,8 @@
 use std::num::{Real, One, abs};
 use std::rand::random;
 use std::cmp::ApproxEq;
-use traits::inv::Inv;
-use traits::rotation::Rotation;
-use traits::indexable::Indexable;
-use traits::transpose::Transpose;
-use traits::norm::Norm;
-use vec::{Vec1, Vec3};
-use mat::{Mat1, Mat2, Mat3, Mat4, Mat5, Mat6};
+use vec::*;
+use mat::*;
 use adaptors::rotmat::Rotmat;
 
 macro_rules! test_inv_mat_impl(
diff --git a/src/tests/vec.rs b/src/tests/vec.rs
index 003efb0f..cf7064ab 100644
--- a/src/tests/vec.rs
+++ b/src/tests/vec.rs
@@ -1,15 +1,8 @@
 use std::num::{Zero, One};
 use std::rand::{random};
 use std::cmp::ApproxEq;
-use vec::{Vec0, Vec1, Vec2, Vec3, Vec4, Vec5, Vec6};
-use mat::Mat3;
-use traits::basis::Basis;
-use traits::cross::Cross;
-use traits::dot::Dot;
-use traits::norm::Norm;
-use traits::iterable::{Iterable, IterableMut};
-use traits::scalar_op::{ScalarAdd, ScalarSub};
-use traits::outer::Outer;
+use vec::*;
+use mat::*;
 
 macro_rules! test_iterator_impl(
     ($t: ty, $n: ty) => (
diff --git a/src/vec.rs b/src/vec.rs
index 5f11ef29..cc70ba6f 100644
--- a/src/vec.rs
+++ b/src/vec.rs
@@ -45,6 +45,7 @@ orderable_impl!(Vec1, x)
 vec_axis_impl!(Vec1, x)
 vec_cast_impl!(Vec1, x)
 indexable_impl!(Vec1, 1)
+at_fast_impl!(Vec1, 1)
 new_repeat_impl!(Vec1, val, x)
 dim_impl!(Vec1, 1)
 container_impl!(Vec1)
@@ -87,6 +88,7 @@ orderable_impl!(Vec2, x, y)
 vec_axis_impl!(Vec2, x, y)
 vec_cast_impl!(Vec2, x, y)
 indexable_impl!(Vec2, 2)
+at_fast_impl!(Vec2, 2)
 new_repeat_impl!(Vec2, val, x, y)
 dim_impl!(Vec2, 2)
 container_impl!(Vec2)
@@ -131,6 +133,7 @@ orderable_impl!(Vec3, x, y, z)
 vec_axis_impl!(Vec3, x, y, z)
 vec_cast_impl!(Vec3, x, y, z)
 indexable_impl!(Vec3, 3)
+at_fast_impl!(Vec3, 3)
 new_repeat_impl!(Vec3, val, x, y, z)
 dim_impl!(Vec3, 3)
 container_impl!(Vec3)
@@ -184,6 +187,7 @@ orderable_impl!(PVec3, x, y, z)
 vec_axis_impl!(PVec3, x, y, z)
 vec_cast_impl!(PVec3, x, y, z)
 indexable_impl!(PVec3, 3)
+at_fast_impl!(PVec3, 3)
 new_repeat_impl!(PVec3, val, x, y, z, _unused)
 dim_impl!(PVec3, 3)
 container_impl!(PVec3)
@@ -232,6 +236,7 @@ orderable_impl!(Vec4, x, y, z, w)
 vec_axis_impl!(Vec4, x, y, z, w)
 vec_cast_impl!(Vec4, x, y, z, w)
 indexable_impl!(Vec4, 4)
+at_fast_impl!(Vec4, 4)
 new_repeat_impl!(Vec4, val, x, y, z, w)
 dim_impl!(Vec4, 4)
 container_impl!(Vec4)
@@ -280,6 +285,7 @@ orderable_impl!(Vec5, x, y, z, w, a)
 vec_axis_impl!(Vec5, x, y, z, w, a)
 vec_cast_impl!(Vec5, x, y, z, w, a)
 indexable_impl!(Vec5, 5)
+at_fast_impl!(Vec5, 5)
 new_repeat_impl!(Vec5, val, x, y, z, w, a)
 dim_impl!(Vec5, 5)
 container_impl!(Vec5)
@@ -330,6 +336,7 @@ orderable_impl!(Vec6, x, y, z, w, a, b)
 vec_axis_impl!(Vec6, x, y, z, w, a, b)
 vec_cast_impl!(Vec6, x, y, z, w, a, b)
 indexable_impl!(Vec6, 6)
+at_fast_impl!(Vec6, 6)
 new_repeat_impl!(Vec6, val, x, y, z, w, a, b)
 dim_impl!(Vec6, 6)
 container_impl!(Vec6)
diff --git a/src/vec_macros.rs b/src/vec_macros.rs
index 3cbb723c..78acbdd0 100644
--- a/src/vec_macros.rs
+++ b/src/vec_macros.rs
@@ -15,6 +15,25 @@ macro_rules! new_impl(
     )
 )
 
+macro_rules! at_fast_impl(
+    ($t: ident, $dim: expr) => (
+        impl<N: Clone> $t<N> {
+            /// Unsafe read access to a vector element by index.
+            #[inline]
+            pub unsafe fn at_fast(&self, i: uint) -> N {
+                (*cast::transmute::<&$t<N>, &[N, ..$dim]>(self)
+                 .unsafe_ref(i)).clone()
+            }
+
+            /// Unsafe write access to a vector element by index.
+            #[inline]
+            pub unsafe fn set_fast(&mut self, i: uint, val: N) {
+                (*cast::transmute::<&mut $t<N>, &mut [N, ..$dim]>(self).unsafe_mut_ref(i)) = val
+            }
+        }
+    )
+)
+
 macro_rules! ord_impl(
     ($t: ident, $comp0: ident $(,$compN: ident)*) => (
         impl<N: Ord> Ord for $t<N> {