/*! @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);
      Matrix<T,2,4> m6(-m1);
      Matrix<T,2,4> m7(-1, -2, -3, -4,
                       -5, -6, -7, -8);
      CPPUNIT_ASSERT_EQUAL(m7, m6);
      CPPUNIT_ASSERT(m1==m2);
      CPPUNIT_ASSERT(m1!=m3);
      CPPUNIT_ASSERT(m1!=m4);
      CPPUNIT_ASSERT(m1!=m5);
      CPPUNIT_ASSERT(m1!=m6);
    }
    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 norm() {
      Matrix<T,3,3> m1(1, -1, 1,
                       -1, 1, -1,
                       1, 1, -1);
      Matrix<T,3,3> m2(-1, 3, 2,
                       -2, 0, -3,
                       2, 1, -2);
      Matrix<T,3,3> m3(0, 3, 0,
                       -1, 0, 2,
                       0, -3, 1);
      CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0, m1.norm(), 0.00001);
      CPPUNIT_ASSERT_DOUBLES_EQUAL(6.0, m2.norm(), 0.00001);
      CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0, m3.norm(3), 0.00001);
    }
    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(initFromArray1<long 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(initFromArray2<long 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(initFromArray3<long 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(initFromOtherMatrix<long 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(access<long 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(equality<long 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_plus<long 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(operator_minus<long 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<long 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(scalar_div<long 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(matrix_mult<long 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(transpose<long 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(apply<long double>);
    CPPUNIT_TEST(norm<int>);
    CPPUNIT_TEST(norm<long>);
    CPPUNIT_TEST(norm<float>);
    CPPUNIT_TEST(norm<double>);
    CPPUNIT_TEST(norm<long double>);
    CPPUNIT_TEST(gauss<int>);
    CPPUNIT_TEST(gauss<long>);
    CPPUNIT_TEST(gauss<float>);
    CPPUNIT_TEST(gauss<double>);
    CPPUNIT_TEST(gauss<long double>);
    CPPUNIT_TEST(det<int>);
    CPPUNIT_TEST(det<long>);
    CPPUNIT_TEST(det<float>);
    CPPUNIT_TEST(det<double>);
    CPPUNIT_TEST(det<long 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(i<long 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(stream<long 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);
      Matrix<T> m5(-m1);
      Matrix<T> m6(2, 4,
                   -1, -2, -3, -4,
                   -5, -6, -7, -8);
      CPPUNIT_ASSERT_EQUAL(m6, m5);
      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> 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 norm() {
      Matrix<T> m1(3, 3,
                   1, -1, 1,
                   -1, 1, -1,
                   1, 1, -1);
      Matrix<T> m2(3, 3,
                   -1, 3, 2,
                   -2, 0, -3,
                   2, 1, -2);
      Matrix<T> m3(3, 3,
                   0, 3, 0,
                   -1, 0, 2,
                   0, -3, 1);
      CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0, m1.norm(), 0.00001);
      CPPUNIT_ASSERT_DOUBLES_EQUAL(6.0, m2.norm(), 0.00001);
      CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0, m3.norm(3), 0.00001);
    }
    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(initFromArray1<long 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(initFromArray2<long 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(initFromOtherMatrix<long 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(access<long 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(equality<long 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_plus<long 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(operator_minus<long 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<long 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(scalar_div<long 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(matrix_mult<long 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(transpose<long 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(apply<long double>);
    CPPUNIT_TEST(norm<int>);
    CPPUNIT_TEST(norm<long>);
    CPPUNIT_TEST(norm<float>);
    CPPUNIT_TEST(norm<double>);
    CPPUNIT_TEST(norm<long double>);
    CPPUNIT_TEST(gauss<int>);
    CPPUNIT_TEST(gauss<long>);
    CPPUNIT_TEST(gauss<float>);
    CPPUNIT_TEST(gauss<double>);
    CPPUNIT_TEST(gauss<long double>);
    CPPUNIT_TEST(det<int>);
    CPPUNIT_TEST(det<long>);
    CPPUNIT_TEST(det<float>);
    CPPUNIT_TEST(det<double>);
    CPPUNIT_TEST(det<long 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(i<long 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(stream<long 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;
 }