From a56d1e339aa7693e92eeb1ada6b25bc65fdbceb9 Mon Sep 17 00:00:00 2001
From: Danny Robson <danny@nerdcruft.net>
Date: Tue, 21 Feb 2017 21:19:28 +1100
Subject: [PATCH] matrix: parameterise on rows+cols+type

---
 matrix.cpp  | 225 ++++++++++++++++++++--------------------------------
 matrix.hpp  | 147 +++++++++++++++++-----------------
 matrix.ipp  | 203 ++++++++++++++++++++++++++---------------------
 matrix2.cpp |  28 +++----
 matrix3.cpp |  30 +++----
 matrix4.cpp |  28 +++----
 6 files changed, 319 insertions(+), 342 deletions(-)

diff --git a/matrix.cpp b/matrix.cpp
index 0bbcc8b0..dea7a429 100644
--- a/matrix.cpp
+++ b/matrix.cpp
@@ -23,41 +23,13 @@
 #include <cstring>
 #include <cmath>
 
-using namespace util;
+using util::matrix;
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-matrix<S,T>
-matrix<S,T>::transposed (void) const
-{
-    matrix<S,T> m;
-
-    for (size_t i = 0; i < S; ++i)
-        for (size_t j = 0; j < S; ++j)
-            m.values[i][j] = values[j][i];
-
-    return m;
-}
-
-
-//-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix<S,T>&
-matrix<S,T>::transpose (void)
-{
-    for (size_t i = 0; i < S; ++i)
-        for (size_t j = i + 1; j < S; ++j)
-            std::swap (values[i][j], values[j][i]);
-
-    return *this;
-}
-
-
-//-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix<S,T>&
-matrix<S,T>::invert (void)
+template <size_t Rows, size_t Cols, typename T>
+matrix<Rows,Cols,T>&
+matrix<Rows,Cols,T>::invert (void)
 {
     return *this = inverse ();
 }
