works perfectly, fully tested

configure.ac

@@ -9,8 +9,9 @@
 ## 45678901234567890123456789012345678901234567890123456789012345678901234567890
 
 m4_define(x_package_name, libmatricxx) # project's name
-m4_define(x_major, 0) # project's major version
+m4_define(x_major, 1) # project's major version
-m4_define(x_minor, 2) # project's minor version
+m4_define(x_minor, 0) # project's minor version
+m4_define(x_least_diff, 12) # 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])

src/matrix.hxx

@@ -6,27 +6,26 @@
 // 45678901234567890123456789012345678901234567890123456789012345678901234567890
 
 #include <iostream>
+#include <sstream>
 #include <cstring>
 #include <cassert>
 #include <type_traits>
 #include <limits>
+#include <cmath>
+#include <cfenv>
+#include <stdexcept>
+#include  <functional>
 
 namespace math {
 
   template<typename TYPE>
       bool almostEqual(TYPE a, TYPE b) {
-    if (a<0)
+    a = std::fabs(a);
-      if (b<0) {
+    b = std::fabs(b);
-        a = -a;
+    TYPE diff(std::fabs(a-b));
-        b = -b;
-      } else return a==b; // different signs
-    else if (b<0) return a==b; // different signs
-    TYPE diff(abs(a-b));
-    a = abs(a);
-    b = abs(b);
     TYPE max(a>b?a:b);
-    if (max<1) return diff<=std::numeric_limits<TYPE>::epsilon();
+    if (max<1) return diff<=1000*std::numeric_limits<TYPE>::epsilon();
-    return diff<=max*std::numeric_limits<TYPE>::epsilon();
+    return diff<=max*1000*std::numeric_limits<TYPE>::epsilon();
   }
 
   template<typename TYPE>
@@ -34,6 +33,11 @@ namespace math {
     return a==b;
   }
 
+  template<>
+      bool equal(const long double& a, const long double& b) {
+    return almostEqual(a, b);
+  }
+
   template<>
       bool equal(const double& a, const double& b) {
     return almostEqual(a, b);
@@ -158,12 +162,6 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
     /// @name operators
     ///{
 
-    MatrixBase& operator=(TYPE oc[]) {
-      assert(sizeof(oc)==MEM_SIZE);
-      memcpy(_c, oc, MEM_SIZE);
-      return *this;
-    }
-
     MatrixBase& operator=(const MatrixBase& o) {
       assert_check(o);
       memcpy(_c, o._c, MEM_SIZE);
@@ -212,11 +210,10 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
     /// @name special operations
     ///{
 
-    bool similar(const MatrixBase& o, const TYPE& diff) const {
+    MatrixBase& apply(std::function<void(TYPE&)> fn) {
-      if (!check(o)) return false;
+      TYPE *to((TYPE*)(_c)+SIZE);
-      TYPE *to((TYPE*)(_c)+SIZE), *from((TYPE*)(o._c)+SIZE);
+      while (to>(TYPE*)(_c)) fn(*--to);
-      while (to>(TYPE*)(_c)) if (abs(*--to - *--from)>diff) return false;
+      return *this;
-      return true;
     }
 
     TYPE det() {
@@ -230,6 +227,10 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
       /// @see http://www.mathebibel.de/determinante-berechnen-nach-gauss
       /// 1. normalize first line to first value
       TYPE lambda(at(0, 0));
+      if (lambda==0) {
+        feraiseexcept(FE_DIVBYZERO);
+        throw std::range_error("gauss calculation failed");
+      }
       at(0, 0) = 1;
       for (TYPE *p((TYPE*)(_c)+COLUMNS); p>(TYPE*)(_c)+1;) *--p/=lambda;
       /// 2. nullify lower triangle
@@ -297,15 +298,34 @@ template<typename TYPE, size_t TROWS=0, size_t TCOLUMNS=0> class Matrix:
     /// @name operators
     ///{
 
-    Matrix& operator=(TYPE oc[TROWS][TCOLUMNS]) {
+    Matrix& operator=(const Matrix& o) {
-      memcpy(Parent::_c, oc, Parent::MEM_SIZE);
+      Parent::operator=(o);
+      return *this;
+    }
+
+    Matrix& operator+=(const Matrix& o) {
+      Parent::operator+=(o);
+      return *this;
+    }
+
+    Matrix& operator-=(const Matrix& o) {
+      Parent::operator-=(o);
+      return *this;
+    }
+
+    Matrix& operator*=(const TYPE& o) {
+      Parent::operator*=(o);
+      return *this;
+    }
+
+    Matrix& operator/=(const TYPE& o) {
+      Parent::operator/=(o);
       return *this;
     }
 
     Matrix operator-() const {
-      Matrix res;
+      Matrix res(*this);
-      for (TYPE *to((TYPE*)(res._c)+Parent::SIZE); to>(TYPE*)(Parent::_c); --to)
+      for (TYPE *to((TYPE*)(res._c)+this->SIZE); to>(TYPE*)(res._c); *--to = -*to);
-        *to = -*to;
       return res;
     }
 
@@ -325,6 +345,11 @@ template<typename TYPE, size_t TROWS=0, size_t TCOLUMNS=0> class Matrix:
     /// @name special operations
     ///{
 
+    Matrix& apply(std::function<void(TYPE&)> fn) {
+      Parent::apply(fn);
+      return *this;
+    }
+
     Matrix<TYPE, TCOLUMNS, TROWS> t() const {
       Matrix<TYPE, TCOLUMNS, TROWS> res;
       for (size_t row(0); row<TROWS; ++row)
@@ -440,10 +465,34 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
     /// @name operators
     ///{
 
+    Matrix& operator=(const Matrix& o) {
+      Parent::operator=(o);
+      return *this;
+    }
+
+    Matrix& operator+=(const Matrix& o) {
+      Parent::operator+=(o);
+      return *this;
+    }
+
+    Matrix& operator-=(const Matrix& o) {
+      Parent::operator-=(o);
+      return *this;
+    }
+
+    Matrix& operator*=(const TYPE& o) {
+      Parent::operator*=(o);
+      return *this;
+    }
+
+    Matrix& operator/=(const TYPE& o) {
+      Parent::operator/=(o);
+      return *this;
+    }
+
     Matrix operator-() const {
-      Matrix res(Parent::COLUMNS, Parent::ROWS);
+      Matrix res(*this);
-      for (TYPE *to((TYPE*)(res._c)+Parent::SIZE); to>(TYPE*)(Parent::_c); --to)
+      for (TYPE *to((TYPE*)(res._c)+this->SIZE); to>(TYPE*)(res._c); *--to = -*to);
-        *to = -*to;
       return res;
     }
 
@@ -462,6 +511,11 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
     ///@name special operations
     ///{
 
+    Matrix& apply(std::function<void(TYPE&)> fn) {
+      Parent::apply(fn);
+      return *this;
+    }
+
     Matrix t() const {
       Matrix res(this->COLUMNS, this->ROWS);
       for (size_t row(0); row<this->ROWS; ++row)
@@ -582,12 +636,61 @@ template<typename TYPE, size_t ROWS, size_t COLUMNS>
    return res;
 }
 
+template<typename TYPE, size_t ROWS, size_t COLUMNS>
+    Matrix<TYPE, ROWS, COLUMNS> operator/(const Matrix<TYPE, ROWS, COLUMNS>& m,
+                                          const TYPE& v) {
+  Matrix<TYPE, ROWS, COLUMNS> res(m);
+  res /= v;
+  return res;
+}
+
+template<typename TYPE, size_t ROWS, size_t COLUMNS>
+    Matrix<TYPE, ROWS, COLUMNS> operator/(const TYPE& v,
+                                          const Matrix<TYPE, ROWS, COLUMNS>& m) {
+  Matrix<TYPE, ROWS, COLUMNS> res(m);
+  res.inv() *= v;
+  return res;
+}
+
+template<typename TYPE, size_t ROWS, size_t COLUMNS>
+    Matrix<TYPE, ROWS, COLUMNS> operator/(const Matrix<TYPE, ROWS, COLUMNS>& m1,
+                                          const Matrix<TYPE, ROWS, COLUMNS>& m2) {
+  Matrix<TYPE, ROWS, COLUMNS> res(m2);
+  return m1 * res.inv();
+}
+
 template<typename TYPE, size_t ROWS, size_t COLUMNS>
     std::ostream& operator<<(std::ostream& s, const Matrix<TYPE, ROWS, COLUMNS>& m) {
-  for (size_t w = 0; w < m.ROWS; ++w) {
+  s<<'['<<ROWS<<'x'<<COLUMNS<<"]{";
-    for (size_t h = 0; h < m.COLUMNS;++h) {
+  for (size_t row = 0; row < m.ROWS; ++row) {
-      s<<m[w][h]<<' ';
+    for (size_t column = 0; column < m.COLUMNS; ++column) {
+      if (row!=0||column!=0) s<<',';
+      s<<m(row, column);
+    }
+  }
+  return s<<'}';
+}
+
+template<typename TYPE, size_t ROWS, size_t COLUMNS>
+    std::istream& operator>>(std::istream& in, Matrix<TYPE, ROWS, COLUMNS>& m) {
+  std::ios_base::failure err("illegal matrix format");
+  char c(0);
+  size_t sz(0);
+  TYPE val(0);
+  std::string s;
+  if (!in.get(c) || c!='[') throw err;
+  if (!std::getline(in, s, 'x') || !(std::stringstream(s)>>sz) || sz!=m.ROWS) throw err;
+  if (!std::getline(in, s, ']') || !(std::stringstream(s)>>sz) || sz!=m.COLUMNS) throw err;
+  if (!in.get(c) || c!='{') throw err;
+  for (size_t row = 0; row < m.ROWS; ++row) {
+    for (size_t column = 0; column < m.COLUMNS; ++column) {
+      if (row==m.ROWS-1&&column==m.COLUMNS) {
+        if (!std::getline(in, s, '}') || !(std::stringstream(s)>>val)) throw err;
+      } else {
+        if (!std::getline(in, s, ',') || !(std::stringstream(s)>>val)) throw err;
+      }
+      m(row, column) = val;
     }
-    s<<'\n';
   }
+  return in;
 }

test/basic.cxx

@@ -19,6 +19,10 @@ class TemplateMatrixTest: public CppUnit::TestFixture {
         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]);
@@ -62,6 +66,9 @@ class TemplateMatrixTest: public CppUnit::TestFixture {
       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]);
@@ -133,6 +140,7 @@ class TemplateMatrixTest: public CppUnit::TestFixture {
       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() {
@@ -148,6 +156,15 @@ class TemplateMatrixTest: public CppUnit::TestFixture {
       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,
@@ -169,6 +186,14 @@ class TemplateMatrixTest: public CppUnit::TestFixture {
                               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,
@@ -211,8 +236,12 @@ class TemplateMatrixTest: public CppUnit::TestFixture {
         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(m.similar(res, 0.0001));
+        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,
@@ -229,7 +258,34 @@ class TemplateMatrixTest: public CppUnit::TestFixture {
         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,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,
@@ -240,6 +296,19 @@ class TemplateMatrixTest: public CppUnit::TestFixture {
         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>);
@@ -295,6 +364,12 @@ class TemplateMatrixTest: public CppUnit::TestFixture {
     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>);
@@ -307,6 +382,12 @@ class TemplateMatrixTest: public CppUnit::TestFixture {
     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>);
@@ -323,6 +404,13 @@ class TemplateMatrixTest: public CppUnit::TestFixture {
     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);
@@ -334,6 +422,11 @@ class VariableMatrixTest: public CppUnit::TestFixture {
       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]);
@@ -345,9 +438,9 @@ class VariableMatrixTest: public CppUnit::TestFixture {
     }
     template<typename T>
         void initFromArray2() {
-        Matrix<T> m(2, 4,
+      Matrix<T> m(2, 4,
-                    1, 2, 3, 4,
+                  1, 2, 3, 4,
-                    5, 6, 7, 8);
+                  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]);
@@ -366,6 +459,9 @@ class VariableMatrixTest: public CppUnit::TestFixture {
       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]);
@@ -445,6 +541,7 @@ class VariableMatrixTest: public CppUnit::TestFixture {
                          -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() {
@@ -461,6 +558,17 @@ class VariableMatrixTest: public CppUnit::TestFixture {
       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,
@@ -487,6 +595,15 @@ class VariableMatrixTest: public CppUnit::TestFixture {
                           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,
@@ -537,8 +654,11 @@ class VariableMatrixTest: public CppUnit::TestFixture {
                             0, 1, 2,
                             -1, 2, 4,
                             -1, 2, 5);
+        Matrix<T> o1(m), o2(m);
         m.inv();
-        CPPUNIT_ASSERT(m.similar(res, 0.0001));
+        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,
@@ -550,6 +670,54 @@ class VariableMatrixTest: public CppUnit::TestFixture {
                             -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);
       }
     }
     CPPUNIT_TEST_SUITE(VariableMatrixTest);
@@ -602,11 +770,12 @@ class VariableMatrixTest: public CppUnit::TestFixture {
     CPPUNIT_TEST(scalar_mult<float>);
     CPPUNIT_TEST(scalar_mult<double>);
     CPPUNIT_TEST(scalar_mult<int>);
-    CPPUNIT_TEST(scalar_mult<long>);
+    CPPUNIT_TEST(scalar_div<int>);
-    CPPUNIT_TEST(scalar_mult<unsigned>);
+    CPPUNIT_TEST(scalar_div<long>);
-    CPPUNIT_TEST(scalar_mult<unsigned long>);
+    CPPUNIT_TEST(scalar_div<unsigned>);
-    CPPUNIT_TEST(scalar_mult<float>);
+    CPPUNIT_TEST(scalar_div<unsigned long>);
-    CPPUNIT_TEST(scalar_mult<double>);
+    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>);
@@ -619,6 +788,12 @@ class VariableMatrixTest: public CppUnit::TestFixture {
     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>);
@@ -635,6 +810,7 @@ class VariableMatrixTest: public CppUnit::TestFixture {
     CPPUNIT_TEST(i<double>);
     CPPUNIT_TEST(inv<float>);
     CPPUNIT_TEST(inv<double>);
+    CPPUNIT_TEST(inv<long double>);
     CPPUNIT_TEST_SUITE_END();
 };
 CPPUNIT_TEST_SUITE_REGISTRATION(VariableMatrixTest);