From 286a620021fbba518e9085dfb6f92cf70a2d7faa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marc=20W=C3=A4ckerlin?= <marc@mrw.sh>
Date: Thu, 25 Aug 2016 19:53:49 +0000
Subject: [PATCH] documentation updated

---
 configure.ac      |   4 +-
 debian/control.in |   2 +-
 src/matrix.hxx    | 129 ++++---
 test/basic.cxx    | 858 ----------------------------------------------
 test/makefile.am  |   4 +-
 test/matricxx.cxx | 845 ++++++++++++++++++++++++++++++++++++++++++++-
 6 files changed, 926 insertions(+), 916 deletions(-)
 delete mode 100644 test/basic.cxx

diff --git a/configure.ac b/configure.ac
index 6ed4f19..5c67173 100644
--- a/configure.ac
+++ b/configure.ac
@@ -10,8 +10,8 @@
 
 m4_define(x_package_name, libmatricxx) # project's name
 m4_define(x_major, 1) # project's major version
-m4_define(x_minor, 0) # project's minor version
-m4_define(x_least_diff, 12) # start at 0
+m4_define(x_minor, 2) # project's minor version
+m4_define(x_least_diff, 16) # start at 0
 m4_include(ax_init_standard_project.m4)
 AC_INIT(x_package_name, x_version, x_bugreport, x_package_name)
 AM_INIT_AUTOMAKE([1.9 tar-pax parallel-tests color-tests])
diff --git a/debian/control.in b/debian/control.in
index d8d0ba2..55c937f 100644
--- a/debian/control.in
+++ b/debian/control.in
@@ -1,7 +1,7 @@
 Source: @PACKAGE_NAME@
 Priority: extra
 Maintainer: @AUTHOR@
-Build-Depends: debhelper, subversion, pkg-config, automake, libtool, autotools-dev, lsb-release , doxygen, graphviz, mscgen, libcppunit-dev, g++
+Build-Depends: debhelper, subversion, pkg-config, automake, libtool, autotools-dev, lsb-release , doxygen, graphviz, mscgen, libcppunit-dev, g++, pandoc
 
 Package: @PACKAGE_NAME@
 Section: libs
diff --git a/src/matrix.hxx b/src/matrix.hxx
index 5c7073e..a12c049 100644
--- a/src/matrix.hxx
+++ b/src/matrix.hxx
@@ -22,35 +22,50 @@
 
     */
 
