diff --git a/util/Matrix.cc b/util/Matrix.cc
index 5a46fe7..aaa0d98 100644
--- a/util/Matrix.cc
+++ b/util/Matrix.cc
@@ -1,6 +1,8 @@
 
 #include "Matrix.h"
 #include <math.h>       /* fabs() */
+#include <iostream>
+#include <iomanip>      /* setprecision() */
 
 #define FP_EQ(x,y) (fabs((x)-(y)) < 0.00001)
 
@@ -209,3 +211,21 @@ bool operator==(const Matrix & m1, const Matrix & m2)
     }
     return true;
 }
+
+std::ostream & operator<<(std::ostream & out, const Matrix & m)
+{
+    out << std::setprecision(3);
+    for (int i = 0; i < 4; i++)
+    {
+        out << (i == 0 ? "[ " : "  ");
+        out << "[ ";
+        for (int j = 0; j < 4; j++)
+        {
+            out << m[i][j];
+            if (j < 3)
+                out << ", ";
+        }
+        out << "]" << std::endl;
+    }
+    return out;
+}
diff --git a/util/Matrix.h b/util/Matrix.h
index 3c0bf80..5b42019 100644
--- a/util/Matrix.h
+++ b/util/Matrix.h
@@ -3,6 +3,7 @@
 #define MATRIX_H MATRIX_H
 
 #include "Vector.h"
+#include <iostream>
 
 class Matrix
 {
@@ -25,6 +26,7 @@ class Matrix
 Matrix operator*(const Matrix & m1, const Matrix & m2);
 Vector operator*(const Matrix & m, const Vector & v);
 bool operator==(const Matrix & m1, const Matrix & m2);
+std::ostream & operator<<(std::ostream & out, const Matrix & m);
 
 #endif