diff --git a/src/basics/basics.hh b/src/basics/basics.hh
index 838bbcd..96b8aa2 100644
--- a/src/basics/basics.hh
+++ b/src/basics/basics.hh
@@ -9,8 +9,8 @@
 #ifndef __BASICS_BASICS_HH__
 #define __BASICS_BASICS_HH__
 
-#include "matrix.hh"
-#include "types.hh"
+#include "basics/matrix.hh"
+#include "basics/types.hh"
 
 
 namespace charles {
diff --git a/src/basics/matrix.hh b/src/basics/matrix.hh
index 88ad95b..84c53a2 100644
--- a/src/basics/matrix.hh
+++ b/src/basics/matrix.hh
@@ -18,27 +18,42 @@ namespace basics {
 
 /**
  * A generic, templated Matrix class taking two template parameters. `N` is the
- * number of rows. `M` is the number of columns.
+ * number of rows. `M` is the number of columns. If `M` is not specified, the
+ * matrix will be square.
  */
-template<uint N, uint M>
+template<uint N, uint M = N>
 struct Matrix
 {
     /** Construct an N x M matrix of zeros. */
     static Matrix<N,M> Zero();
 
-    /** Construct an N x M identity matrix. */
+    /** 
+     * Construct an N x M identity matrix. Identity matrices are always square.
+     * It is a (compile time) error to call Identity on a Matrix class where
+     * N != M.
+     */
     static Matrix<N,M> Identity();
 
+    Matrix();
+    Matrix(const Double data[N*M]);
+    Matrix(const Matrix<N,M>& rhs);
+
+    Matrix<N,M>& operator=(const Matrix<N,M>& rhs);
+
+    bool operator==(const Matrix<N,M>& rhs);
+    bool operator!=(const Matrix<N,M>& rhs);
+
     /** Value accessor. Get the ij'th item. */
     Double& operator(uint i, uint j);
 
     /** Scalar multiplication */
-    Matrix<N,M> operator*(const Double& lhs) const;
+    Matrix<N,M> operator*(const Double& rhs) const;
 
     /** Matrix multiplication */
     template<uint P>
-    Matrix<N,P> operator*(Matrix<M,P> lhs) const;
+    Matrix<N,P> operator*(Matrix<M,P> rhs) const;
 
+    /** Get the underlying C array */
     const Double* CArray() const;
 
 private:
@@ -47,11 +62,83 @@ private:
 };
 
 
+typedef Matrix<4> Matrix4;
+
+
 /** Scalar multiplication, scalar factor on the left. */
 template<uint N, uint M>
 Matrix<N,M> operator*(const Double& lhs, const Matrix<N,M>& rhs);
 
 
+/*
+ * charles::basics::Matrix<>::Matrix --
+ */
+template<uint N, uint M>
+Matrix<N,M>::Matrix()
+    : mData()
+{ }
+
+
+/*
+ * charles::basics::Matrix<>::Matrix --
+ */
+template<uint N, uint M>
+Matrix<N,M>::Matrix(const Double data[N*M])
+{
+    memcpy(mData, data, sizeof(Double) * N * M);
+}
+
+
+/*
+ * charles::basics::Matrix<>::Matrix --
+ */
+template<uint N, uint M>
+Matrix<N,M>::Matrix(const Matrix<N,M>& rhs)
+    : Matrix(rhs.mData)
+{ }
+
+
+/*
+ * charles::basics::Matrix<>::operator= --
+ */
+template<uint N, uint M>
+Matrix<N,M>&
+Matrix<N,M>::operator=(const Matrix<N,M>& rhs)
+{
+    memcpy(mData, rhs.mData, sizeof(Double) * N * M);
+    return *this;
+}
+
+
+/*
+ * charles::basics::Matrix<>::operator== --
+ */
+template<uint N, uint M>
+bool
+Matrix<N,M>::operator==(const Matrix<N,M>& rhs)
+    const
+{
+    for (int i = 0; i < N*M; i++) {
+        if (mData[i] != rhs.mData[i]) {
+            return false;
+        }
+    }
+    return true;
+}
+
+
+/*
+ * charles::basics::Matrix<>::operator!= --
+ */
+template<uint N, uint M>
+bool
+Matrix<N,M>::operator!=(const Matrix<N,M>& rhs)
+    const
+{
+    return !(*this == rhs);
+}
+
+
 /*
  * charles::basics::Matrix<>::Zero --
  */
@@ -103,12 +190,12 @@ Matrix<N,M>::operator()(uint i, uint j)
  */
 template<uint N, uint M>
 Matrix<N,M>
-Matrix<N,M>::operator*(const Double& lhs)
+Matrix<N,M>::operator*(const Double& rhs)
     const
 {
     Matrix<N,M> result;
     for (int i = 0; i < N*M; i++) {
-        result.mData = mData[i] * lhs;
+        result.mData = mData[i] * rhs;
     }
     return result;
 }
@@ -120,7 +207,7 @@ Matrix<N,M>::operator*(const Double& lhs)
 template<uint N, uint M>
 template<uint P>
 Matrix<N,P>
-Matrix<N,M>::operator*(Matrix<M,P> lhs)
+Matrix<N,M>::operator*(Matrix<M,P> rhs)
     const
 {
     Matrix<N,P> result;