-/// Mathematical Classes and Functions
-/** */
+/// Auxiliary Mathematical Functions
 namespace math {
 
-  template<typename TYPE>
-      bool almostEqual(TYPE a, TYPE b) {
-    a = std::fabs(a);
-    b = std::fabs(b);
-    TYPE diff(std::fabs(a-b));
-    TYPE max(a>b?a:b);
-    if (max<1) return diff<=1000*std::numeric_limits<TYPE>::epsilon();
-    return diff<=max*1000*std::numeric_limits<TYPE>::epsilon();
-  }
-
+    /// Compare Floating Points Whether They Are Almost Equal
+    /** Floating points such as @c float and @c double are not 100%
+        exact, because the numbers are represented by a limited number
+        of bits. That's why floating points should not be compared
+        with normal equality operator @c ==, but use function
+        math::equal. This function detects floating points and then
+        calls almostEqual instead of @c ==. */
+    template<typename TYPE>
+        bool almostEqual(TYPE a, TYPE b) {
+      a = std::fabs(a);
+      b = std::fabs(b);
+      TYPE diff(std::fabs(a-b));
+      TYPE max(a>b?a:b);
+      if (max<1) return diff<=1000*std::numeric_limits<TYPE>::epsilon();
+      return diff<=max*1000*std::numeric_limits<TYPE>::epsilon();
+    }
+
+  /// Check Two Values For Equality
+  /** If the values are floating point variables, it calls
+      math::aux::almostEqual. */
   template<typename TYPE>
       bool equal(const TYPE& a, const TYPE& b) {
     return a==b;
   }
 
+  /// Check if Two <code>long double</code> Values are Nearly Equal
+  /** calls math::aux::almostEqual. */
   template<>
       bool equal(const long double& a, const long double& b) {
     return almostEqual(a, b);
   }
-
+  
+  /// Check if Two @c double Values are Nearly Equal
+  /** calls math::aux::almostEqual. */
   template<>
       bool equal(const double& a, const double& b) {
     return almostEqual(a, b);
   }
-
+  
+  /// Check if Two @c float Values are Nearly Equal
+  /** calls math::aux::almostEqual. */
   template<>
       bool equal(const float& a, const float& b) {
     return almostEqual(a, b);
@@ -58,6 +73,8 @@ namespace math {
 
 }
 
+/** Base class with common functions for Matrix and
+    Matrix<TYPE,0,0>. Implements generic common methods. */
 template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
 
     //..............................................................variables
@@ -72,7 +89,7 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
   public:
 
     /// @name Auxiliary Classes
-    ///{
+    ///@{
 
     /// Return One Row as Vector, internally used for element access
     /** Only used to access values:
@@ -83,6 +100,7 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
         @endcode */
     class RowVector {
       public:
+        /// Get Column given a Matrix Row
         TYPE& operator[](size_t column) {
           assert(column<_m.COLUMNS);
           return _v[column];
@@ -98,6 +116,7 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
     /// Same as RowVector, but in a constant environment.
     class ConstRowVector {
       public:
+        /// Get Column given a Matrix Row
         const TYPE& operator[](size_t column) const {
           assert(column<_m.COLUMNS);
           return _v[column];
@@ -110,13 +129,13 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
         const TYPE *_v;
     };
 
-    ///}
+    ///@}
 
     //................................................................methods
   public:
 
     /// @name construction
-    ///{
+    ///@{
 
     MatrixBase(size_t rows, size_t columns):
         ROWS(rows), COLUMNS(columns),
@@ -130,60 +149,84 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
             _c{std::forward<TYPE>(t)...} {
     }
 
-    ///}
+    ///@}
 
     /// @name element access
-    ///{
+    ///@{
 
+    /// Access Matrix Element at Given Row and Column
+    /** You have three possibilities to access an element of a
+        matrix:
+
+        @code
+           Matrix<int,3,3> m;
+           int a21 = m[2][1];    // use bracket operator
+           int b21 = m(2, 1);    // use function operator
+           int c21 = m.at(2, 1); // use at
+        @endcode */
     TYPE& at(size_t row, size_t column) {
       assert(row<ROWS);
       assert(column<COLUMNS);
       return *((TYPE*)_c+row*COLUMNS+column);
     }
 
+    /// Access Matrix Element at Given Row and Column
+    /** @copydoc at */
     const TYPE& at(size_t row, size_t column) const {
       assert(row<ROWS);
       assert(column<COLUMNS);
       return *((TYPE*)_c+row*COLUMNS+column);
     }
 
+    /// Access Matrix Element at Given Row and Column
+    /** @copydoc at */
     TYPE& operator()(size_t row, size_t column) {
       return at(row, column);
     }
 
+    /// Access Matrix Element at Given Row and Column
+    /** @copydoc at */
     const TYPE& operator()(size_t row, size_t column) const {
       return at(row, column);
     }
 
+    /// Access Matrix Element at Given Row and Column
+    /** @copydoc at */
     RowVector operator[](size_t row) {
       assert(row<ROWS);
       return RowVector(*this, (TYPE*)_c+row*COLUMNS);
     }
 
+    /// Access Matrix Element at Given Row and Column
+    /** @copydoc at */
     const ConstRowVector operator[](size_t row) const {
       assert(row<ROWS);
       return ConstRowVector(*this, (TYPE*)_c+row*COLUMNS);
     }
 
+    /// Get Number Of Rows
     size_t rows() const {
       return ROWS;
     }
 
+    /// Get Number Of Columns
     size_t columns() const {
       return COLUMNS;
     }
 
-    ///}
+    ///@}
 
     /// @name operators
-    ///{
+    ///@{
 
+    /// Assign Other Matrix Of Same Size
     MatrixBase& operator=(const MatrixBase& o) {
       assert_check(o);
       memcpy(_c, o._c, MEM_SIZE);
       return *this;
     }
 
+    /// Compare To Other Matrix
     bool operator==(const MatrixBase& o) const {
       if (!check(o)) return false;
       TYPE *to((TYPE*)(_c)+SIZE), *from((TYPE*)(o._c)+SIZE);
@@ -191,10 +234,12 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
       return true;
     }
 
+    /// Compare To Other Matrix
     bool operator!=(const MatrixBase& o) const {
       return !operator==(o);
     }
 
+    /// Add Other Matrix
     MatrixBase& operator+=(const MatrixBase& o) {
       assert_check(o);
       TYPE *to((TYPE*)(_c)+SIZE), *from((TYPE*)(o._c)+SIZE);
@@ -202,6 +247,7 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
       return *this;
     }
 
+    /// Subtract Other Matrix
     MatrixBase& operator-=(const MatrixBase& o) {
       assert_check(o);
       TYPE *to((TYPE*)(_c)+SIZE), *from((TYPE*)(o._c)+SIZE);
@@ -209,35 +255,42 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
       return *this;
     }
 
+    /// Multiply Matrix With Scalar
     MatrixBase& operator*=(const TYPE& o) {
       TYPE *res((TYPE*)(_c)+SIZE);
       while (res>(TYPE*)(_c)) *--res *= o;
       return *this;
     }
 
+    /// Divide Matrix By Scalar
     MatrixBase& operator/=(const TYPE& o) {
       TYPE *res((TYPE*)(_c)+SIZE);
       while (res>(TYPE*)(_c)) *--res /= o;
       return *this;
     }
 
-    ///}
+    ///@}
 
     /// @name special operations
-    ///{
+    ///@{
 
+    /// Apply Any External Function To Each Element
     MatrixBase& apply(std::function<void(TYPE&)> fn) {
       TYPE *to((TYPE*)(_c)+SIZE);
       while (to>(TYPE*)(_c)) fn(*--to);
       return *this;
     }
 
+    /// Calculate Determinant Of The Matrix
+    /** The Matrix is replaced by it's gaussian representation. */
     TYPE det() {
       TYPE res(gauss());
       for (TYPE *p((TYPE*)(_c)+SIZE); --p>=(TYPE*)(_c); p-=COLUMNS) res *= *p;
       return res;
     }
 
+    /// Calculate Gaussian Representation
+    /** The Matrix is replaced by it's gaussian representation. */
     TYPE gauss() {
       /// calculate using gauss algorithmus
       /// @see http://www.mathebibel.de/determinante-berechnen-nach-gauss
@@ -263,8 +316,6 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
       return lambda;
     }
 
-    ///}
-
     //................................................................methods
   protected:
 
@@ -294,7 +345,7 @@ template<typename TYPE, size_t TROWS=0, size_t TCOLUMNS=0> class Matrix:
   public:
 
     /// @name construction
-    ///{
+    ///@{
 
     Matrix(): Parent(TROWS, TCOLUMNS) {
       memset(Parent::_c, 0, Parent::MEM_SIZE);
@@ -309,10 +360,10 @@ template<typename TYPE, size_t TROWS=0, size_t TCOLUMNS=0> class Matrix:
       static_assert(sizeof...(t)==TROWS*TCOLUMNS, "wrong array size");
     }
 
-    ///}
+    ///@}
 
     /// @name operators
-    ///{
+    ///@{
 
     Matrix& operator=(const Matrix& o) {
       Parent::operator=(o);
@@ -356,10 +407,10 @@ template<typename TYPE, size_t TROWS=0, size_t TCOLUMNS=0> class Matrix:
       return res;
     }
 
-    ///}
+    ///@}
 
     /// @name special operations
-    ///{
+    ///@{
 
     Matrix& apply(std::function<void(TYPE&)> fn) {
       Parent::apply(fn);
@@ -428,7 +479,7 @@ template<typename TYPE, size_t TROWS=0, size_t TCOLUMNS=0> class Matrix:
       return *this;
     }
 
-    ///}
+    ///@}
 };
 
 //==============================================================================
@@ -444,7 +495,7 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
   public:
 
     /// @name construction
-    ///{
+    ///@{
 
     Matrix() = delete;
 
@@ -467,19 +518,19 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
       copy_args(Parent::_c, t...);
     }
 
-    ///}
+    ///@}
 
     /// @name destruction
-    ///{
+    ///@{
 
     virtual ~Matrix() {
       delete[] Parent::_c;
     }
 
-    ///}
+    ///@}
 
     /// @name operators
-    ///{
+    ///@{
 
     Matrix& operator=(const Matrix& o) {
       Parent::operator=(o);
@@ -522,10 +573,10 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
       return res;
     }
 
-    ///}
+    ///@}
 
     ///@name special operations
-    ///{
+    ///@{
 
     Matrix& resize(size_t rows, size_t columns) {
       if (rows!=this->ROWS||columns!=this->COLUMNS) {
@@ -608,7 +659,7 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
       return *this;
     }
 
-    ///}
+    ///@}
 
     //................................................................methods
   protected:
diff --git a/test/basic.cxx b/test/basic.cxx
deleted file mode 100644
index a6233a6..0000000
--- a/test/basic.cxx
+++ /dev/null
@@ -1,858 +0,0 @@
-/*! @file
-
-    @id $Id$
-*/
-//       1         2         3         4         5         6         7         8
-// 45678901234567890123456789012345678901234567890123456789012345678901234567890
-
-#include <matrix.hxx>
-#include <cppunit/TestFixture.h>
-#include <cppunit/ui/text/TestRunner.h>
-#include <cppunit/extensions/HelperMacros.h>
-#include <cppunit/extensions/TestFactoryRegistry.h>
-#include <cppunit/XmlOutputter.h>
-#include <fstream>
-
-class TemplateMatrixTest: public CppUnit::TestFixture { 
-  public:
-    template<typename T>
-        void initFromArray1() {
-      const Matrix<T,2,4> m {1, 2, 3, 4,
-                             5, 6, 7, 8};
-      Matrix<T,2,4> m2;
-      m2 = {1, 2, 3, 4,
-            5, 6, 7, 8};
-      CPPUNIT_ASSERT_EQUAL(m, m2);
-      CPPUNIT_ASSERT_EQUAL((T)1, m[0][0]);
-      CPPUNIT_ASSERT_EQUAL((T)2, m[0][1]);
-      CPPUNIT_ASSERT_EQUAL((T)3, m[0][2]);
-      CPPUNIT_ASSERT_EQUAL((T)4, m[0][3]);
-      CPPUNIT_ASSERT_EQUAL((T)5, m[1][0]);
-      CPPUNIT_ASSERT_EQUAL((T)6, m[1][1]);
-      CPPUNIT_ASSERT_EQUAL((T)7, m[1][2]);
-      CPPUNIT_ASSERT_EQUAL((T)8, m[1][3]);
-    }
-    template<typename T>
-        void initFromArray2() {
-      const Matrix<T,2,4> m(1, 2, 3, 4,
-                            5, 6, 7, 8);
-      CPPUNIT_ASSERT_EQUAL((T)1, m[0][0]);
-      CPPUNIT_ASSERT_EQUAL((T)2, m[0][1]);
-      CPPUNIT_ASSERT_EQUAL((T)3, m[0][2]);
-      CPPUNIT_ASSERT_EQUAL((T)4, m[0][3]);
-      CPPUNIT_ASSERT_EQUAL((T)5, m[1][0]);
-      CPPUNIT_ASSERT_EQUAL((T)6, m[1][1]);
-      CPPUNIT_ASSERT_EQUAL((T)7, m[1][2]);
-      CPPUNIT_ASSERT_EQUAL((T)8, m[1][3]);
-    }
-    template<typename T>
-        void initFromArray3() {
-      const Matrix<T,2,4> m({1, 2, 3, 4,
-                             5, 6, 7, 8});
-      CPPUNIT_ASSERT_EQUAL((T)1, m[0][0]);
-      CPPUNIT_ASSERT_EQUAL((T)2, m[0][1]);
-      CPPUNIT_ASSERT_EQUAL((T)3, m[0][2]);
-      CPPUNIT_ASSERT_EQUAL((T)4, m[0][3]);
-      CPPUNIT_ASSERT_EQUAL((T)5, m[1][0]);
-      CPPUNIT_ASSERT_EQUAL((T)6, m[1][1]);
-      CPPUNIT_ASSERT_EQUAL((T)7, m[1][2]);
-      CPPUNIT_ASSERT_EQUAL((T)8, m[1][3]);
-    }
-    template<typename T>
-        void initFromOtherMatrix() {
-      Matrix<T,2,4> m1(1, 2, 3, 4,
-                       5, 6, 7, 8);
-      m1[1][2] = 13;
-      Matrix<T,2,4> m2(m1);
-      m1[0][2] = 16;
-      m2[0][0] = 0;
-      Matrix<T,2,4> m3;
-      m3 = m2;
-      CPPUNIT_ASSERT_EQUAL(m2, m3);
-      CPPUNIT_ASSERT_EQUAL((T)1, m1[0][0]);
-      CPPUNIT_ASSERT_EQUAL((T)2, m1[0][1]);
-      CPPUNIT_ASSERT_EQUAL((T)16, m1[0][2]);
-      CPPUNIT_ASSERT_EQUAL((T)4, m1[0][3]);
-      CPPUNIT_ASSERT_EQUAL((T)5, m1[1][0]);
-      CPPUNIT_ASSERT_EQUAL((T)6, m1[1][1]);
-      CPPUNIT_ASSERT_EQUAL((T)13, m1[1][2]);
-      CPPUNIT_ASSERT_EQUAL((T)8, m1[1][3]);
-      CPPUNIT_ASSERT_EQUAL((T)0, m2[0][0]);
-      CPPUNIT_ASSERT_EQUAL((T)2, m2[0][1]);
-      CPPUNIT_ASSERT_EQUAL((T)3, m2[0][2]);
-      CPPUNIT_ASSERT_EQUAL((T)4, m2[0][3]);
-      CPPUNIT_ASSERT_EQUAL((T)5, m2[1][0]);
-      CPPUNIT_ASSERT_EQUAL((T)6, m2[1][1]);
-      CPPUNIT_ASSERT_EQUAL((T)13, m2[1][2]);
-      CPPUNIT_ASSERT_EQUAL((T)8, m2[1][3]);
-    }
-    template<typename T>
-        void initFromDefault() {
-      const Matrix<T,2,2> m;
-      CPPUNIT_ASSERT_EQUAL((T)0, m[0][0]);
-      CPPUNIT_ASSERT_EQUAL((T)0, m[0][1]);
-      CPPUNIT_ASSERT_EQUAL((T)0, m[1][0]);
-      CPPUNIT_ASSERT_EQUAL((T)0, m[1][1]);
-    }
-    template<typename T>
-        void access() {
-      Matrix<T,2,4> m1(1, 2, 3, 4,
-                       5, 6, 7, 8);
-      for (size_t row(0); row<2; ++row)
-        for (size_t column(0); column<4; ++column)
-          CPPUNIT_ASSERT_EQUAL(m1(row, column), m1[row][column]);
-    }
-    template<typename T>
-        void equality() {
-      Matrix<T,2,4> m1(1, 2, 3, 4,
-                       5, 6, 7, 8);
-      Matrix<T,2,4> m2(1, 2, 3, 4,
-                       5, 6, 7, 8);
-      Matrix<T,2,4> m3(1, 2, 3, 4,
-                       5, 6, 7, 9);
-      Matrix<T,2,4> m4(9, 2, 3, 4,
-                       5, 6, 7, 8);
-      Matrix<T,2,4> m5(1, 2, 0, 4,
-                       5, 6, 7, 8);
-      CPPUNIT_ASSERT(m1==m2);
-      CPPUNIT_ASSERT(m1!=m3);
-      CPPUNIT_ASSERT(m1!=m4);
-      CPPUNIT_ASSERT(m1!=m5);
-    }
-    template<typename T>
-        void operator_plus() {
-      const Matrix<T,2,4> m1(1, 2, 3, 4,
-                             5, 6, 7, 8);
-      const Matrix<T,2,4> m2(2, 4, 6, 8,
-                             1, 3, 5, 7);
-      const Matrix<T,2,4> m(m1+m2);
-      const Matrix<T,2,4> res(3, 6, 9, 12,
-                              6, 9, 12, 15);
-      CPPUNIT_ASSERT_EQUAL(res, m);
-    }
-    template<typename T>
-        void operator_minus() {
-      const Matrix<T,2,4> m1(1, 2, 3, 4,
-                             5, 6, 7, 8);
-      const Matrix<T,2,4> m2(2, 4, 6, 8,
-                             1, 3, 5, 7);
-      const Matrix<T,2,4> m(m1-m2);
-      const Matrix<T,2,4> res(-1, -2, -3, -4,
-                              4, 3, 2, 1);
-      CPPUNIT_ASSERT_EQUAL(res, m);
-      CPPUNIT_ASSERT_EQUAL(res, -m2+m1);
-    }
-    template<typename T>
-        void scalar_mult() {
-      const Matrix<T,2,4> m1(1, 2, 3, 4,
-                             5, 6, 7, 8);
-      T two(2), three(3), four(4), big(32);
-      CPPUNIT_ASSERT_EQUAL(m1*two, m1+m1);
-      CPPUNIT_ASSERT_EQUAL(m1*three, m1+m1+m1);
-      CPPUNIT_ASSERT_EQUAL(m1*four, m1+m1+m1+m1);
-      CPPUNIT_ASSERT_EQUAL(two*m1, m1*two);
-      CPPUNIT_ASSERT_EQUAL(three*m1, m1*three);
-      CPPUNIT_ASSERT_EQUAL(four*m1, m1*four);
-      CPPUNIT_ASSERT_EQUAL(two*m1*two, m1*four);
-      CPPUNIT_ASSERT_EQUAL(big*m1*four, m1*four*big);
-    }
-    template<typename T>
-        void scalar_div() {
-      const Matrix<T,2,4> m1(2, 4, 6, 8,
-                             10, 12, 14, 16);
-      const Matrix<T,2,4> m2(1, 2, 3, 4,
-                             5, 6, 7, 8);
-      CPPUNIT_ASSERT_EQUAL(m2, m1/(T)2);
-      CPPUNIT_ASSERT_EQUAL(m2, (T)3*m1/(T)6);
-    }
-    template<typename T>
-        void matrix_mult() {
-      const Matrix<T,2,3> m1(1, 2, 3,
-                             4, 5, 6);
-      const Matrix<T,3,4> m2(1, 4, 7, 10,
-                             2, 5, 8, 11,
-                             3, 6, 9, 12);
-      const Matrix<T,2,4> m(m1*m2);
-      const Matrix<T,2,4> res(14, 32,  50,  68,
-                              32, 77, 122, 167);
-      CPPUNIT_ASSERT_EQUAL(res, m);
-    }
-    template<typename T>
-        void transpose() {
-      const Matrix<T,2,3> m(1, 2, 3,
-                            4, 5, 6);
-      const Matrix<T,3,2> res(1, 4,
-                              2, 5,
-                              3, 6);
-      CPPUNIT_ASSERT_EQUAL(res, m.t());
-    }
-    template<typename T>
-        void apply() {
-      Matrix<T,3,3> m(2, -2, 4,
-                      -2, 1, -6,
-                      1, 0, -2);
-      Matrix<T,3,3> o(m);
-      CPPUNIT_ASSERT_EQUAL((T)3*o, m.apply([](T& t){t*=3;}));
-    }
-    template<typename T>
-        void gauss() {
-      Matrix<T,3,3> m(2, -2, 4,
-                      -2, 1, -6,
-                      1, 0, -2);
-      const Matrix<T,3,3> res(1, -1, 2,
-                              0, -1, -2,
-                              0, 0, -6);
-      T lambda(m.gauss());
-      CPPUNIT_ASSERT_EQUAL(res, m);
-      CPPUNIT_ASSERT_EQUAL((T)2, lambda);
-    }
-    template<typename T>
-        void det() {
-      Matrix<T,3,3> m(2, -2, 4,
-                      -2, 1, -6,
-                      1, 0, -2);
-      CPPUNIT_ASSERT_EQUAL((T)12, m.det());
-    }
-    template<typename T>
-        void i() {
-      const Matrix<T,3,3> m1(1, 0, 0,
-                             0, 1, 0,
-                             0, 0, 1);
-      const Matrix<T,3,2> m2(1, 0,
-                             0, 1,
-                             0, 0);
-      const Matrix<T,2,3> m3(1, 0, 0,
-                             0, 1, 0);
-      CPPUNIT_ASSERT_EQUAL(m1, m1.i());
-      CPPUNIT_ASSERT_EQUAL(m2, m2.i());
-      CPPUNIT_ASSERT_EQUAL(m3, m3.i());
-    }
-    template<typename T>
-        void inv() {
-      {
-        Matrix<T,3,3> m(2, -1, 0,
-                        1, 2, -2,
-                        0, -1, 1);
-        const Matrix<T,3,3> res(0, 1, 2,
-                                -1, 2, 4,
-                                -1, 2, 5);
-        Matrix<T,3,3> o1(m), o2(m);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL((T)2*res, (T)2/o1);
-        CPPUNIT_ASSERT_EQUAL(o1.i(), o1/o2);
-        CPPUNIT_ASSERT_EQUAL(res, m);
-        
-      } {
-        Matrix<T,3,3> m(1, 2, 3,
-                        0, 1, 4,
-                        5, 6, 0);
-        const Matrix<T,3,3> res(-24, 18, 5,
-                                20, -15, -4,
-                                -5, 4, 1);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL(res, m);
-      } {
-        Matrix<T,3,3> m(1, 2, 3,
-                        0, 4, 5,
-                        1, 0, 6);
-        const Matrix<T,3,3> res((T)12/11, (T)-6/11, (T)-1/11,
-                                (T)5/22, (T)3/22, (T)-5/22,
-                                (T)-2/11, (T)1/11, (T)2/11);
-        Matrix<T,3,3> o1(m), o2(m);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL((T)2*res, (T)2/o1);
-        CPPUNIT_ASSERT_EQUAL(o1.i(), o1/o2);
-        CPPUNIT_ASSERT_EQUAL(res, m);
-      } {
-        Matrix<T,3,3> m(-2, 0, 1,
-                        9, 2, -3,
-                        5, 1, -2);
-        const Matrix<T,3,3> res(-1, 1, -2,
-                                3, -1, 3,
-                                -1, 2, -4);
-        Matrix<T,3,3> o(m);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL(m.i(), m*o);
-        CPPUNIT_ASSERT_EQUAL(res, m);
-      } {
-        Matrix<T,3,3> m(1, 3, 1,
-                        1, 1, 2,
-                        2, 3, 4);
-        const Matrix<T,3,3> res(2, 9, -5,
-                                0, -2, 1,
-                                -1, -3, 2);
-        Matrix<T,3,3> o(m);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL(m.i(), m*o);
-        CPPUNIT_ASSERT_EQUAL(res, m);
-      } {
-        Matrix<T,4,4> m(2, 1, 4, 1,
-                        -1, 1, 0, 2,
-                        0, 0, 2, 4,
-                        2, -2, 0, 1);
-        const Matrix<T,4,4> res((T)-1/3, (T)13/15, (T)-2/3, 0.6,
-                                (T)1/3, (T)16/15, (T)-2/3, 0.2,
-                                0, -0.8, 0.5, -0.4,
-                                0, 0.4, 0, 0.2);
-        Matrix<T,4,4> o(m);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL(m.i(), m*o);
-        CPPUNIT_ASSERT_EQUAL(res, m);
-      } {
-        Matrix<T,2,2> m(4, 3,
-                        3, 2);
-        const Matrix<T,2,2> res(-2, 3,
-                                3, -4);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL(res, m);
-      }
-    }
-    template<typename T>
-        void stream() {
-      const Matrix<T,3,4> m1(1, 2, 3, 4,
-                             5, 6, 7, 8,
-                             1, 4, 2, 8);
-      Matrix<T,3,4> m2;
-      std::string res("[3x4]{1,2,3,4,5,6,7,8,1,4,2,8}");
-      std::stringstream ss;
-      ss<<m1;
-      CPPUNIT_ASSERT_EQUAL(res, ss.str());
-      ss>>m2;
-      CPPUNIT_ASSERT_EQUAL(m1, m2);
-    }
-    CPPUNIT_TEST_SUITE(TemplateMatrixTest);
-    CPPUNIT_TEST(initFromArray1<int>);
-    CPPUNIT_TEST(initFromArray1<long>);
-    CPPUNIT_TEST(initFromArray1<unsigned>);
-    CPPUNIT_TEST(initFromArray1<unsigned long>);
-    CPPUNIT_TEST(initFromArray1<float>);
-    CPPUNIT_TEST(initFromArray1<double>);
-    CPPUNIT_TEST(initFromArray2<int>);
-    CPPUNIT_TEST(initFromArray2<long>);
-    CPPUNIT_TEST(initFromArray2<unsigned>);
-    CPPUNIT_TEST(initFromArray2<unsigned long>);
-    CPPUNIT_TEST(initFromArray2<float>);
-    CPPUNIT_TEST(initFromArray2<double>);
-    CPPUNIT_TEST(initFromArray3<int>);
-    CPPUNIT_TEST(initFromArray3<long>);
-    CPPUNIT_TEST(initFromArray3<unsigned>);
-    CPPUNIT_TEST(initFromArray3<unsigned long>);
-    CPPUNIT_TEST(initFromArray3<float>);
-    CPPUNIT_TEST(initFromArray3<double>);
-    CPPUNIT_TEST(initFromOtherMatrix<int>);
-    CPPUNIT_TEST(initFromOtherMatrix<long>);
-    CPPUNIT_TEST(initFromOtherMatrix<unsigned>);
-    CPPUNIT_TEST(initFromOtherMatrix<unsigned long>);
-    CPPUNIT_TEST(initFromOtherMatrix<float>);
-    CPPUNIT_TEST(initFromOtherMatrix<double>);
-    CPPUNIT_TEST(access<int>);
-    CPPUNIT_TEST(access<long>);
-    CPPUNIT_TEST(access<unsigned>);
-    CPPUNIT_TEST(access<unsigned long>);
-    CPPUNIT_TEST(access<float>);
-    CPPUNIT_TEST(access<double>);
-    CPPUNIT_TEST(equality<int>);
-    CPPUNIT_TEST(equality<long>);
-    CPPUNIT_TEST(equality<unsigned>);
-    CPPUNIT_TEST(equality<unsigned long>);
-    CPPUNIT_TEST(equality<float>);
-    CPPUNIT_TEST(equality<double>);
-    CPPUNIT_TEST(operator_plus<int>);
-    CPPUNIT_TEST(operator_plus<long>);
-    CPPUNIT_TEST(operator_plus<unsigned>);
-    CPPUNIT_TEST(operator_plus<unsigned long>);
-    CPPUNIT_TEST(operator_plus<float>);
-    CPPUNIT_TEST(operator_plus<double>);
-    CPPUNIT_TEST(operator_minus<int>);
-    CPPUNIT_TEST(operator_minus<long>);
-    CPPUNIT_TEST(operator_minus<unsigned>);
-    CPPUNIT_TEST(operator_minus<unsigned long>);
-    CPPUNIT_TEST(operator_minus<float>);
-    CPPUNIT_TEST(operator_minus<double>);
-    CPPUNIT_TEST(scalar_mult<int>);
-    CPPUNIT_TEST(scalar_mult<long>);
-    CPPUNIT_TEST(scalar_mult<unsigned>);
-    CPPUNIT_TEST(scalar_mult<unsigned long>);
-    CPPUNIT_TEST(scalar_mult<float>);
-    CPPUNIT_TEST(scalar_mult<double>);
-    CPPUNIT_TEST(scalar_div<int>);
-    CPPUNIT_TEST(scalar_div<long>);
-    CPPUNIT_TEST(scalar_div<unsigned>);
-    CPPUNIT_TEST(scalar_div<unsigned long>);
-    CPPUNIT_TEST(scalar_div<float>);
-    CPPUNIT_TEST(scalar_div<double>);
-    CPPUNIT_TEST(matrix_mult<int>);
-    CPPUNIT_TEST(matrix_mult<long>);
-    CPPUNIT_TEST(matrix_mult<unsigned>);
-    CPPUNIT_TEST(matrix_mult<unsigned long>);
-    CPPUNIT_TEST(matrix_mult<float>);
-    CPPUNIT_TEST(matrix_mult<double>);
-    CPPUNIT_TEST(transpose<int>);
-    CPPUNIT_TEST(transpose<long>);
-    CPPUNIT_TEST(transpose<unsigned>);
-    CPPUNIT_TEST(transpose<unsigned long>);
-    CPPUNIT_TEST(transpose<float>);
-    CPPUNIT_TEST(transpose<double>);
-    CPPUNIT_TEST(apply<int>);
-    CPPUNIT_TEST(apply<long>);
-    CPPUNIT_TEST(apply<unsigned>);
-    CPPUNIT_TEST(apply<unsigned long>);
-    CPPUNIT_TEST(apply<float>);
-    CPPUNIT_TEST(apply<double>);
-    CPPUNIT_TEST(gauss<int>);
-    CPPUNIT_TEST(gauss<long>);
-    CPPUNIT_TEST(gauss<float>);
-    CPPUNIT_TEST(gauss<double>);
-    CPPUNIT_TEST(det<int>);
-    CPPUNIT_TEST(det<long>);
-    CPPUNIT_TEST(det<float>);
-    CPPUNIT_TEST(det<double>);
-    CPPUNIT_TEST(i<int>);
-    CPPUNIT_TEST(i<long>);
-    CPPUNIT_TEST(i<unsigned>);
-    CPPUNIT_TEST(i<unsigned long>);
-    CPPUNIT_TEST(i<float>);
-    CPPUNIT_TEST(i<double>);
-    CPPUNIT_TEST(inv<float>);
-    CPPUNIT_TEST(inv<double>);
-    CPPUNIT_TEST(inv<long double>);
-    CPPUNIT_TEST(stream<int>);
-    CPPUNIT_TEST(stream<long>);
-    CPPUNIT_TEST(stream<unsigned>);
-    CPPUNIT_TEST(stream<unsigned long>);
-    CPPUNIT_TEST(stream<float>);
-    CPPUNIT_TEST(stream<double>);
-    CPPUNIT_TEST_SUITE_END();
-};
-CPPUNIT_TEST_SUITE_REGISTRATION(TemplateMatrixTest);
-
-class VariableMatrixTest: public CppUnit::TestFixture { 
-  public:
-    template<typename T>
-        void initFromArray1() {
-      Matrix<T> m(2,4,
-                  1, 2, 3, 4,
-                  5, 6, 7, 8);
-      Matrix<T> m2(2, 4);
-      m2 = {2, 4,
-            1, 2, 3, 4,
-            5, 6, 7, 8};
-      CPPUNIT_ASSERT_EQUAL(m, m2);
-      CPPUNIT_ASSERT_EQUAL((T)1, m[0][0]);
-      CPPUNIT_ASSERT_EQUAL((T)2, m[0][1]);
-      CPPUNIT_ASSERT_EQUAL((T)3, m[0][2]);
-      CPPUNIT_ASSERT_EQUAL((T)4, m[0][3]);
-      CPPUNIT_ASSERT_EQUAL((T)5, m[1][0]);
-      CPPUNIT_ASSERT_EQUAL((T)6, m[1][1]);
-      CPPUNIT_ASSERT_EQUAL((T)7, m[1][2]);
-      CPPUNIT_ASSERT_EQUAL((T)8, m[1][3]);
-    }
-    template<typename T>
-        void initFromArray2() {
-      Matrix<T> m(2, 4,
-                  1, 2, 3, 4,
-                  5, 6, 7, 8);
-      CPPUNIT_ASSERT_EQUAL((T)1, m[0][0]);
-      CPPUNIT_ASSERT_EQUAL((T)2, m[0][1]);
-      CPPUNIT_ASSERT_EQUAL((T)3, m[0][2]);
-      CPPUNIT_ASSERT_EQUAL((T)4, m[0][3]);
-      CPPUNIT_ASSERT_EQUAL((T)5, m[1][0]);
-      CPPUNIT_ASSERT_EQUAL((T)6, m[1][1]);
-      CPPUNIT_ASSERT_EQUAL((T)7, m[1][2]);
-      CPPUNIT_ASSERT_EQUAL((T)8, m[1][3]);
-    }
-    template<typename T>
-        void initFromOtherMatrix() {
-      Matrix<T> m1(2, 4,
-                   1, 2, 3, 4,
-                   5, 6, 7, 8);
-      m1[1][2] = 13;
-      Matrix<T> m2(m1);
-      m1[0][2] = 16;
-      m2[0][0] = 0;
-      Matrix<T> m3(2, 4);
-      m3 = m2;
-      CPPUNIT_ASSERT_EQUAL(m2, m3);
-      CPPUNIT_ASSERT_EQUAL((T)1, m1[0][0]);
-      CPPUNIT_ASSERT_EQUAL((T)2, m1[0][1]);
-      CPPUNIT_ASSERT_EQUAL((T)16, m1[0][2]);
-      CPPUNIT_ASSERT_EQUAL((T)4, m1[0][3]);
-      CPPUNIT_ASSERT_EQUAL((T)5, m1[1][0]);
-      CPPUNIT_ASSERT_EQUAL((T)6, m1[1][1]);
-      CPPUNIT_ASSERT_EQUAL((T)13, m1[1][2]);
-      CPPUNIT_ASSERT_EQUAL((T)8, m1[1][3]);
-      CPPUNIT_ASSERT_EQUAL((T)0, m2[0][0]);
-      CPPUNIT_ASSERT_EQUAL((T)2, m2[0][1]);
-      CPPUNIT_ASSERT_EQUAL((T)3, m2[0][2]);
-      CPPUNIT_ASSERT_EQUAL((T)4, m2[0][3]);
-      CPPUNIT_ASSERT_EQUAL((T)5, m2[1][0]);
-      CPPUNIT_ASSERT_EQUAL((T)6, m2[1][1]);
-      CPPUNIT_ASSERT_EQUAL((T)13, m2[1][2]);
-      CPPUNIT_ASSERT_EQUAL((T)8, m2[1][3]);
-    }
-    template<typename T>
-        void initFromDefault() {
-      Matrix<T> m(2, 4);
-      CPPUNIT_ASSERT_EQUAL((T)0, m[0][0]);
-      CPPUNIT_ASSERT_EQUAL((T)0, m[0][1]);
-      CPPUNIT_ASSERT_EQUAL((T)0, m[1][0]);
-      CPPUNIT_ASSERT_EQUAL((T)0, m[1][1]);
-    }
-    template<typename T>
-        void access() {
-      Matrix<T> m1(2, 4,
-                   1, 2, 3, 4,
-                   5, 6, 7, 8);
-      for (size_t row(0); row<2; ++row)
-        for (size_t column(0); column<4; ++column)
-          CPPUNIT_ASSERT_EQUAL(m1(row, column), m1[row][column]);
-    }
-    template<typename T>
-        void equality() {
-      Matrix<T> m1(2, 4,
-                   1, 2, 3, 4,
-                   5, 6, 7, 8);
-      Matrix<T> m2(2, 4,
-                   1, 2, 3, 4,
-                   5, 6, 7, 8);
-      Matrix<T> m3(2, 4,
-                   1, 2, 3, 4,
-                   5, 6, 7, 9);
-      Matrix<T> m4(2, 4,
-                   9, 2, 3, 4,
-                   5, 6, 7, 8);
-      CPPUNIT_ASSERT(m1==m2);
-      CPPUNIT_ASSERT(m1!=m3);
-      CPPUNIT_ASSERT(m1!=m4);
-    }
-    template<typename T>
-        void operator_plus() {
-      const Matrix<T> m1(2, 4,
-                         1, 2, 3, 4,
-                         5, 6, 7, 8);
-      const Matrix<T> m2(2, 4,
-                         2, 4, 6, 8,
-                         1, 3, 5, 7);
-      const Matrix<T> m(m1+m2);
-      const Matrix<T> res(2, 4,
-                          3, 6, 9, 12,
-                          6, 9, 12, 15);
-      CPPUNIT_ASSERT_EQUAL(res, m);
-    }
-    template<typename T>
-        void operator_minus() {
-      const Matrix<T> m1(2, 4,
-                         1, 2, 3, 4,
-                         5, 6, 7, 8);
-      const Matrix<T> m2(2, 4,
-                         2, 4, 6, 8,
-                         1, 3, 5, 7);
-      const Matrix<T> m(m1-m2);
-      const Matrix<T> res(2, 4,
-                         -1, -2, -3, -4,
-                          4, 3, 2, 1);
-      CPPUNIT_ASSERT_EQUAL(res, m);
-      CPPUNIT_ASSERT_EQUAL(res, -m2+m1);
-    }
-    template<typename T>
-        void scalar_mult() {
-      const Matrix<T> m1(2, 4,
-                         1, 2, 3, 4,
-                         5, 6, 7, 8);
-      T two(2), three(3), four(4), big(32);
-      CPPUNIT_ASSERT_EQUAL(m1*two, m1+m1);
-      CPPUNIT_ASSERT_EQUAL(m1*three, m1+m1+m1);
-      CPPUNIT_ASSERT_EQUAL(m1*four, m1+m1+m1+m1);
-      CPPUNIT_ASSERT_EQUAL(two*m1, m1*two);
-      CPPUNIT_ASSERT_EQUAL(three*m1, m1*three);
-      CPPUNIT_ASSERT_EQUAL(four*m1, m1*four);
-      CPPUNIT_ASSERT_EQUAL(two*m1*two, m1*four);
-      CPPUNIT_ASSERT_EQUAL(big*m1*four, m1*four*big);
-    }
-    template<typename T>
-        void scalar_div() {
-      const Matrix<T> m1(2, 4,
-                         2, 4, 6, 8,
-                         10, 12, 14, 16);
-      const Matrix<T> m2(2, 4,
-                         1, 2, 3, 4,
-                         5, 6, 7, 8);
-      CPPUNIT_ASSERT_EQUAL(m2, m1/(T)2);
-      CPPUNIT_ASSERT_EQUAL(m2, (T)3*m1/(T)6);
-    }
-    template<typename T>
-        void matrix_mult() {
-      const Matrix<T> m1(2, 3,
-                         1, 2, 3,
-                         4, 5, 6);
-      const Matrix<T> m2(3, 4,
-                         1, 4, 7, 10,
-                         2, 5, 8, 11,
-                         3, 6, 9, 12);
-      const Matrix<T> m(m1*m2);
-      const Matrix<T> res(2, 4,
-                          14, 32,  50,  68,
-                          32, 77, 122, 167);
-      CPPUNIT_ASSERT_EQUAL(res, m);
-    }
-    template<typename T>
-        void transpose() {
-      const Matrix<T> m(2, 3,
-                        1, 2, 3,
-                        4, 5, 6);
-      const Matrix<T> res(3, 2,
-                          1, 4,
-                          2, 5,
-                          3, 6);
-      CPPUNIT_ASSERT_EQUAL(res, m.t());
-    }
-    template<typename T>
-        void apply() {
-          Matrix<T> m(3, 3,
-                      2, -2, 4,
-                      -2, 1, -6,
-                      1, 0, -2);
-          Matrix<T> o(m);
-          CPPUNIT_ASSERT_EQUAL((T)3*o, m.apply([](T& t){t*=3;}));
-    }
-    template<typename T>
-        void gauss() {
-      Matrix<T> m(3, 3,
-                  2, -2, 4,
-                  -2, 1, -6,
-                  1, 0, -2);
-      const Matrix<T> res(3,3,
-                          1, -1, 2,
-                          0, -1, -2,
-                          0, 0, -6);
-      T lambda(m.gauss());
-      CPPUNIT_ASSERT_EQUAL(res, m);
-      CPPUNIT_ASSERT_EQUAL((T)2, lambda);
-    }
-    template<typename T>
-        void det() {
-      Matrix<T> m(3, 3,
-                  2, -2, 4,
-                  -2, 1, -6,
-                  1, 0, -2);
-      CPPUNIT_ASSERT_EQUAL((T)12, m.det());
-    }
-    template<typename T>
-        void i() {
-      const Matrix<T> m1(3, 3,
-                         1, 0, 0,
-                         0, 1, 0,
-                         0, 0, 1);
-      const Matrix<T> m2(3, 2,
-                         1, 0,
-                         0, 1,
-                         0, 0);
-      const Matrix<T> m3(2, 3,
-                         1, 0, 0,
-                         0, 1, 0);
-      CPPUNIT_ASSERT_EQUAL(m1, m1.i());
-      CPPUNIT_ASSERT_EQUAL(m2, m2.i());
-      CPPUNIT_ASSERT_EQUAL(m3, m3.i());
-    }
-    template<typename T>
-        void inv() {
-      {
-        Matrix<T> m(3, 3,
-                    2, -1, 0,
-                    1, 2, -2,
-                    0, -1, 1);
-        const Matrix<T> res(3, 3,
-                            0, 1, 2,
-                            -1, 2, 4,
-                            -1, 2, 5);
-        Matrix<T> o1(m), o2(m);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL((T)2*res, (T)2/o1);
-        CPPUNIT_ASSERT_EQUAL(o1.i(), o1/o2);
-        CPPUNIT_ASSERT_EQUAL(res, m);
-      } {
-        Matrix<T> m(3, 3,
-                    1, 2, 3,
-                    0, 1, 4,
-                    5, 6, 0);
-        const Matrix<T> res(3, 3,
-                            -24, 18, 5,
-                            20, -15, -4,
-                            -5, 4, 1);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL(res, m);
-      } {
-        Matrix<T> m(3, 3,
-                    1, 2, 3,
-                    0, 4, 5,
-                    1, 0, 6);
-        const Matrix<T> res(3, 3,
-                            (T)12/11, (T)-6/11, (T)-1/11,
-                            (T)5/22, (T)3/22, (T)-5/22,
-                            (T)-2/11, (T)1/11, (T)2/11);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL(res, m);
-      } {
-        Matrix<T> m(3, 3,
-                    -2, 0, 1,
-                    9, 2, -3,
-                    5, 1, -2);
-        const Matrix<T> res(3, 3,
-                            -1, 1, -2,
-                            3, -1, 3,
-                            -1, 2, -4);
-        Matrix<T> o(m);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL(m.i(), m*o);
-        CPPUNIT_ASSERT_EQUAL(res, m);
-      } {
-        Matrix<T> m(4, 4,
-                    2, 1, 4, 1,
-                    -1, 1, 0, 2,
-                    0, 0, 2, 4,
-                    2, -2, 0, 1);
-        const Matrix<T> res(4, 4,
-                            (T)-1/3, (T)13/15, (T)-2/3, 0.6,
-                            (T)1/3, (T)16/15, (T)-2/3, 0.2,
-                            0, -0.8, 0.5, -0.4,
-                            0, 0.4, 0, 0.2);
-        Matrix<T> o(m);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL(m.i(), m*o);
-        CPPUNIT_ASSERT_EQUAL(res, m);
-      } {
-        Matrix<T> m(2, 2,
-                    4, 3,
-                    3, 2);
-        const Matrix<T> res(2, 2,
-                            -2, 3,
-                            3, -4);
-        m.inv();
-        CPPUNIT_ASSERT_EQUAL(res, m);
-      }
-    }
-    template<typename T>
-        void stream() {
-      const Matrix<T> m1(3, 4,
-                         1, 2, 3, 4,
-                         5, 6, 7, 8,
-                         1, 4, 2, 8);
-      Matrix<T> m2(1, 1);
-      std::string res("[3x4]{1,2,3,4,5,6,7,8,1,4,2,8}");
-      std::stringstream ss;
-      ss<<m1;
-      CPPUNIT_ASSERT_EQUAL(res, ss.str());
-      ss>>m2;
-      CPPUNIT_ASSERT_EQUAL(m1, m2);
-    }
-    CPPUNIT_TEST_SUITE(VariableMatrixTest);
-    CPPUNIT_TEST(initFromArray1<int>);
-    CPPUNIT_TEST(initFromArray1<long>);
-    CPPUNIT_TEST(initFromArray1<unsigned>);
-    CPPUNIT_TEST(initFromArray1<unsigned long>);
-    CPPUNIT_TEST(initFromArray1<float>);
-    CPPUNIT_TEST(initFromArray1<double>);
-    CPPUNIT_TEST(initFromArray2<int>);
-    CPPUNIT_TEST(initFromArray2<long>);
-    CPPUNIT_TEST(initFromArray2<unsigned>);
-    CPPUNIT_TEST(initFromArray2<unsigned long>);
-    CPPUNIT_TEST(initFromArray2<float>);
-    CPPUNIT_TEST(initFromArray2<double>);
-    CPPUNIT_TEST(initFromOtherMatrix<int>);
-    CPPUNIT_TEST(initFromOtherMatrix<long>);
-    CPPUNIT_TEST(initFromOtherMatrix<unsigned>);
-    CPPUNIT_TEST(initFromOtherMatrix<unsigned long>);
-    CPPUNIT_TEST(initFromOtherMatrix<float>);
-    CPPUNIT_TEST(initFromOtherMatrix<double>);
-    CPPUNIT_TEST(access<int>);
-    CPPUNIT_TEST(access<long>);
-    CPPUNIT_TEST(access<unsigned>);
-    CPPUNIT_TEST(access<unsigned long>);
-    CPPUNIT_TEST(access<float>);
-    CPPUNIT_TEST(access<double>);
-    CPPUNIT_TEST(equality<int>);
-    CPPUNIT_TEST(equality<long>);
-    CPPUNIT_TEST(equality<unsigned>);
-    CPPUNIT_TEST(equality<unsigned long>);
-    CPPUNIT_TEST(equality<float>);
-    CPPUNIT_TEST(equality<double>);
-    CPPUNIT_TEST(operator_plus<int>);
-    CPPUNIT_TEST(operator_plus<long>);
-    CPPUNIT_TEST(operator_plus<unsigned>);
-    CPPUNIT_TEST(operator_plus<unsigned long>);
-    CPPUNIT_TEST(operator_plus<float>);
-    CPPUNIT_TEST(operator_plus<double>);
-    CPPUNIT_TEST(operator_minus<int>);
-    CPPUNIT_TEST(operator_minus<long>);
-    CPPUNIT_TEST(operator_minus<unsigned>);
-    CPPUNIT_TEST(operator_minus<unsigned long>);
-    CPPUNIT_TEST(operator_minus<float>);
-    CPPUNIT_TEST(operator_minus<double>);
-    CPPUNIT_TEST(scalar_mult<int>);
-    CPPUNIT_TEST(scalar_mult<long>);
-    CPPUNIT_TEST(scalar_mult<unsigned>);
-    CPPUNIT_TEST(scalar_mult<unsigned long>);
-    CPPUNIT_TEST(scalar_mult<float>);
-    CPPUNIT_TEST(scalar_mult<double>);
-    CPPUNIT_TEST(scalar_mult<int>);
-    CPPUNIT_TEST(scalar_div<int>);
-    CPPUNIT_TEST(scalar_div<long>);
-    CPPUNIT_TEST(scalar_div<unsigned>);
-    CPPUNIT_TEST(scalar_div<unsigned long>);
-    CPPUNIT_TEST(scalar_div<float>);
-    CPPUNIT_TEST(scalar_div<double>);
-    CPPUNIT_TEST(matrix_mult<int>);
-    CPPUNIT_TEST(matrix_mult<long>);
-    CPPUNIT_TEST(matrix_mult<unsigned>);
-    CPPUNIT_TEST(matrix_mult<unsigned long>);
-    CPPUNIT_TEST(matrix_mult<float>);
-    CPPUNIT_TEST(matrix_mult<double>);
-    CPPUNIT_TEST(transpose<int>);
-    CPPUNIT_TEST(transpose<long>);
-    CPPUNIT_TEST(transpose<unsigned>);
-    CPPUNIT_TEST(transpose<unsigned long>);
-    CPPUNIT_TEST(transpose<float>);
-    CPPUNIT_TEST(transpose<double>);
-    CPPUNIT_TEST(apply<int>);
-    CPPUNIT_TEST(apply<long>);
-    CPPUNIT_TEST(apply<unsigned>);
-    CPPUNIT_TEST(apply<unsigned long>);
-    CPPUNIT_TEST(apply<float>);
-    CPPUNIT_TEST(apply<double>);
-    CPPUNIT_TEST(gauss<int>);
-    CPPUNIT_TEST(gauss<long>);
-    CPPUNIT_TEST(gauss<float>);
-    CPPUNIT_TEST(gauss<double>);
-    CPPUNIT_TEST(det<int>);
-    CPPUNIT_TEST(det<long>);
-    CPPUNIT_TEST(det<float>);
-    CPPUNIT_TEST(det<double>);
-    CPPUNIT_TEST(i<int>);
-    CPPUNIT_TEST(i<long>);
-    CPPUNIT_TEST(i<unsigned>);
-    CPPUNIT_TEST(i<unsigned long>);
-    CPPUNIT_TEST(i<float>);
-    CPPUNIT_TEST(i<double>);
-    CPPUNIT_TEST(inv<float>);
-    CPPUNIT_TEST(inv<double>);
-    CPPUNIT_TEST(inv<long double>);
-    CPPUNIT_TEST(stream<int>);
-    CPPUNIT_TEST(stream<long>);
-    CPPUNIT_TEST(stream<unsigned>);
-    CPPUNIT_TEST(stream<unsigned long>);
-    CPPUNIT_TEST(stream<float>);
-    CPPUNIT_TEST(stream<double>);
-    CPPUNIT_TEST_SUITE_END();
-};
-CPPUNIT_TEST_SUITE_REGISTRATION(VariableMatrixTest);
-
-int main(int argc, char** argv) try {
-  std::ofstream ofs((*argv+std::string(".xml")).c_str());
-  CppUnit::TextUi::TestRunner runner;
-  runner.setOutputter(new CppUnit::XmlOutputter(&runner.result(), ofs));
-  runner.addTest(CppUnit::TestFactoryRegistry::getRegistry().makeTest());
-  return runner.run() ? 0 : 1;
- } catch (std::exception& e) {
-  std::cerr<<"***Exception: "<<e.what()<<std::endl;
-  return 1;
- }
diff --git a/test/makefile.am b/test/makefile.am
index 8904e48..da3c735 100644
--- a/test/makefile.am
+++ b/test/makefile.am
@@ -12,10 +12,10 @@ AM_CPPFLAGS = -I${top_srcdir}/src -I${top_builddir}/src
 AM_LDFLAGS = -L${abs_top_builddir}/src/.libs
 LDADD = -lcppunit -lmatricxx
 
-check_PROGRAMS = basic
+check_PROGRAMS = matricxx
 TESTS = ${check_PROGRAMS}
 
-basic_SOURCES = basic.cxx
+matricxx_SOURCES = matricxx.cxx
 
 CLEANFILES = ${check_PROGRAMS:%=%.xml}
 MAINTAINERCLEANFILES = makefile.in
diff --git a/test/matricxx.cxx b/test/matricxx.cxx
index 2dc5fd0..a6233a6 100644
--- a/test/matricxx.cxx
+++ b/test/matricxx.cxx
@@ -1,15 +1,11 @@
-/** @id $Id$
-
-   This file has been added:
-     - by bootstrap.sh
-     - on Thu, 25 August 2016 18:18:47 +0200
+/*! @file
 
+    @id $Id$
 */
 //       1         2         3         4         5         6         7         8
 // 45678901234567890123456789012345678901234567890123456789012345678901234567890
 
-
-
+#include <matrix.hxx>
 #include <cppunit/TestFixture.h>
 #include <cppunit/ui/text/TestRunner.h>
 #include <cppunit/extensions/HelperMacros.h>
@@ -17,17 +13,838 @@
 #include <cppunit/XmlOutputter.h>
 #include <fstream>
 
-/// @todo Rename DummyTest and DummyTest::dummy()
-/// @todo Write test cases
-class DummyTest: public CppUnit::TestFixture { 
+class TemplateMatrixTest: public CppUnit::TestFixture { 
   public:
-    void dummy() {
+    template<typename T>
+        void initFromArray1() {
+      const Matrix<T,2,4> m {1, 2, 3, 4,
+                             5, 6, 7, 8};
+      Matrix<T,2,4> m2;
+      m2 = {1, 2, 3, 4,
+            5, 6, 7, 8};
+      CPPUNIT_ASSERT_EQUAL(m, m2);
+      CPPUNIT_ASSERT_EQUAL((T)1, m[0][0]);
+      CPPUNIT_ASSERT_EQUAL((T)2, m[0][1]);
+      CPPUNIT_ASSERT_EQUAL((T)3, m[0][2]);
+      CPPUNIT_ASSERT_EQUAL((T)4, m[0][3]);
+      CPPUNIT_ASSERT_EQUAL((T)5, m[1][0]);
+      CPPUNIT_ASSERT_EQUAL((T)6, m[1][1]);
+      CPPUNIT_ASSERT_EQUAL((T)7, m[1][2]);
+      CPPUNIT_ASSERT_EQUAL((T)8, m[1][3]);
+    }
+    template<typename T>
+        void initFromArray2() {
+      const Matrix<T,2,4> m(1, 2, 3, 4,
+                            5, 6, 7, 8);
+      CPPUNIT_ASSERT_EQUAL((T)1, m[0][0]);
+      CPPUNIT_ASSERT_EQUAL((T)2, m[0][1]);
+      CPPUNIT_ASSERT_EQUAL((T)3, m[0][2]);
+      CPPUNIT_ASSERT_EQUAL((T)4, m[0][3]);
+      CPPUNIT_ASSERT_EQUAL((T)5, m[1][0]);
+      CPPUNIT_ASSERT_EQUAL((T)6, m[1][1]);
+      CPPUNIT_ASSERT_EQUAL((T)7, m[1][2]);
+      CPPUNIT_ASSERT_EQUAL((T)8, m[1][3]);
+    }
+    template<typename T>
+        void initFromArray3() {
+      const Matrix<T,2,4> m({1, 2, 3, 4,
+                             5, 6, 7, 8});
+      CPPUNIT_ASSERT_EQUAL((T)1, m[0][0]);
+      CPPUNIT_ASSERT_EQUAL((T)2, m[0][1]);
+      CPPUNIT_ASSERT_EQUAL((T)3, m[0][2]);
+      CPPUNIT_ASSERT_EQUAL((T)4, m[0][3]);
+      CPPUNIT_ASSERT_EQUAL((T)5, m[1][0]);
+      CPPUNIT_ASSERT_EQUAL((T)6, m[1][1]);
+      CPPUNIT_ASSERT_EQUAL((T)7, m[1][2]);
+      CPPUNIT_ASSERT_EQUAL((T)8, m[1][3]);
+    }
+    template<typename T>
+        void initFromOtherMatrix() {
+      Matrix<T,2,4> m1(1, 2, 3, 4,
+                       5, 6, 7, 8);
+      m1[1][2] = 13;
+      Matrix<T,2,4> m2(m1);
+      m1[0][2] = 16;
+      m2[0][0] = 0;
+      Matrix<T,2,4> m3;
+      m3 = m2;
+      CPPUNIT_ASSERT_EQUAL(m2, m3);
+      CPPUNIT_ASSERT_EQUAL((T)1, m1[0][0]);
+      CPPUNIT_ASSERT_EQUAL((T)2, m1[0][1]);
+      CPPUNIT_ASSERT_EQUAL((T)16, m1[0][2]);
+      CPPUNIT_ASSERT_EQUAL((T)4, m1[0][3]);
+      CPPUNIT_ASSERT_EQUAL((T)5, m1[1][0]);
+      CPPUNIT_ASSERT_EQUAL((T)6, m1[1][1]);
+      CPPUNIT_ASSERT_EQUAL((T)13, m1[1][2]);
+      CPPUNIT_ASSERT_EQUAL((T)8, m1[1][3]);
+      CPPUNIT_ASSERT_EQUAL((T)0, m2[0][0]);
+      CPPUNIT_ASSERT_EQUAL((T)2, m2[0][1]);
+      CPPUNIT_ASSERT_EQUAL((T)3, m2[0][2]);
+      CPPUNIT_ASSERT_EQUAL((T)4, m2[0][3]);
+      CPPUNIT_ASSERT_EQUAL((T)5, m2[1][0]);
+      CPPUNIT_ASSERT_EQUAL((T)6, m2[1][1]);
+      CPPUNIT_ASSERT_EQUAL((T)13, m2[1][2]);
+      CPPUNIT_ASSERT_EQUAL((T)8, m2[1][3]);
+    }
+    template<typename T>
+        void initFromDefault() {
+      const Matrix<T,2,2> m;
+      CPPUNIT_ASSERT_EQUAL((T)0, m[0][0]);
+      CPPUNIT_ASSERT_EQUAL((T)0, m[0][1]);
+      CPPUNIT_ASSERT_EQUAL((T)0, m[1][0]);
+      CPPUNIT_ASSERT_EQUAL((T)0, m[1][1]);
+    }
+    template<typename T>
+        void access() {
+      Matrix<T,2,4> m1(1, 2, 3, 4,
+                       5, 6, 7, 8);
+      for (size_t row(0); row<2; ++row)
+        for (size_t column(0); column<4; ++column)
+          CPPUNIT_ASSERT_EQUAL(m1(row, column), m1[row][column]);
+    }
+    template<typename T>
+        void equality() {
+      Matrix<T,2,4> m1(1, 2, 3, 4,
+                       5, 6, 7, 8);
+      Matrix<T,2,4> m2(1, 2, 3, 4,
+                       5, 6, 7, 8);
+      Matrix<T,2,4> m3(1, 2, 3, 4,
+                       5, 6, 7, 9);
+      Matrix<T,2,4> m4(9, 2, 3, 4,
+                       5, 6, 7, 8);
+      Matrix<T,2,4> m5(1, 2, 0, 4,
+                       5, 6, 7, 8);
+      CPPUNIT_ASSERT(m1==m2);
+      CPPUNIT_ASSERT(m1!=m3);
+      CPPUNIT_ASSERT(m1!=m4);
+      CPPUNIT_ASSERT(m1!=m5);
+    }
+    template<typename T>
+        void operator_plus() {
+      const Matrix<T,2,4> m1(1, 2, 3, 4,
+                             5, 6, 7, 8);
+      const Matrix<T,2,4> m2(2, 4, 6, 8,
+                             1, 3, 5, 7);
+      const Matrix<T,2,4> m(m1+m2);
+      const Matrix<T,2,4> res(3, 6, 9, 12,
+                              6, 9, 12, 15);
+      CPPUNIT_ASSERT_EQUAL(res, m);
+    }
+    template<typename T>
+        void operator_minus() {
+      const Matrix<T,2,4> m1(1, 2, 3, 4,
+                             5, 6, 7, 8);
+      const Matrix<T,2,4> m2(2, 4, 6, 8,
+                             1, 3, 5, 7);
+      const Matrix<T,2,4> m(m1-m2);
+      const Matrix<T,2,4> res(-1, -2, -3, -4,
+                              4, 3, 2, 1);
+      CPPUNIT_ASSERT_EQUAL(res, m);
+      CPPUNIT_ASSERT_EQUAL(res, -m2+m1);
+    }
+    template<typename T>
+        void scalar_mult() {
+      const Matrix<T,2,4> m1(1, 2, 3, 4,
+                             5, 6, 7, 8);
+      T two(2), three(3), four(4), big(32);
+      CPPUNIT_ASSERT_EQUAL(m1*two, m1+m1);
+      CPPUNIT_ASSERT_EQUAL(m1*three, m1+m1+m1);
+      CPPUNIT_ASSERT_EQUAL(m1*four, m1+m1+m1+m1);
+      CPPUNIT_ASSERT_EQUAL(two*m1, m1*two);
+      CPPUNIT_ASSERT_EQUAL(three*m1, m1*three);
+      CPPUNIT_ASSERT_EQUAL(four*m1, m1*four);
+      CPPUNIT_ASSERT_EQUAL(two*m1*two, m1*four);
+      CPPUNIT_ASSERT_EQUAL(big*m1*four, m1*four*big);
+    }
+    template<typename T>
+        void scalar_div() {
+      const Matrix<T,2,4> m1(2, 4, 6, 8,
+                             10, 12, 14, 16);
+      const Matrix<T,2,4> m2(1, 2, 3, 4,
+                             5, 6, 7, 8);
+      CPPUNIT_ASSERT_EQUAL(m2, m1/(T)2);
+      CPPUNIT_ASSERT_EQUAL(m2, (T)3*m1/(T)6);
+    }
+    template<typename T>
+        void matrix_mult() {
+      const Matrix<T,2,3> m1(1, 2, 3,
+                             4, 5, 6);
+      const Matrix<T,3,4> m2(1, 4, 7, 10,
+                             2, 5, 8, 11,
+                             3, 6, 9, 12);
+      const Matrix<T,2,4> m(m1*m2);
+      const Matrix<T,2,4> res(14, 32,  50,  68,
+                              32, 77, 122, 167);
+      CPPUNIT_ASSERT_EQUAL(res, m);
+    }
+    template<typename T>
+        void transpose() {
+      const Matrix<T,2,3> m(1, 2, 3,
+                            4, 5, 6);
+      const Matrix<T,3,2> res(1, 4,
+                              2, 5,
+                              3, 6);
+      CPPUNIT_ASSERT_EQUAL(res, m.t());
+    }
+    template<typename T>
+        void apply() {
+      Matrix<T,3,3> m(2, -2, 4,
+                      -2, 1, -6,
+                      1, 0, -2);
+      Matrix<T,3,3> o(m);
+      CPPUNIT_ASSERT_EQUAL((T)3*o, m.apply([](T& t){t*=3;}));
+    }
+    template<typename T>
+        void gauss() {
+      Matrix<T,3,3> m(2, -2, 4,
+                      -2, 1, -6,
+                      1, 0, -2);
+      const Matrix<T,3,3> res(1, -1, 2,
+                              0, -1, -2,
+                              0, 0, -6);
+      T lambda(m.gauss());
+      CPPUNIT_ASSERT_EQUAL(res, m);
+      CPPUNIT_ASSERT_EQUAL((T)2, lambda);
+    }
+    template<typename T>
+        void det() {
+      Matrix<T,3,3> m(2, -2, 4,
+                      -2, 1, -6,
+                      1, 0, -2);
+      CPPUNIT_ASSERT_EQUAL((T)12, m.det());
+    }
+    template<typename T>
+        void i() {
+      const Matrix<T,3,3> m1(1, 0, 0,
+                             0, 1, 0,
+                             0, 0, 1);
+      const Matrix<T,3,2> m2(1, 0,
+                             0, 1,
+                             0, 0);
+      const Matrix<T,2,3> m3(1, 0, 0,
+                             0, 1, 0);
+      CPPUNIT_ASSERT_EQUAL(m1, m1.i());
+      CPPUNIT_ASSERT_EQUAL(m2, m2.i());
+      CPPUNIT_ASSERT_EQUAL(m3, m3.i());
+    }
+    template<typename T>
+        void inv() {
+      {
+        Matrix<T,3,3> m(2, -1, 0,
+                        1, 2, -2,
+                        0, -1, 1);
+        const Matrix<T,3,3> res(0, 1, 2,
+                                -1, 2, 4,
+                                -1, 2, 5);
+        Matrix<T,3,3> o1(m), o2(m);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL((T)2*res, (T)2/o1);
+        CPPUNIT_ASSERT_EQUAL(o1.i(), o1/o2);
+        CPPUNIT_ASSERT_EQUAL(res, m);
+        
+      } {
+        Matrix<T,3,3> m(1, 2, 3,
+                        0, 1, 4,
+                        5, 6, 0);
+        const Matrix<T,3,3> res(-24, 18, 5,
+                                20, -15, -4,
+                                -5, 4, 1);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL(res, m);
+      } {
+        Matrix<T,3,3> m(1, 2, 3,
+                        0, 4, 5,
+                        1, 0, 6);
+        const Matrix<T,3,3> res((T)12/11, (T)-6/11, (T)-1/11,
+                                (T)5/22, (T)3/22, (T)-5/22,
+                                (T)-2/11, (T)1/11, (T)2/11);
+        Matrix<T,3,3> o1(m), o2(m);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL((T)2*res, (T)2/o1);
+        CPPUNIT_ASSERT_EQUAL(o1.i(), o1/o2);
+        CPPUNIT_ASSERT_EQUAL(res, m);
+      } {
+        Matrix<T,3,3> m(-2, 0, 1,
+                        9, 2, -3,
+                        5, 1, -2);
+        const Matrix<T,3,3> res(-1, 1, -2,
+                                3, -1, 3,
+                                -1, 2, -4);
+        Matrix<T,3,3> o(m);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL(m.i(), m*o);
+        CPPUNIT_ASSERT_EQUAL(res, m);
+      } {
+        Matrix<T,3,3> m(1, 3, 1,
+                        1, 1, 2,
+                        2, 3, 4);
+        const Matrix<T,3,3> res(2, 9, -5,
+                                0, -2, 1,
+                                -1, -3, 2);
+        Matrix<T,3,3> o(m);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL(m.i(), m*o);
+        CPPUNIT_ASSERT_EQUAL(res, m);
+      } {
+        Matrix<T,4,4> m(2, 1, 4, 1,
+                        -1, 1, 0, 2,
+                        0, 0, 2, 4,
+                        2, -2, 0, 1);
+        const Matrix<T,4,4> res((T)-1/3, (T)13/15, (T)-2/3, 0.6,
+                                (T)1/3, (T)16/15, (T)-2/3, 0.2,
+                                0, -0.8, 0.5, -0.4,
+                                0, 0.4, 0, 0.2);
+        Matrix<T,4,4> o(m);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL(m.i(), m*o);
+        CPPUNIT_ASSERT_EQUAL(res, m);
+      } {
+        Matrix<T,2,2> m(4, 3,
+                        3, 2);
+        const Matrix<T,2,2> res(-2, 3,
+                                3, -4);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL(res, m);
+      }
+    }
+    template<typename T>
+        void stream() {
+      const Matrix<T,3,4> m1(1, 2, 3, 4,
+                             5, 6, 7, 8,
+                             1, 4, 2, 8);
+      Matrix<T,3,4> m2;
+      std::string res("[3x4]{1,2,3,4,5,6,7,8,1,4,2,8}");
+      std::stringstream ss;
+      ss<<m1;
+      CPPUNIT_ASSERT_EQUAL(res, ss.str());
+      ss>>m2;
+      CPPUNIT_ASSERT_EQUAL(m1, m2);
+    }
+    CPPUNIT_TEST_SUITE(TemplateMatrixTest);
+    CPPUNIT_TEST(initFromArray1<int>);
+    CPPUNIT_TEST(initFromArray1<long>);
+    CPPUNIT_TEST(initFromArray1<unsigned>);
+    CPPUNIT_TEST(initFromArray1<unsigned long>);
+    CPPUNIT_TEST(initFromArray1<float>);
+    CPPUNIT_TEST(initFromArray1<double>);
+    CPPUNIT_TEST(initFromArray2<int>);
+    CPPUNIT_TEST(initFromArray2<long>);
+    CPPUNIT_TEST(initFromArray2<unsigned>);
+    CPPUNIT_TEST(initFromArray2<unsigned long>);
+    CPPUNIT_TEST(initFromArray2<float>);
+    CPPUNIT_TEST(initFromArray2<double>);
+    CPPUNIT_TEST(initFromArray3<int>);
+    CPPUNIT_TEST(initFromArray3<long>);
+    CPPUNIT_TEST(initFromArray3<unsigned>);
+    CPPUNIT_TEST(initFromArray3<unsigned long>);
+    CPPUNIT_TEST(initFromArray3<float>);
+    CPPUNIT_TEST(initFromArray3<double>);
+    CPPUNIT_TEST(initFromOtherMatrix<int>);
+    CPPUNIT_TEST(initFromOtherMatrix<long>);
+    CPPUNIT_TEST(initFromOtherMatrix<unsigned>);
+    CPPUNIT_TEST(initFromOtherMatrix<unsigned long>);
+    CPPUNIT_TEST(initFromOtherMatrix<float>);
+    CPPUNIT_TEST(initFromOtherMatrix<double>);
+    CPPUNIT_TEST(access<int>);
+    CPPUNIT_TEST(access<long>);
+    CPPUNIT_TEST(access<unsigned>);
+    CPPUNIT_TEST(access<unsigned long>);
+    CPPUNIT_TEST(access<float>);
+    CPPUNIT_TEST(access<double>);
+    CPPUNIT_TEST(equality<int>);
+    CPPUNIT_TEST(equality<long>);
+    CPPUNIT_TEST(equality<unsigned>);
+    CPPUNIT_TEST(equality<unsigned long>);
+    CPPUNIT_TEST(equality<float>);
+    CPPUNIT_TEST(equality<double>);
+    CPPUNIT_TEST(operator_plus<int>);
+    CPPUNIT_TEST(operator_plus<long>);
+    CPPUNIT_TEST(operator_plus<unsigned>);
+    CPPUNIT_TEST(operator_plus<unsigned long>);
+    CPPUNIT_TEST(operator_plus<float>);
+    CPPUNIT_TEST(operator_plus<double>);
+    CPPUNIT_TEST(operator_minus<int>);
+    CPPUNIT_TEST(operator_minus<long>);
+    CPPUNIT_TEST(operator_minus<unsigned>);
+    CPPUNIT_TEST(operator_minus<unsigned long>);
+    CPPUNIT_TEST(operator_minus<float>);
+    CPPUNIT_TEST(operator_minus<double>);
+    CPPUNIT_TEST(scalar_mult<int>);
+    CPPUNIT_TEST(scalar_mult<long>);
+    CPPUNIT_TEST(scalar_mult<unsigned>);
+    CPPUNIT_TEST(scalar_mult<unsigned long>);
+    CPPUNIT_TEST(scalar_mult<float>);
+    CPPUNIT_TEST(scalar_mult<double>);
+    CPPUNIT_TEST(scalar_div<int>);
+    CPPUNIT_TEST(scalar_div<long>);
+    CPPUNIT_TEST(scalar_div<unsigned>);
+    CPPUNIT_TEST(scalar_div<unsigned long>);
+    CPPUNIT_TEST(scalar_div<float>);
+    CPPUNIT_TEST(scalar_div<double>);
+    CPPUNIT_TEST(matrix_mult<int>);
+    CPPUNIT_TEST(matrix_mult<long>);
+    CPPUNIT_TEST(matrix_mult<unsigned>);
+    CPPUNIT_TEST(matrix_mult<unsigned long>);
+    CPPUNIT_TEST(matrix_mult<float>);
+    CPPUNIT_TEST(matrix_mult<double>);
+    CPPUNIT_TEST(transpose<int>);
+    CPPUNIT_TEST(transpose<long>);
+    CPPUNIT_TEST(transpose<unsigned>);
+    CPPUNIT_TEST(transpose<unsigned long>);
+    CPPUNIT_TEST(transpose<float>);
+    CPPUNIT_TEST(transpose<double>);
+    CPPUNIT_TEST(apply<int>);
+    CPPUNIT_TEST(apply<long>);
+    CPPUNIT_TEST(apply<unsigned>);
+    CPPUNIT_TEST(apply<unsigned long>);
+    CPPUNIT_TEST(apply<float>);
+    CPPUNIT_TEST(apply<double>);
+    CPPUNIT_TEST(gauss<int>);
+    CPPUNIT_TEST(gauss<long>);
+    CPPUNIT_TEST(gauss<float>);
+    CPPUNIT_TEST(gauss<double>);
+    CPPUNIT_TEST(det<int>);
+    CPPUNIT_TEST(det<long>);
+    CPPUNIT_TEST(det<float>);
+    CPPUNIT_TEST(det<double>);
+    CPPUNIT_TEST(i<int>);
+    CPPUNIT_TEST(i<long>);
+    CPPUNIT_TEST(i<unsigned>);
+    CPPUNIT_TEST(i<unsigned long>);
+    CPPUNIT_TEST(i<float>);
+    CPPUNIT_TEST(i<double>);
+    CPPUNIT_TEST(inv<float>);
+    CPPUNIT_TEST(inv<double>);
+    CPPUNIT_TEST(inv<long double>);
+    CPPUNIT_TEST(stream<int>);
+    CPPUNIT_TEST(stream<long>);
+    CPPUNIT_TEST(stream<unsigned>);
+    CPPUNIT_TEST(stream<unsigned long>);
+    CPPUNIT_TEST(stream<float>);
+    CPPUNIT_TEST(stream<double>);
+    CPPUNIT_TEST_SUITE_END();
+};
+CPPUNIT_TEST_SUITE_REGISTRATION(TemplateMatrixTest);
+
+class VariableMatrixTest: public CppUnit::TestFixture { 
+  public:
+    template<typename T>
+        void initFromArray1() {
+      Matrix<T> m(2,4,
+                  1, 2, 3, 4,
+                  5, 6, 7, 8);
+      Matrix<T> m2(2, 4);
+      m2 = {2, 4,
+            1, 2, 3, 4,
+            5, 6, 7, 8};
+      CPPUNIT_ASSERT_EQUAL(m, m2);
+      CPPUNIT_ASSERT_EQUAL((T)1, m[0][0]);
+      CPPUNIT_ASSERT_EQUAL((T)2, m[0][1]);
+      CPPUNIT_ASSERT_EQUAL((T)3, m[0][2]);
+      CPPUNIT_ASSERT_EQUAL((T)4, m[0][3]);
+      CPPUNIT_ASSERT_EQUAL((T)5, m[1][0]);
+      CPPUNIT_ASSERT_EQUAL((T)6, m[1][1]);
+      CPPUNIT_ASSERT_EQUAL((T)7, m[1][2]);
+      CPPUNIT_ASSERT_EQUAL((T)8, m[1][3]);
+    }
+    template<typename T>
+        void initFromArray2() {
+      Matrix<T> m(2, 4,
+                  1, 2, 3, 4,
+                  5, 6, 7, 8);
+      CPPUNIT_ASSERT_EQUAL((T)1, m[0][0]);
+      CPPUNIT_ASSERT_EQUAL((T)2, m[0][1]);
+      CPPUNIT_ASSERT_EQUAL((T)3, m[0][2]);
+      CPPUNIT_ASSERT_EQUAL((T)4, m[0][3]);
+      CPPUNIT_ASSERT_EQUAL((T)5, m[1][0]);
+      CPPUNIT_ASSERT_EQUAL((T)6, m[1][1]);
+      CPPUNIT_ASSERT_EQUAL((T)7, m[1][2]);
+      CPPUNIT_ASSERT_EQUAL((T)8, m[1][3]);
+    }
+    template<typename T>
+        void initFromOtherMatrix() {
+      Matrix<T> m1(2, 4,
+                   1, 2, 3, 4,
+                   5, 6, 7, 8);
+      m1[1][2] = 13;
+      Matrix<T> m2(m1);
+      m1[0][2] = 16;
+      m2[0][0] = 0;
+      Matrix<T> m3(2, 4);
+      m3 = m2;
+      CPPUNIT_ASSERT_EQUAL(m2, m3);
+      CPPUNIT_ASSERT_EQUAL((T)1, m1[0][0]);
+      CPPUNIT_ASSERT_EQUAL((T)2, m1[0][1]);
+      CPPUNIT_ASSERT_EQUAL((T)16, m1[0][2]);
+      CPPUNIT_ASSERT_EQUAL((T)4, m1[0][3]);
+      CPPUNIT_ASSERT_EQUAL((T)5, m1[1][0]);
+      CPPUNIT_ASSERT_EQUAL((T)6, m1[1][1]);
+      CPPUNIT_ASSERT_EQUAL((T)13, m1[1][2]);
+      CPPUNIT_ASSERT_EQUAL((T)8, m1[1][3]);
+      CPPUNIT_ASSERT_EQUAL((T)0, m2[0][0]);
+      CPPUNIT_ASSERT_EQUAL((T)2, m2[0][1]);
+      CPPUNIT_ASSERT_EQUAL((T)3, m2[0][2]);
+      CPPUNIT_ASSERT_EQUAL((T)4, m2[0][3]);
+      CPPUNIT_ASSERT_EQUAL((T)5, m2[1][0]);
+      CPPUNIT_ASSERT_EQUAL((T)6, m2[1][1]);
+      CPPUNIT_ASSERT_EQUAL((T)13, m2[1][2]);
+      CPPUNIT_ASSERT_EQUAL((T)8, m2[1][3]);
+    }
+    template<typename T>
+        void initFromDefault() {
+      Matrix<T> m(2, 4);
+      CPPUNIT_ASSERT_EQUAL((T)0, m[0][0]);
+      CPPUNIT_ASSERT_EQUAL((T)0, m[0][1]);
+      CPPUNIT_ASSERT_EQUAL((T)0, m[1][0]);
+      CPPUNIT_ASSERT_EQUAL((T)0, m[1][1]);
+    }
+    template<typename T>
+        void access() {
+      Matrix<T> m1(2, 4,
+                   1, 2, 3, 4,
+                   5, 6, 7, 8);
+      for (size_t row(0); row<2; ++row)
+        for (size_t column(0); column<4; ++column)
+          CPPUNIT_ASSERT_EQUAL(m1(row, column), m1[row][column]);
+    }
+    template<typename T>
+        void equality() {
+      Matrix<T> m1(2, 4,
+                   1, 2, 3, 4,
+                   5, 6, 7, 8);
+      Matrix<T> m2(2, 4,
+                   1, 2, 3, 4,
+                   5, 6, 7, 8);
+      Matrix<T> m3(2, 4,
+                   1, 2, 3, 4,
+                   5, 6, 7, 9);
+      Matrix<T> m4(2, 4,
+                   9, 2, 3, 4,
+                   5, 6, 7, 8);
+      CPPUNIT_ASSERT(m1==m2);
+      CPPUNIT_ASSERT(m1!=m3);
+      CPPUNIT_ASSERT(m1!=m4);
+    }
+    template<typename T>
+        void operator_plus() {
+      const Matrix<T> m1(2, 4,
+                         1, 2, 3, 4,
+                         5, 6, 7, 8);
+      const Matrix<T> m2(2, 4,
+                         2, 4, 6, 8,
+                         1, 3, 5, 7);
+      const Matrix<T> m(m1+m2);
+      const Matrix<T> res(2, 4,
+                          3, 6, 9, 12,
+                          6, 9, 12, 15);
+      CPPUNIT_ASSERT_EQUAL(res, m);
+    }
+    template<typename T>
+        void operator_minus() {
+      const Matrix<T> m1(2, 4,
+                         1, 2, 3, 4,
+                         5, 6, 7, 8);
+      const Matrix<T> m2(2, 4,
+                         2, 4, 6, 8,
+                         1, 3, 5, 7);
+      const Matrix<T> m(m1-m2);
+      const Matrix<T> res(2, 4,
+                         -1, -2, -3, -4,
+                          4, 3, 2, 1);
+      CPPUNIT_ASSERT_EQUAL(res, m);
+      CPPUNIT_ASSERT_EQUAL(res, -m2+m1);
+    }
+    template<typename T>
+        void scalar_mult() {
+      const Matrix<T> m1(2, 4,
+                         1, 2, 3, 4,
+                         5, 6, 7, 8);
+      T two(2), three(3), four(4), big(32);
+      CPPUNIT_ASSERT_EQUAL(m1*two, m1+m1);
+      CPPUNIT_ASSERT_EQUAL(m1*three, m1+m1+m1);
+      CPPUNIT_ASSERT_EQUAL(m1*four, m1+m1+m1+m1);
+      CPPUNIT_ASSERT_EQUAL(two*m1, m1*two);
+      CPPUNIT_ASSERT_EQUAL(three*m1, m1*three);
+      CPPUNIT_ASSERT_EQUAL(four*m1, m1*four);
+      CPPUNIT_ASSERT_EQUAL(two*m1*two, m1*four);
+      CPPUNIT_ASSERT_EQUAL(big*m1*four, m1*four*big);
+    }
+    template<typename T>
+        void scalar_div() {
+      const Matrix<T> m1(2, 4,
+                         2, 4, 6, 8,
+                         10, 12, 14, 16);
+      const Matrix<T> m2(2, 4,
+                         1, 2, 3, 4,
+                         5, 6, 7, 8);
+      CPPUNIT_ASSERT_EQUAL(m2, m1/(T)2);
+      CPPUNIT_ASSERT_EQUAL(m2, (T)3*m1/(T)6);
+    }
+    template<typename T>
+        void matrix_mult() {
+      const Matrix<T> m1(2, 3,
+                         1, 2, 3,
+                         4, 5, 6);
+      const Matrix<T> m2(3, 4,
+                         1, 4, 7, 10,
+                         2, 5, 8, 11,
+                         3, 6, 9, 12);
+      const Matrix<T> m(m1*m2);
+      const Matrix<T> res(2, 4,
+                          14, 32,  50,  68,
+                          32, 77, 122, 167);
+      CPPUNIT_ASSERT_EQUAL(res, m);
+    }
+    template<typename T>
+        void transpose() {
+      const Matrix<T> m(2, 3,
+                        1, 2, 3,
+                        4, 5, 6);
+      const Matrix<T> res(3, 2,
+                          1, 4,
+                          2, 5,
+                          3, 6);
+      CPPUNIT_ASSERT_EQUAL(res, m.t());
+    }
+    template<typename T>
+        void apply() {
+          Matrix<T> m(3, 3,
+                      2, -2, 4,
+                      -2, 1, -6,
+                      1, 0, -2);
+          Matrix<T> o(m);
+          CPPUNIT_ASSERT_EQUAL((T)3*o, m.apply([](T& t){t*=3;}));
+    }
+    template<typename T>
+        void gauss() {
+      Matrix<T> m(3, 3,
+                  2, -2, 4,
+                  -2, 1, -6,
+                  1, 0, -2);
+      const Matrix<T> res(3,3,
+                          1, -1, 2,
+                          0, -1, -2,
+                          0, 0, -6);
+      T lambda(m.gauss());
+      CPPUNIT_ASSERT_EQUAL(res, m);
+      CPPUNIT_ASSERT_EQUAL((T)2, lambda);
+    }
+    template<typename T>
+        void det() {
+      Matrix<T> m(3, 3,
+                  2, -2, 4,
+                  -2, 1, -6,
+                  1, 0, -2);
+      CPPUNIT_ASSERT_EQUAL((T)12, m.det());
+    }
+    template<typename T>
+        void i() {
+      const Matrix<T> m1(3, 3,
+                         1, 0, 0,
+                         0, 1, 0,
+                         0, 0, 1);
+      const Matrix<T> m2(3, 2,
+                         1, 0,
+                         0, 1,
+                         0, 0);
+      const Matrix<T> m3(2, 3,
+                         1, 0, 0,
+                         0, 1, 0);
+      CPPUNIT_ASSERT_EQUAL(m1, m1.i());
+      CPPUNIT_ASSERT_EQUAL(m2, m2.i());
+      CPPUNIT_ASSERT_EQUAL(m3, m3.i());
+    }
+    template<typename T>
+        void inv() {
+      {
+        Matrix<T> m(3, 3,
+                    2, -1, 0,
+                    1, 2, -2,
+                    0, -1, 1);
+        const Matrix<T> res(3, 3,
+                            0, 1, 2,
+                            -1, 2, 4,
+                            -1, 2, 5);
+        Matrix<T> o1(m), o2(m);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL((T)2*res, (T)2/o1);
+        CPPUNIT_ASSERT_EQUAL(o1.i(), o1/o2);
+        CPPUNIT_ASSERT_EQUAL(res, m);
+      } {
+        Matrix<T> m(3, 3,
+                    1, 2, 3,
+                    0, 1, 4,
+                    5, 6, 0);
+        const Matrix<T> res(3, 3,
+                            -24, 18, 5,
+                            20, -15, -4,
+                            -5, 4, 1);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL(res, m);
+      } {
+        Matrix<T> m(3, 3,
+                    1, 2, 3,
+                    0, 4, 5,
+                    1, 0, 6);
+        const Matrix<T> res(3, 3,
+                            (T)12/11, (T)-6/11, (T)-1/11,
+                            (T)5/22, (T)3/22, (T)-5/22,
+                            (T)-2/11, (T)1/11, (T)2/11);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL(res, m);
+      } {
+        Matrix<T> m(3, 3,
+                    -2, 0, 1,
+                    9, 2, -3,
+                    5, 1, -2);
+        const Matrix<T> res(3, 3,
+                            -1, 1, -2,
+                            3, -1, 3,
+                            -1, 2, -4);
+        Matrix<T> o(m);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL(m.i(), m*o);
+        CPPUNIT_ASSERT_EQUAL(res, m);
+      } {
+        Matrix<T> m(4, 4,
+                    2, 1, 4, 1,
+                    -1, 1, 0, 2,
+                    0, 0, 2, 4,
+                    2, -2, 0, 1);
+        const Matrix<T> res(4, 4,
+                            (T)-1/3, (T)13/15, (T)-2/3, 0.6,
+                            (T)1/3, (T)16/15, (T)-2/3, 0.2,
+                            0, -0.8, 0.5, -0.4,
+                            0, 0.4, 0, 0.2);
+        Matrix<T> o(m);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL(m.i(), m*o);
+        CPPUNIT_ASSERT_EQUAL(res, m);
+      } {
+        Matrix<T> m(2, 2,
+                    4, 3,
+                    3, 2);
+        const Matrix<T> res(2, 2,
+                            -2, 3,
+                            3, -4);
+        m.inv();
+        CPPUNIT_ASSERT_EQUAL(res, m);
+      }
+    }
+    template<typename T>
+        void stream() {
+      const Matrix<T> m1(3, 4,
+                         1, 2, 3, 4,
+                         5, 6, 7, 8,
+                         1, 4, 2, 8);
+      Matrix<T> m2(1, 1);
+      std::string res("[3x4]{1,2,3,4,5,6,7,8,1,4,2,8}");
+      std::stringstream ss;
+      ss<<m1;
+      CPPUNIT_ASSERT_EQUAL(res, ss.str());
+      ss>>m2;
+      CPPUNIT_ASSERT_EQUAL(m1, m2);
     }
-    CPPUNIT_TEST_SUITE(DummyTest);
-    CPPUNIT_TEST(dummy);
+    CPPUNIT_TEST_SUITE(VariableMatrixTest);
+    CPPUNIT_TEST(initFromArray1<int>);
+    CPPUNIT_TEST(initFromArray1<long>);
+    CPPUNIT_TEST(initFromArray1<unsigned>);
+    CPPUNIT_TEST(initFromArray1<unsigned long>);
+    CPPUNIT_TEST(initFromArray1<float>);
+    CPPUNIT_TEST(initFromArray1<double>);
+    CPPUNIT_TEST(initFromArray2<int>);
+    CPPUNIT_TEST(initFromArray2<long>);
+    CPPUNIT_TEST(initFromArray2<unsigned>);
+    CPPUNIT_TEST(initFromArray2<unsigned long>);
+    CPPUNIT_TEST(initFromArray2<float>);
+    CPPUNIT_TEST(initFromArray2<double>);
+    CPPUNIT_TEST(initFromOtherMatrix<int>);
+    CPPUNIT_TEST(initFromOtherMatrix<long>);
+    CPPUNIT_TEST(initFromOtherMatrix<unsigned>);
+    CPPUNIT_TEST(initFromOtherMatrix<unsigned long>);
+    CPPUNIT_TEST(initFromOtherMatrix<float>);
+    CPPUNIT_TEST(initFromOtherMatrix<double>);
+    CPPUNIT_TEST(access<int>);
+    CPPUNIT_TEST(access<long>);
+    CPPUNIT_TEST(access<unsigned>);
+    CPPUNIT_TEST(access<unsigned long>);
+    CPPUNIT_TEST(access<float>);
+    CPPUNIT_TEST(access<double>);
+    CPPUNIT_TEST(equality<int>);
+    CPPUNIT_TEST(equality<long>);
+    CPPUNIT_TEST(equality<unsigned>);
+    CPPUNIT_TEST(equality<unsigned long>);
+    CPPUNIT_TEST(equality<float>);
+    CPPUNIT_TEST(equality<double>);
+    CPPUNIT_TEST(operator_plus<int>);
+    CPPUNIT_TEST(operator_plus<long>);
+    CPPUNIT_TEST(operator_plus<unsigned>);
+    CPPUNIT_TEST(operator_plus<unsigned long>);
+    CPPUNIT_TEST(operator_plus<float>);
+    CPPUNIT_TEST(operator_plus<double>);
+    CPPUNIT_TEST(operator_minus<int>);
+    CPPUNIT_TEST(operator_minus<long>);
+    CPPUNIT_TEST(operator_minus<unsigned>);
+    CPPUNIT_TEST(operator_minus<unsigned long>);
+    CPPUNIT_TEST(operator_minus<float>);
+    CPPUNIT_TEST(operator_minus<double>);
+    CPPUNIT_TEST(scalar_mult<int>);
+    CPPUNIT_TEST(scalar_mult<long>);
+    CPPUNIT_TEST(scalar_mult<unsigned>);
+    CPPUNIT_TEST(scalar_mult<unsigned long>);
+    CPPUNIT_TEST(scalar_mult<float>);
+    CPPUNIT_TEST(scalar_mult<double>);
+    CPPUNIT_TEST(scalar_mult<int>);
+    CPPUNIT_TEST(scalar_div<int>);
+    CPPUNIT_TEST(scalar_div<long>);
+    CPPUNIT_TEST(scalar_div<unsigned>);
+    CPPUNIT_TEST(scalar_div<unsigned long>);
+    CPPUNIT_TEST(scalar_div<float>);
+    CPPUNIT_TEST(scalar_div<double>);
+    CPPUNIT_TEST(matrix_mult<int>);
+    CPPUNIT_TEST(matrix_mult<long>);
+    CPPUNIT_TEST(matrix_mult<unsigned>);
+    CPPUNIT_TEST(matrix_mult<unsigned long>);
+    CPPUNIT_TEST(matrix_mult<float>);
+    CPPUNIT_TEST(matrix_mult<double>);
+    CPPUNIT_TEST(transpose<int>);
+    CPPUNIT_TEST(transpose<long>);
+    CPPUNIT_TEST(transpose<unsigned>);
+    CPPUNIT_TEST(transpose<unsigned long>);
+    CPPUNIT_TEST(transpose<float>);
+    CPPUNIT_TEST(transpose<double>);
+    CPPUNIT_TEST(apply<int>);
+    CPPUNIT_TEST(apply<long>);
+    CPPUNIT_TEST(apply<unsigned>);
+    CPPUNIT_TEST(apply<unsigned long>);
+    CPPUNIT_TEST(apply<float>);
+    CPPUNIT_TEST(apply<double>);
+    CPPUNIT_TEST(gauss<int>);
+    CPPUNIT_TEST(gauss<long>);
+    CPPUNIT_TEST(gauss<float>);
+    CPPUNIT_TEST(gauss<double>);
+    CPPUNIT_TEST(det<int>);
+    CPPUNIT_TEST(det<long>);
+    CPPUNIT_TEST(det<float>);
+    CPPUNIT_TEST(det<double>);
+    CPPUNIT_TEST(i<int>);
+    CPPUNIT_TEST(i<long>);
+    CPPUNIT_TEST(i<unsigned>);
+    CPPUNIT_TEST(i<unsigned long>);
+    CPPUNIT_TEST(i<float>);
+    CPPUNIT_TEST(i<double>);
+    CPPUNIT_TEST(inv<float>);
+    CPPUNIT_TEST(inv<double>);
+    CPPUNIT_TEST(inv<long double>);
+    CPPUNIT_TEST(stream<int>);
+    CPPUNIT_TEST(stream<long>);
+    CPPUNIT_TEST(stream<unsigned>);
+    CPPUNIT_TEST(stream<unsigned long>);
+    CPPUNIT_TEST(stream<float>);
+    CPPUNIT_TEST(stream<double>);
     CPPUNIT_TEST_SUITE_END();
 };
-CPPUNIT_TEST_SUITE_REGISTRATION(DummyTest);
+CPPUNIT_TEST_SUITE_REGISTRATION(VariableMatrixTest);
 
 int main(int argc, char** argv) try {
   std::ofstream ofs((*argv+std::string(".xml")).c_str());