diff --git a/src/mat.rs b/src/mat.rs
index 6fcd9ac..644cf35 100644
--- a/src/mat.rs
+++ b/src/mat.rs
@@ -98,6 +98,7 @@ pub trait NumericMatrixNxN<T, ColRow>: MatrixNxN<T, ColRow>,
     static pure fn identity() -> self;
     
     pure fn mul_m(other: &self) -> self;
+    pure fn dot(other: &self) -> T;
     
     pure fn det() -> T;
     pure fn trace() -> T;
@@ -262,6 +263,10 @@ pub impl<T:Copy Num NumCast DefaultEq> Mat2<T>: NumericMatrixNxN<T, Vec2<T>> {
         Mat2::new(self.row(0).dot(&other.col(0)), self.row(1).dot(&other.col(0)),
                   self.row(0).dot(&other.col(1)), self.row(1).dot(&other.col(1)))
     }
+
+    pure fn dot(other: &Mat2<T>) -> T {
+        other.transpose().mul_m(&self).trace()
+    }
     
     pure fn det() -> T {
        self[0][0] * self[1][1] - self[1][0] * self[0][1]
@@ -509,6 +514,10 @@ pub impl<T:Copy Num NumCast DefaultEq> Mat3<T>: NumericMatrixNxN<T, Vec3<T>> {
                   self.row(0).dot(&other.col(2)), self.row(1).dot(&other.col(2)), self.row(2).dot(&other.col(2)))
     }
     
+    pure fn dot(other: &Mat3<T>) -> T {
+        other.transpose().mul_m(&self).trace()
+    }
+
     pure fn det() -> T {
         self.col(0).dot(&self.col(1).cross(&self.col(2)))
     }
@@ -837,6 +846,10 @@ pub impl<T:Copy Num NumCast DefaultEq Signed Ord> Mat4<T>: NumericMatrixNxN<T, V
                   self.row(0).dot(&other.col(3)), self.row(1).dot(&other.col(3)), self.row(2).dot(&other.col(3)), self.row(3).dot(&other.col(3)))
     }
     
+    pure fn dot(other: &Mat4<T>) -> T {
+        other.transpose().mul_m(&self).trace()
+    }
+
     pure fn det() -> T {
         self[0][0]*Mat3::new(self[1][1], self[2][1], self[3][1],
                              self[1][2], self[2][2], self[3][2],
diff --git a/src/test/test_mat.rs b/src/test/test_mat.rs
index 1892618..9a1abae 100644
--- a/src/test/test_mat.rs
+++ b/src/test/test_mat.rs
@@ -47,6 +47,7 @@ fn test_Mat2() {
                                     -1f, -1f);
     assert a.mul_m(&b) == Mat2::new(10.0, 22.0,
                                     13.0, 29.0);
+    assert a.dot(&b) == 40f; 
     
     assert a.transpose() == Mat2::new(1f, 2f,
                                       3f, 4f);
@@ -157,6 +158,7 @@ fn test_Mat3() {
     assert a.mul_m(&b) == Mat3::new(36f,  81f, 126f,
                                     42f,  96f, 150f,
                                     48f, 111f, 174f);
+    assert a.dot(&b) == 330f;
     
     assert a.transpose() == Mat3::new(1f, 2f, 3f,
                                       4f, 5f, 6f,
@@ -297,6 +299,7 @@ fn test_Mat4() {
                                     110f, 254f, 398f, 542f,
                                     120f, 280f, 440f, 600f,
                                     130f, 306f, 482f, 658f);
+    assert a.dot(&b) == 1632f;
     
     assert a.transpose() == Mat4::new( 1f,  2f,  3f,  4f,
                                        5f,  6f,  7f,  8f,