@@ -113,32 +85,32 @@ matrix<S,T>::invert (void)
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 T
-util::matrix<S,T>::determinant (void) const
+util::matrix<Rows,Cols,T>::determinant (void) const
 {
     return util::determinant (*this);
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-util::matrix<S,T>
-util::matrix<S,T>::inverse (void) const
+template <size_t Rows, size_t Cols, typename T>
+util::matrix<Rows,Cols,T>
+util::matrix<Rows,Cols,T>::inverse (void) const
 {
     return util::inverse (*this);
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-matrix<S,T>
-util::transposed (const matrix<S,T> &m)
+template <size_t Rows, size_t Cols, typename T>
+matrix<Cols,Rows,T>
+util::transposed (const matrix<Rows,Cols,T> &m)
 {
-    util::matrix<S,T> res;
+    util::matrix<Cols,Rows,T> res;
 
-    for (size_t y = 0; y < S; ++y)
-        for (size_t x = 0; x < S; ++x)
+    for (size_t y = 0; y < Rows; ++y)
+        for (size_t x = 0; x < Cols; ++x)
             res[y][x] = m[x][y];
 
     return res;
@@ -146,47 +118,18 @@ util::transposed (const matrix<S,T> &m)
 
 
 //-----------------------------------------------------------------------------
-template matrix3f util::transposed (const matrix3f&);
-template matrix4f util::transposed (const matrix4f&);
+template util::matrix3f util::transposed (const matrix3f&);
+template util::matrix4f util::transposed (const matrix4f&);
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-matrix<S,T>
-matrix<S,T>::operator* (const matrix<S,T> &rhs) const
+template <size_t Rows, size_t Cols, typename T>
+util::vector<Rows,T>
+matrix<Rows,Cols,T>::operator* (const vector<Rows,T> &rhs) const
 {
-    matrix<S,T> m;
+    vector<Rows,T> out;
 
-    for (unsigned row = 0; row < S; ++row) {
-        for (unsigned col = 0; col < S; ++col) {
-            m.values[row][col] = T {0};
-
-            for (unsigned inner = 0; inner < S; ++inner)
-                m.values[row][col] += values[row][inner] * rhs.values[inner][col];
-        }
-    }
-
-    return m;
-}
-
-
-//-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix<S,T>&
-matrix<S,T>::operator*=(const matrix<S,T> &rhs)
-{
-    return *this = *this * rhs;
-}
-
-
-///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-vector<S,T>
-matrix<S,T>::operator* (const vector<S,T> &rhs) const
-{
-    vector<S,T> out;
-
-    for (size_t i = 0; i < S; ++i)
+    for (size_t i = 0; i < Rows; ++i)
         out[i] = dot (rhs, values[i]);
 
     return out;
@@ -194,13 +137,13 @@ matrix<S,T>::operator* (const vector<S,T> &rhs) const
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-point<S,T>
-matrix<S,T>::operator* (const point<S,T> &rhs) const
+template <size_t Rows, size_t Cols, typename T>
+util::point<Rows,T>
+matrix<Rows,Cols,T>::operator* (const point<Rows,T> &rhs) const
 {
-    point<S,T> out;
+    point<Rows,T> out;
 
-    for (size_t i = 0; i < S; ++i)
+    for (size_t i = 0; i < Rows; ++i)
         out[i] = dot (rhs, values[i]);
 
     return out;
@@ -208,24 +151,27 @@ matrix<S,T>::operator* (const point<S,T> &rhs) const
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 bool
-matrix<S,T>::is_affine (void) const
+matrix<Rows,Cols,T>::is_affine (void) const
 {
-    for (size_t i = 0; i < S - 1; ++i)
-        if (!exactly_zero (values[S-1][i]))
+    if (Rows != Cols)
+        return false;
+
+    for (size_t i = 0; i < Rows - 1; ++i)
+        if (!exactly_zero (values[Rows-1][i]))
             return false;
 
-    return exactly_equal (values[S-1][S-1], T{1});
+    return exactly_equal (values[Rows-1][Rows-1], T{1});
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix4<T>
-matrix<S,T>::ortho (T left, T right,
-                    T bottom, T top,
-                    T near, T far)
+template <size_t Rows, size_t Cols, typename T>
+util::matrix4<T>
+matrix<Rows,Cols,T>::ortho (T left,   T right,
+                            T bottom, T top,
+                            T near,   T far)
 {
     CHECK_GT (far, near);
 
@@ -247,19 +193,19 @@ matrix<S,T>::ortho (T left, T right,
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix4<T>
-matrix<S,T>::ortho2D (T left, T right,
-                      T bottom, T top)
+template <size_t Rows, size_t Cols, typename T>
+util::matrix4<T>
+matrix<Rows, Cols,T>::ortho2D (T left  , T right,
+                               T bottom, T top)
 {
     return ortho (left, right, bottom, top, -1, 1);
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix4<T>
-matrix<S,T>::perspective (T fov, T aspect, range<T> Z)
+template <size_t Rows, size_t Cols, typename T>
+util::matrix4<T>
+matrix<Rows,Cols,T>::perspective (T fov, T aspect, range<T> Z)
 {
     CHECK_GE (Z.lo, 0);
     CHECK_GE (Z.hi, 0);
@@ -288,11 +234,11 @@ matrix<S,T>::perspective (T fov, T aspect, range<T> Z)
 // Implemented for right handed world coordinates.
 //
 // Assumes 'up' is normalised.
-template <size_t S, typename T>
-matrix4<T>
-matrix<S,T>::look_at (const util::point<3,T> eye,
-                      const util::point<3,T> centre,
-                      const util::vector<3,T> up)
+template <size_t Rows, size_t Cols, typename T>
+util::matrix4<T>
+matrix<Rows,Cols,T>::look_at (const util::point<3,T> eye,
+                              const util::point<3,T> centre,
+                              const util::vector<3,T> up)
 {
     CHECK (is_normalised (up));
 
@@ -312,18 +258,18 @@ matrix<S,T>::look_at (const util::point<3,T> eye,
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix4<T>
-matrix<S,T>::translation (util::vector<2,T> v)
+template <size_t Rows, size_t Cols, typename T>
+util::matrix4<T>
+matrix<Rows,Cols,T>::translation (util::vector<2,T> v)
 {
     return translation ({v.x, v.y, 0});
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix4<T>
-matrix<S,T>::translation (util::vector<3,T> v)
+template <size_t Rows, size_t Cols, typename T>
+util::matrix4<T>
+matrix<Rows,Cols,T>::translation (util::vector<3,T> v)
 {
     return { {
         { 1.f, 0.f, 0.f, v.x },
@@ -335,17 +281,17 @@ matrix<S,T>::translation (util::vector<3,T> v)
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix4<T>
-matrix<S,T>::scale (T mag)
+template <size_t Rows, size_t Cols, typename T>
+util::matrix4<T>
+matrix<Rows,Cols,T>::scale (T mag)
 {
     return scale (vector<3,T> (mag));
 }
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix4<T>
-matrix<S,T>::scale (util::vector<3,T> v)
+template <size_t Rows, size_t Cols, typename T>
+util::matrix4<T>
+matrix<Rows,Cols,T>::scale (util::vector<3,T> v)
 {
     return { {
         { v.x, 0.f, 0.f, 0.f },
@@ -357,9 +303,9 @@ matrix<S,T>::scale (util::vector<3,T> v)
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix4<T>
-matrix<S,T>::rotation (T angle, util::vector<3,T> about)
+template <size_t Rows, size_t Cols, typename T>
+util::matrix4<T>
+matrix<Rows,Cols,T>::rotation (T angle, util::vector<3,T> about)
 {
     CHECK (is_normalised (about));
 
@@ -394,25 +340,26 @@ matrix<S,T>::rotation (T angle, util::vector<3,T> about)
 
 
 //-----------------------------------------------------------------------------
-template struct util::matrix<2,float>;
-template struct util::matrix<2,double>;
+template struct util::matrix<2,2,float>;
+template struct util::matrix<2,2,double>;
 
-template struct util::matrix<3,float>;
-template struct util::matrix<3,double>;
+template struct util::matrix<3,3,float>;
+template struct util::matrix<3,3,double>;
 
-template struct util::matrix<4,float>;
-template struct util::matrix<4,double>;
+template struct util::matrix<4,4,float>;
+template struct util::matrix<4,4,double>;
 
 
 ///////////////////////////////////////////////////////////////////////////////
 // Uses the algorithm from:
 //    "Extracting Euler Angles from a Rotation Matrix" by
 //    Mike Day, Insomniac Games.
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 util::vector<3,T>
-util::to_euler (const matrix<S,T> &m)
+util::to_euler (const matrix<Rows,Cols,T> &m)
 {
-    static_assert (S == 3 || S == 4, "only defined for 3d affine transforms");
+    static_assert (Rows == Cols && (Rows == 3 || Rows == 4),
+                   "only defined for 3d affine transforms");
 
     const auto theta0 = std::atan2 (m[2][1], m[2][2]);
 
@@ -431,21 +378,21 @@ util::to_euler (const matrix<S,T> &m)
 
 
 //-----------------------------------------------------------------------------
-template util::vector<3,float> util::to_euler (const matrix<3,float>&);
-template util::vector<3,float> util::to_euler (const matrix<4,float>&);
+template util::vector<3,float> util::to_euler (const matrix<3,3,float>&);
+template util::vector<3,float> util::to_euler (const matrix<4,4,float>&);
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 std::ostream&
-util::operator<< (std::ostream &os, const matrix<S,T> &m)
+util::operator<< (std::ostream &os, const matrix<Rows,Cols,T> &m)
 {
     os << "{ ";
 
-    for (size_t i = 0; i < S; ++i) {
+    for (size_t i = 0; i < Rows; ++i) {
         os << "{ ";
-        std::copy (m[i], m[i]+S, util::infix_iterator<float> (os, ", "));
-        os << ((i == S - 1) ? " }" : " }, ");
+        std::copy_n (m[i], Cols, util::infix_iterator<float> (os, ", "));
+        os << ((i == Rows - 1) ? " }" : " }, ");
     }
 
     return os << " }";
@@ -453,5 +400,5 @@ util::operator<< (std::ostream &os, const matrix<S,T> &m)
 
 
 //-----------------------------------------------------------------------------
-template std::ostream& util::operator<< (std::ostream&, const matrix<4,float>&);
-template std::ostream& util::operator<< (std::ostream&, const matrix<4,double>&);
+template std::ostream& util::operator<< (std::ostream&, const matrix<4,4,float>&);
+template std::ostream& util::operator<< (std::ostream&, const matrix<4,4,double>&);
diff --git a/matrix.hpp b/matrix.hpp
index 47e0bca1..05463b5b 100644
--- a/matrix.hpp
+++ b/matrix.hpp
@@ -21,14 +21,15 @@
 #include "range.hpp"
 
 #include <ostream>
+#include <cstdlib>
 
 namespace util {
-    template <size_t S, typename T>
+    template <size_t Rows, size_t Cols, typename T>
     struct matrix {
-        T values[S][S];
+        static constexpr auto rows = Rows;
+        static constexpr auto cols = Cols;
 
-        static constexpr size_t rows = S;
-        static constexpr size_t cols = S;
+        T values[Rows][Cols];
 
         // index operators return a pointer into the data array so that
         // multidimensional array syntax can be used transparently on this
@@ -46,9 +47,6 @@ namespace util {
         const T* cbegin (void) const;
         const T* cend   (void) const;
 
-        matrix& transpose  (void);
-        matrix  transposed (void) const;
-
         T determinant (void) const;
 
         matrix  inverse (void) const;
@@ -56,29 +54,27 @@ namespace util {
         matrix  inverse_affine (void) const;
         matrix& invert_affine  (void);
 
-        matrix   operator* (const matrix&) const;
-        matrix&  operator*=(const matrix&);
-
-        vector<S,T> operator* (const vector<S,T>&) const;
-        point<S,T>  operator* (const point<S,T> &) const;
+        vector<Rows,T> operator* (const vector<Rows,T>&) const;
+        point<Rows,T>  operator* (const point<Rows,T> &) const;
 
         bool is_affine (void) const;
 
         template <typename U>
-        matrix<S,U> cast (void) const;
+        matrix<Rows,Cols,U>
+        cast (void) const;
 
         // Perspective matrices
-        static matrix<4,T> ortho   (T left, T right, T bottom, T top, T near, T far);
-        static matrix<4,T> ortho2D (T left, T right, T bottom, T top);
-        static matrix<4,T> perspective (T fov, T aspect, range<T> Z);
-        static matrix<4,T> look_at (point<3,T> eye, point<3,T> target, vector<3,T> up);
+        static matrix<4,4,T> ortho   (T left, T right, T bottom, T top, T near, T far);
+        static matrix<4,4,T> ortho2D (T left, T right, T bottom, T top);
+        static matrix<4,4,T> perspective (T fov, T aspect, range<T> Z);
+        static matrix<4,4,T> look_at (point<3,T> eye, point<3,T> target, vector<3,T> up);
 
         // Affine matrices
-        static matrix<4,T> translation (util::vector<2,T>);
-        static matrix<4,T> translation (util::vector<3,T>);
-        static matrix<4,T> scale       (util::vector<3,T>);
-        static matrix<4,T> scale       (T);
-        static matrix<4,T> rotation    (T angle, util::vector<3,T> about);
+        static matrix<4,4,T> translation (util::vector<2,T>);
+        static matrix<4,4,T> translation (util::vector<3,T>);
+        static matrix<4,4,T> scale       (util::vector<3,T>);
+        static matrix<4,4,T> scale       (T);
+        static matrix<4,4,T> rotation    (T angle, util::vector<3,T> about);
 
         // Constant matrices
         static constexpr matrix identity ();
@@ -91,101 +87,110 @@ namespace util {
     //
     // Results are undefined if the matrix is not purely a rotation matrix,
     // or if the dimension is not 3x3 or 4x4.
-    template <size_t S, typename T>
+    template <size_t Rows, size_t Cols, typename T>
     vector<3,T>
-    to_euler (const matrix<S,T>&);
+    to_euler (const matrix<Rows, Cols, T>&);
 
 
     ///////////////////////////////////////////////////////////////////////////
     // logical operations
-    template <size_t S, typename T>
+    template <size_t Rows, size_t Cols, typename T>
     constexpr
     bool
-    operator== (const matrix<S,T>&, const matrix<S,T>&);
+    operator== (const matrix<Rows,Cols,T>&, const matrix<Rows,Cols,T>&);
 
 
-    template <size_t S, typename T>
+    template <size_t Rows, size_t Cols, typename T>
     constexpr
     bool
-    operator!= (const matrix<S,T>&, const matrix<S,T>&);
+    operator!= (const matrix<Rows,Cols,T>&, const matrix<Rows,Cols,T>&);
 
 
     ///////////////////////////////////////////////////////////////////////////
     // element operations
-    template <size_t S, typename T>
+    template <size_t Rows, size_t Cols, typename T>
     constexpr
-    matrix<S,T>
-    operator+ (const matrix<S,T>&, const matrix<S,T>&);
+    matrix<Rows,Cols,T>
+    operator+ (const matrix<Rows,Cols,T>&, const matrix<Rows,Cols,T>&);
 
-    template <size_t S, typename T>
+    template <size_t Rows, size_t Cols, typename T>
     constexpr
-    matrix<S,T>
-    operator- (const matrix<S,T>&, const matrix<S,T>&);
+    matrix<Rows,Cols,T>
+    operator- (const matrix<Rows,Cols,T>&, const matrix<Rows,Cols,T>&);
 
 
     ///////////////////////////////////////////////////////////////////////////
     // scalar operations
-    template <size_t S, typename T> constexpr matrix<S,T> operator* (const matrix<S,T>&, T);
-    template <size_t S, typename T> constexpr matrix<S,T> operator/ (const matrix<S,T>&, T);
-    template <size_t S, typename T> constexpr matrix<S,T> operator+ (const matrix<S,T>&, T);
-    template <size_t S, typename T> constexpr matrix<S,T> operator- (const matrix<S,T>&, T);
+    template <size_t R, size_t C, typename T> constexpr matrix<R,C,T> operator* (const matrix<R,C,T>&, T);
+    template <size_t R, size_t C, typename T> constexpr matrix<R,C,T> operator/ (const matrix<R,C,T>&, T);
+    template <size_t R, size_t C, typename T> constexpr matrix<R,C,T> operator+ (const matrix<R,C,T>&, T);
+    template <size_t R, size_t C, typename T> constexpr matrix<R,C,T> operator- (const matrix<R,C,T>&, T);
 
-    template <size_t S, typename T> constexpr matrix<S,T> operator* (T, const matrix<S,T>&);
-    template <size_t S, typename T> constexpr matrix<S,T> operator/ (T, const matrix<S,T>&);
-    template <size_t S, typename T> constexpr matrix<S,T> operator+ (T, const matrix<S,T>&);
-    template <size_t S, typename T> constexpr matrix<S,T> operator- (T, const matrix<S,T>&);
+    template <size_t R, size_t C, typename T> constexpr matrix<R,C,T> operator* (T, const matrix<R,C,T>&);
+    template <size_t R, size_t C, typename T> constexpr matrix<R,C,T> operator/ (T, const matrix<R,C,T>&);
+    template <size_t R, size_t C, typename T> constexpr matrix<R,C,T> operator+ (T, const matrix<R,C,T>&);
+    template <size_t R, size_t C, typename T> constexpr matrix<R,C,T> operator- (T, const matrix<R,C,T>&);
 
-    template <size_t S, typename T> constexpr matrix<S,T>& operator*= (matrix<S,T>&, T);
-    template <size_t S, typename T> constexpr matrix<S,T>& operator/= (matrix<S,T>&, T);
-    template <size_t S, typename T> constexpr matrix<S,T>& operator+= (matrix<S,T>&, T);
-    template <size_t S, typename T> constexpr matrix<S,T>& operator-= (matrix<S,T>&, T);
+    template <size_t R, size_t C, typename T> constexpr matrix<R,C,T>& operator*= (matrix<R,C,T>&, T);
+    template <size_t R, size_t C, typename T> constexpr matrix<R,C,T>& operator/= (matrix<R,C,T>&, T);
+    template <size_t R, size_t C, typename T> constexpr matrix<R,C,T>& operator+= (matrix<R,C,T>&, T);
+    template <size_t R, size_t C, typename T> constexpr matrix<R,C,T>& operator-= (matrix<R,C,T>&, T);
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <size_t S, typename T>
+    // matrix operations
+    template <
+        size_t R1, size_t C1,
+        size_t R2, size_t C2,
+        typename T
+    > constexpr matrix<R1,C2,T> operator* (const matrix<R1,C1,T>&, const matrix<R2,C2,T>&);
+
+
+    ///////////////////////////////////////////////////////////////////////////
+    template <size_t Rows, size_t Cols, typename T>
     T
-    determinant (const matrix<S,T>&);
+    determinant (const matrix<Rows,Cols,T>&);
 
-    template <size_t S, typename T>
-    matrix<S,T>
-    inverse (const matrix<S,T>&);
+    template <size_t Rows, size_t Cols, typename T>
+    matrix<Rows,Cols,T>
+    inverse (const matrix<Rows,Cols,T>&);
 
-    template <size_t S, typename T>
-    matrix<S,T>
-    transposed (const matrix<S,T>&);
+    template <size_t Rows, size_t Cols, typename T>
+    matrix<Cols,Rows,T>
+    transposed (const matrix<Rows,Cols,T>&);
 
     ///////////////////////////////////////////////////////////////////////////
 
-    template <size_t S, typename T>
-    matrix<S,T>
-    abs (const matrix<S,T>&);
+    template <size_t Rows, size_t Cols, typename T>
+    matrix<Rows,Cols,T>
+    abs (const matrix<Rows,Cols,T>&);
 
-    template <size_t S, typename T>
+    template <size_t Rows, size_t Cols, typename T>
     constexpr
     T
-    sum (const matrix<S,T>&);
+    sum (const matrix<Rows,Cols,T>&);
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T> using matrix3 = matrix<3,T>;
-    template <typename T> using matrix4 = matrix<4,T>;
+    template <typename T> using matrix3 = matrix<3,3,T>;
+    template <typename T> using matrix4 = matrix<4,4,T>;
 
-    template <size_t S> using matrixf = matrix<S,float>;
-    template <size_t S> using matrixd = matrix<S,double>;
+    template <size_t Rows, size_t Cols> using matrixf = matrix<Rows,Cols,float>;
+    template <size_t Rows, size_t Cols> using matrixd = matrix<Rows,Cols,double>;
 
-    typedef matrix<2,float> matrix2f;
-    typedef matrix<2,double> matrix2d;
+    typedef matrix<2,2,float> matrix2f;
+    typedef matrix<2,2,double> matrix2d;
 
-    typedef matrix<3,float> matrix3f;
-    typedef matrix<3,double> matrix3d;
+    typedef matrix<3,3,float> matrix3f;
+    typedef matrix<3,3,double> matrix3d;
 
-    typedef matrix<4,float> matrix4f;
-    typedef matrix<4,double> matrix4d;
+    typedef matrix<4,4,float> matrix4f;
+    typedef matrix<4,4,double> matrix4d;
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <size_t S, typename T>
-    std::ostream& operator<< (std::ostream&, const matrix<S,T>&);
+    template <size_t Rows, size_t Cols, typename T>
+    std::ostream& operator<< (std::ostream&, const matrix<Rows,Cols,T>&);
 }
 
 #include "matrix.ipp"
diff --git a/matrix.ipp b/matrix.ipp
index 0ab88a2f..1f3f53f8 100644
--- a/matrix.ipp
+++ b/matrix.ipp
@@ -23,102 +23,102 @@
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 T*
-util::matrix<S,T>::operator[] (size_t idx)
+util::matrix<Rows,Cols,T>::operator[] (size_t idx)
 {
     return this->values[idx];
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 const T*
-util::matrix<S,T>::operator[] (size_t idx) const
+util::matrix<Rows,Cols,T>::operator[] (size_t idx) const
 {
     return this->values[idx];
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 T*
-util::matrix<S,T>::data (void)
+util::matrix<Rows,Cols,T>::data (void)
 {
     return begin ();
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 const T*
-util::matrix<S,T>::data (void) const
+util::matrix<Rows,Cols,T>::data (void) const
 {
     return begin ();
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 const T*
-util::matrix<S,T>::begin (void) const
+util::matrix<Rows,Cols,T>::begin (void) const
 {
     return &(*this)[0][0];
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 const T*
-util::matrix<S,T>::end (void) const
+util::matrix<Rows,Cols,T>::end (void) const
 {
-    return &(*this)[S][0];
+    return &(*this)[Rows][0];
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 const T*
-util::matrix<S,T>::cbegin (void) const
+util::matrix<Rows,Cols,T>::cbegin (void) const
 {
     return begin ();
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 const T*
-util::matrix<S,T>::cend (void) const
+util::matrix<Rows,Cols,T>::cend (void) const
 {
     return end ();
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 T*
-util::matrix<S,T>::begin (void)
+util::matrix<Rows,Cols,T>::begin (void)
 {
     return &(*this)[0][0];
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 T*
-util::matrix<S,T>::end (void)
+util::matrix<Rows,Cols,T>::end (void)
 {
-    return &(*this)[S][0];
+    return &(*this)[Rows][0];
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 template <typename U>
-util::matrix<S,U>
-util::matrix<S,T>::cast (void) const
+util::matrix<Rows,Cols,U>
+util::matrix<Rows,Cols,T>::cast (void) const
 {
-    util::matrix<S,U> out;
+    util::matrix<Rows,Cols,T> out;
     std::copy (cbegin (), cend (), std::begin (out));
     return out;
 }
@@ -126,14 +126,14 @@ util::matrix<S,T>::cast (void) const
 
 ///////////////////////////////////////////////////////////////////////////////
 #define MATRIX_ELEMENT_OP(OP)                                       \
-template <size_t S, typename T>                                     \
+template <size_t Rows, size_t Cols, typename T>                     \
 constexpr                                                           \
-util::matrix<S,T>                                                   \
+util::matrix<Rows,Cols,T>                                           \
 util::operator OP (                                                 \
-    const util::matrix<S,T> &a,                                     \
-    const util::matrix<S,T> &b)                                     \
+    const util::matrix<Rows,Cols,T> &a,                             \
+    const util::matrix<Rows,Cols,T> &b)                             \
 {                                                                   \
-    util::matrix<S,T> res {};                                       \
+    util::matrix<Rows,Cols,T> res {};                               \
                                                                     \
     for (size_t i = 0; i < a.rows; ++i)                             \
         for (size_t j = 0; j < a.cols; ++j)                         \
@@ -149,47 +149,47 @@ MATRIX_ELEMENT_OP(+)
 
 
 ///////////////////////////////////////////////////////////////////////////////
-#define MATRIX_SCALAR_OP(OP)                                \
-template <size_t S, typename T>                             \
-constexpr                                                   \
-util::matrix<S,T>                                           \
-util::operator OP (const util::matrix<S,T> &m, const T t)   \
-{                                                           \
-    util::matrix<S,T> res {};                               \
-                                                            \
-    std::transform (                                        \
-        std::cbegin (m),                                    \
-        std::cend   (m),                                    \
-        std::begin  (res),                                  \
-        [&t] (auto x) { return x OP t; }                    \
-    );                                                      \
-                                                            \
-    return res;                                             \
-}                                                           \
-                                                            \
-                                                            \
-template <size_t S, typename T>                             \
-constexpr                                                   \
-util::matrix<S,T>                                           \
-util::operator OP (const T t, const util::matrix<S,T> &m)   \
-{                                                           \
-    return m OP t;                                          \
-}                                                           \
-                                                            \
-                                                            \
-template <size_t S, typename T>                             \
-constexpr                                                   \
-util::matrix<S,T>&                                          \
-util::operator OP##= (util::matrix<S,T> &m, T t)            \
-{                                                           \
-    std::transform (                                        \
-        std::cbegin (m),                                    \
-        std::cend   (m),                                    \
-        std::begin  (m),                                    \
-        [&t] (auto x) { return x OP t; }                    \
-    );                                                      \
-                                                            \
-    return m;                                               \
+#define MATRIX_SCALAR_OP(OP)                                        \
+template <size_t Rows, size_t Cols, typename T>                     \
+constexpr                                                           \
+util::matrix<Rows,Cols,T>                                           \
+util::operator OP (const util::matrix<Rows,Cols,T> &m, const T t)   \
+{                                                                   \
+    util::matrix<Rows,Cols,T> res {};                               \
+                                                                    \
+    std::transform (                                                \
+        std::cbegin (m),                                            \
+        std::cend   (m),                                            \
+        std::begin  (res),                                          \
+        [&t] (auto x) { return x OP t; }                            \
+    );                                                              \
+                                                                    \
+    return res;                                                     \
+}                                                                   \
+                                                                    \
+                                                                    \
+template <size_t Rows, size_t Cols, typename T>                     \
+constexpr                                                           \
+util::matrix<Rows,Cols,T>                                           \
+util::operator OP (const T t, const util::matrix<Rows,Cols,T> &m)   \
+{                                                                   \
+    return m OP t;                                                  \
+}                                                                   \
+                                                                    \
+                                                                    \
+template <size_t Rows, size_t Cols, typename T>                     \
+constexpr                                                           \
+util::matrix<Rows,Cols,T>&                                          \
+util::operator OP##= (util::matrix<Rows,Cols,T> &m, T t)            \
+{                                                                   \
+    std::transform (                                                \
+        std::cbegin (m),                                            \
+        std::cend   (m),                                            \
+        std::begin  (m),                                            \
+        [&t] (auto x) { return x OP t; }                            \
+    );                                                              \
+                                                                    \
+    return m;                                                       \
 }
 
 
@@ -202,64 +202,89 @@ MATRIX_SCALAR_OP(-)
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
+template <
+    std::size_t R1, std::size_t C1,
+    std::size_t R2, std::size_t C2,
+    typename T
+>
 constexpr
-util::matrix<S,T>
-util::matrix<S,T>::zeroes (void)
+util::matrix<R1,C2,T>
+util::operator* (const matrix<R1,C1,T> &a, const matrix<R2,C2,T> &b)
+{
+    static_assert (R2 == C1);
+
+    matrix<R1,C2,T> res {0};
+
+    for (size_t r = 0; r < R1; ++r)
+        for (size_t c = 0; c < C2; ++c)
+            for (size_t i = 0; i < R2; ++i)
+                res[r][c] += a[r][i] * b[i][c];
+
+    return res;
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+template <size_t Rows, size_t Cols, typename T>
+constexpr
+util::matrix<Rows,Cols,T>
+util::matrix<Rows,Cols,T>::zeroes (void)
 {
     return {0};
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 constexpr
-util::matrix<S,T>
-util::matrix<S,T>::identity (void)
+util::matrix<Rows,Cols,T>
+util::matrix<Rows,Cols,T>::identity (void)
 {
+    static_assert (Rows == Cols);
+
     auto m = zeroes ();
-    for (size_t i = 0; i < S; ++i)
+    for (size_t i = 0; i < Rows; ++i)
         m[i][i] = 1;
     return m;
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 constexpr
 bool
-util::operator== (const matrix<S,T> &a, const matrix<S,T> &b)
+util::operator== (const matrix<Rows,Cols,T> &a, const matrix<Rows,Cols,T> &b)
 {
     return std::equal (std::cbegin (a), std::cend (a), std::cbegin (b));
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 constexpr
 bool
-util::operator!= (const matrix<S,T> &a, const matrix<S,T> &b)
+util::operator!= (const matrix<Rows,Cols,T> &a, const matrix<Rows,Cols,T> &b)
 {
     return !(a == b);
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-util::matrix<S,T>
-util::abs (const util::matrix<S,T> &src)
+template <size_t Rows, size_t Cols, typename T>
+util::matrix<Rows,Cols,T>
+util::abs (const util::matrix<Rows,Cols,T> &src)
 {
-    util::matrix<S,T> dst;
+    util::matrix<Rows,Cols,T> dst;
     std::transform (std::cbegin (src), std::cend (src), std::begin (dst), util::abs<T>);
     return dst;
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 constexpr
 T
-util::sum (const util::matrix<S,T> &src)
+util::sum (const util::matrix<Rows, Cols,T> &src)
 {
     return sum (std::cbegin (src), std::cend (src));
 }
diff --git a/matrix2.cpp b/matrix2.cpp
index 4bb7f2f4..393066a4 100644
--- a/matrix2.cpp
+++ b/matrix2.cpp
@@ -20,27 +20,27 @@ using util::matrix;
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 T
-util::determinant (const matrix<S,T> &m)
+util::determinant (const matrix<Rows,Cols,T> &m)
 {
-    static_assert (S == 2, "partial specialization for 2 dimensions");
+    static_assert (Rows == 2 && Cols == 2);
 
     return m[0][0] * m[1][1] - m[0][1] * m[1][0];
 }
 
-template float  util::determinant (const matrix<2,float>&);
-template double util::determinant (const matrix<2,double>&);
+template float  util::determinant (const matrix<2,2,float>&);
+template double util::determinant (const matrix<2,2,double>&);
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-matrix<S,T>
-util::inverse (const matrix<S,T> &m)
+template <size_t Rows, size_t Cols, typename T>
+matrix<Rows,Cols,T>
+util::inverse (const matrix<Rows,Cols,T> &m)
 {
-    static_assert (S == 2, "partial specialization for 2 dimensions");
+    static_assert (Rows == 2 && Cols == 2);
 
-    return matrix<S,T> {
+    return matrix<2,2,T> {
          m[1][1], -m[0][1],
         -m[1][0],  m[0][0]
     } / determinant (m);
@@ -48,10 +48,10 @@ util::inverse (const matrix<S,T> &m)
 
 
 //-----------------------------------------------------------------------------
-template util::matrix<2,float>  util::inverse (const matrix<2,float>&);
-template util::matrix<2,double> util::inverse (const matrix<2,double>&);
+template util::matrix<2,2,float>  util::inverse (const matrix<2,2,float>&);
+template util::matrix<2,2,double> util::inverse (const matrix<2,2,double>&);
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template struct util::matrix<2,float>;
-template struct util::matrix<2,double>;
+template struct util::matrix<2,2,float>;
+template struct util::matrix<2,2,double>;
diff --git a/matrix3.cpp b/matrix3.cpp
index 4e9d232c..64b0b678 100644
--- a/matrix3.cpp
+++ b/matrix3.cpp
@@ -20,29 +20,29 @@ using util::matrix;
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 T
-util::determinant (const matrix<S,T>& m)
+util::determinant (const matrix<Rows,Cols,T>& m)
 {
-    static_assert (S == 3, "hard coded 3x3 specialisation");
+    static_assert (Rows == 3 && Cols == 3);
 
     return m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2]) -
            m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) +
            m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
 }
 
-template float  util::determinant (const matrix<3,float>&);
-template double util::determinant (const matrix<3,double>&);
+template float  util::determinant (const matrix<3,3,float>&);
+template double util::determinant (const matrix<3,3,double>&);
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix<S,T>
-util::inverse (const matrix<S,T> &m)
+template <size_t Rows, size_t Cols, typename T>
+matrix<Rows,Cols,T>
+util::inverse (const matrix<Rows,Cols,T> &m)
 {
-    static_assert (S == 3, "hard coded 3x3 specialisation");
-    
-    return matrix<S,T> {
+    static_assert (Rows == 3 && Cols == 3);
+
+    return matrix<3,3,T> {
         m[1][1] * m[2][2] - m[2][1] * m[1][2],
         m[0][2] * m[2][1] - m[0][1] * m[2][2],
         m[0][1] * m[1][2] - m[0][2] * m[1][1],
@@ -57,10 +57,10 @@ util::inverse (const matrix<S,T> &m)
     } / determinant (m);
 }
 
-template util::matrix<3,float>  util::inverse (const matrix<3,float>&);
-template util::matrix<3,double> util::inverse (const matrix<3,double>&);
+template util::matrix<3,3,float>  util::inverse (const matrix<3,3,float>&);
+template util::matrix<3,3,double> util::inverse (const matrix<3,3,double>&);
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template struct util::matrix<3,float>;
-template struct util::matrix<3,double>;
+template struct util::matrix<3,3,float>;
+template struct util::matrix<3,3,double>;
diff --git a/matrix4.cpp b/matrix4.cpp
index efa1b9d5..bbaca1be 100644
--- a/matrix4.cpp
+++ b/matrix4.cpp
@@ -20,11 +20,11 @@ using util::matrix;
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
+template <size_t Rows, size_t Cols, typename T>
 T
-util::determinant (const matrix<S,T> &m)
+util::determinant (const matrix<Rows,Cols,T> &m)
 {
-    static_assert (S == 4, "hard coded 4x4 specialisation");
+    static_assert (Rows == 4 && Cols == 4);
 
     return m[0][3] * m[1][2] * m[2][1] * m[3][0] - m[0][2] * m[1][3] * m[2][1] * m[3][0] -
            m[0][3] * m[1][1] * m[2][2] * m[3][0] + m[0][1] * m[1][3] * m[2][2] * m[3][0] +
@@ -40,18 +40,18 @@ util::determinant (const matrix<S,T> &m)
            m[0][1] * m[1][0] * m[2][2] * m[3][3] + m[0][0] * m[1][1] * m[2][2] * m[3][3];
 }
 
-template float  util::determinant (const matrix<4,float>&);
-template double util::determinant (const matrix<4,double>&);
+template float  util::determinant (const matrix<4,4,float>&);
+template double util::determinant (const matrix<4,4,double>&);
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-matrix<S,T>
-util::inverse (const matrix<S,T> &m)
+template <size_t Rows, size_t Cols, typename T>
+matrix<Rows,Cols,T>
+util::inverse (const matrix<Rows,Cols,T> &m)
 {
-    static_assert (S == 4, "hard coded 4x4 specialisation");
+    static_assert (Rows == 4 && Cols == 4);
 
-    return matrix<S,T> {
+    return matrix<4,4,T> {
         m[1][2]*m[2][3]*m[3][1] - m[1][3]*m[2][2]*m[3][1] + m[1][3]*m[2][1]*m[3][2] - m[1][1]*m[2][3]*m[3][2] - m[1][2]*m[2][1]*m[3][3] + m[1][1]*m[2][2]*m[3][3],
         m[0][3]*m[2][2]*m[3][1] - m[0][2]*m[2][3]*m[3][1] - m[0][3]*m[2][1]*m[3][2] + m[0][1]*m[2][3]*m[3][2] + m[0][2]*m[2][1]*m[3][3] - m[0][1]*m[2][2]*m[3][3],
         m[0][2]*m[1][3]*m[3][1] - m[0][3]*m[1][2]*m[3][1] + m[0][3]*m[1][1]*m[3][2] - m[0][1]*m[1][3]*m[3][2] - m[0][2]*m[1][1]*m[3][3] + m[0][1]*m[1][2]*m[3][3],
@@ -71,10 +71,10 @@ util::inverse (const matrix<S,T> &m)
     } / determinant (m);
 }
 
-template util::matrix<4,float> util::inverse  (const matrix<4,float>&);
-template util::matrix<4,double> util::inverse (const matrix<4,double>&);
+template util::matrix<4,4,float>  util::inverse (const matrix<4,4,float>&);
+template util::matrix<4,4,double> util::inverse (const matrix<4,4,double>&);
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template struct util::matrix<4,float>;
-template struct util::matrix<4,double>;
+template struct util::matrix<4,4,float>;
+template struct util::matrix<4,4,double>;