diff --git a/ext/cvmat.cpp b/ext/cvmat.cpp
index 89e031b..7ddfb8d 100644
--- a/ext/cvmat.cpp
+++ b/ext/cvmat.cpp
@@ -233,7 +233,8 @@ void define_ruby_class()
   rb_define_method(rb_klass, "sub", RUBY_METHOD_FUNC(rb_sub), -1);
   rb_define_alias(rb_klass, "-", "sub");
   rb_define_method(rb_klass, "mul", RUBY_METHOD_FUNC(rb_mul), -1);
-  rb_define_alias(rb_klass, "*", "mul");
+  rb_define_method(rb_klass, "mat_mul", RUBY_METHOD_FUNC(rb_mat_mul), -1);
+  rb_define_alias(rb_klass, "*", "mat_mul");
   rb_define_method(rb_klass, "div", RUBY_METHOD_FUNC(rb_div), -1);
   rb_define_alias(rb_klass, "/", "div");
   rb_define_method(rb_klass, "and", RUBY_METHOD_FUNC(rb_and), -1);
@@ -1693,7 +1694,7 @@ rb_mul(int argc, VALUE *argv, VALUE self)
   if (rb_scan_args(argc, argv, "11", &val, &scale) < 2)
     scale = rb_float_new(1.0);
   dest = new_object(cvGetSize(CVARR(self)), cvGetElemType(CVARR(self)));
-  if (rb_obj_is_kind_of(val, rb_klass)) {        
+  if (rb_obj_is_kind_of(val, rb_klass)) {
     cvMul(CVARR(self), CVARR(val), CVARR(dest), NUM2DBL(scale));
   }else{
     CvScalar scl = VALUE_TO_CVSCALAR(val);
@@ -1704,6 +1705,28 @@ rb_mul(int argc, VALUE *argv, VALUE self)
   return dest;
 }
 
+/*
+ * call-seq:
+ *   mat_mul(<i>val[,shiftvec]</i>) -> cvmat
+ * Performs matrix multiplication
+ *   dst = src1 * src2 + shiftvec
+ * <i>val</i> and <i>shiftvec</i> should be CvMat
+ * All the matrices should have the same data type and coordinated sizes. 
+ * Real or complex floating-point matrices are supported.
+ */
+VALUE
+rb_mat_mul(int argc, VALUE *argv, VALUE self)
+{
+  VALUE val, shiftvec, dest;
+  rb_scan_args(argc, argv, "11", &val, &shiftvec);
+  dest = new_object(CVMAT(self)->rows, CVMAT(val)->cols, cvGetElemType(CVARR(self)));
+  if (NIL_P(shiftvec))
+    cvMatMul(CVARR(self), CVARR(val), CVARR(dest));
+  else
+    cvMatMulAdd(CVARR(self), CVARR(val), CVARR(shiftvec), CVARR(dest));
+  return dest;
+}
+
 /*
  * call-seq:
  *   div(<i>val[,scale = 1.0]</i>) -> cvmat
diff --git a/ext/cvmat.h b/ext/cvmat.h
index 02c865f..2857446 100644
--- a/ext/cvmat.h
+++ b/ext/cvmat.h
@@ -93,6 +93,7 @@ VALUE rb_convert_scale_abs(VALUE self, VALUE hash);
 VALUE rb_add(int argc, VALUE *argv, VALUE self);
 VALUE rb_sub(int argc, VALUE *argv, VALUE self);    
 VALUE rb_mul(int argc, VALUE *argv, VALUE self);
+VALUE rb_mat_mul(int argc, VALUE *argv, VALUE self);
 VALUE rb_div(int argc, VALUE *argv, VALUE self);    
 VALUE rb_and(int argc, VALUE *argv, VALUE self);
 VALUE rb_or(int argc, VALUE *argv, VALUE self);
diff --git a/test/test_cvmat.rb b/test/test_cvmat.rb
index d3dd9d9..5295036 100755
--- a/test/test_cvmat.rb
+++ b/test/test_cvmat.rb
@@ -1041,7 +1041,47 @@ class TestCvMat < OpenCVTestCase
   end
 
   def test_mat_mul
-    flunk('FIXME: CvMat#mat_mul is not implemented yet.')
+    m0 = create_cvmat(3, 3, :cv32f, 1) { |j, i, c|
+      CvScalar.new(c * 0.1)
+    }
+    m1 = create_cvmat(3, 3, :cv32f, 1) { |j, i, c|
+      CvScalar.new(c)
+    }
+    m2 = create_cvmat(3, 3, :cv32f, 1) { |j, i, c|
+      CvScalar.new(c + 1)
+    }
+
+    m3 = m0.mat_mul(m1)
+    m4 = m0 * m1
+
+    [m3, m4].each { |m|
+      assert_equal(m1.width, m.width)
+      assert_equal(m1.height, m.height)
+      assert_in_delta(1.5, m[0, 0][0], 0.001)
+      assert_in_delta(1.8, m[0, 1][0], 0.001)
+      assert_in_delta(2.1, m[0, 2][0], 0.001)
+      assert_in_delta(4.2, m[1, 0][0], 0.001)
+      assert_in_delta(5.4, m[1, 1][0], 0.001)
+      assert_in_delta(6.6, m[1, 2][0], 0.001)
+      assert_in_delta(6.9, m[2, 0][0], 0.001)
+      assert_in_delta(9, m[2, 1][0], 0.001)
+      assert_in_delta(11.1, m[2, 2][0], 0.001)
+    }
+
+    m5 = m0.mat_mul(m1, m2)
+    [m5].each { |m|
+      assert_equal(m1.width, m.width)
+      assert_equal(m1.height, m.height)
+      assert_in_delta(2.5, m[0, 0][0], 0.001)
+      assert_in_delta(3.8, m[0, 1][0], 0.001)
+      assert_in_delta(5.1, m[0, 2][0], 0.001)
+      assert_in_delta(8.2, m[1, 0][0], 0.001)
+      assert_in_delta(10.4, m[1, 1][0], 0.001)
+      assert_in_delta(12.6, m[1, 2][0], 0.001)
+      assert_in_delta(13.9, m[2, 0][0], 0.001)
+      assert_in_delta(17, m[2, 1][0], 0.001)
+      assert_in_delta(20.1, m[2, 2][0], 0.001)
+    }
   end
 
   def test_div