Math matrix, mainly for graphics

Question

I wrote a matrix library to replace my old matrix library. I was using my old matrix library to do basic operations for 3D rendering but that library grew old and now I'm seeking to replace it. I'm mostly looking for comfort while writing code and I want to know what you think about the coding style and the simplicity in usage. I'm also expecting it to perform, mainly low memory usage and speed but that isn't my biggest concern.

Matrix.h:

#ifndef MATRIX_H_INCLUDED
#define MATRIX_H_INCLUDED

#include <iostream>
#include <type_traits>
#include <vector>
#include <functional>
#include <cmath>

namespace MathLib {
    template <typename, typename>
    constexpr bool is_same_template{false}; 

    template <
        template <int, int, typename> class T,
        int l1, int c1, typename A,
        int l2, int c2, typename B
    >
    constexpr bool is_same_template <T<l1, c1, A>, T<l2, c2, B>> {true};
    constexpr static const int MAX_MATRIX_SIZE = 10000;

    template <int row_count, int column_count, typename Type, bool bigMatrix>
    class MatrixContainer {
    public: 
        Type **matrix;

        MatrixContainer() {
            matrix = new Type*[row_count]{0};

            for (int i = 0; i < row_count; i++)
                matrix[i] = new Type[column_count] {Type(0)};
        }

        void freeMem() {
            for (int i = 0; i < row_count; i++) {
                delete [] matrix[i];
            }

            delete [] matrix;
        }

        Type *operator [] (int index) {
            return matrix[index];
        }

        ~MatrixContainer() {
            freeMem();
        }
    };

    template <int row_count, int column_count, typename Type>
    class MatrixContainer <row_count, column_count, Type, false> {
    public: 
        Type matrix[row_count][column_count];

        MatrixContainer() {
            for (int i = 0; i < row_count; i++)
                for (int j = 0; j < column_count; j++)
                    matrix[i][j] = Type(0);
        }

        Type *operator [] (int index) {
            return matrix[index];
        }

        ~MatrixContainer() {}
    };

    template <int rows, int cols, typename Type>
    class MatrixToVectorContainer {
    public:
        using MatCont = MatrixContainer<rows, cols, Type, ((rows * cols) > MAX_MATRIX_SIZE)>;
        MatCont matrix;

        MatrixToVectorContainer() : matrix() {}
        ~MatrixToVectorContainer() {}
    };

    template <typename Type>
    struct MatrixToVectorContainer <1, 1, Type> { 
        using MatCont = MatrixContainer<1, 1, Type, ((1 * 1) > MAX_MATRIX_SIZE)>;
        union {
            MatCont matrix;
            Type array[1];  /// for compatibility with opengl
            union {
                Type x;
                Type r;
            };
        };

        MatrixToVectorContainer() : matrix() {}
        ~MatrixToVectorContainer() {}
    };

    template <typename Type>
    struct MatrixToVectorContainer <2, 1, Type> { 
        using MatCont = MatrixContainer<2, 1, Type, ((2 * 1) > MAX_MATRIX_SIZE)>;
        union {
            struct {
                MatCont matrix;
                Type array[2];  /// for compatibility with opengl
                union {
                    Type x;
                    Type r;
                };
                union {
                    Type y;
                    Type g;
                };
            };
        };

        MatrixToVectorContainer() : matrix() {}
        ~MatrixToVectorContainer() {}
    };

    template <typename Type>
    struct MatrixToVectorContainer <3, 1, Type> { 
        using MatCont = MatrixContainer<3, 1, Type, ((3 * 1) > MAX_MATRIX_SIZE)>;
        union {
            MatCont matrix;
            Type array[3];  /// for compatibility with opengl
            struct {
                union {
                    Type x;
                    Type r;
                };
                union {
                    Type y;
                    Type g;
                };
                union {
                    Type z;
                    Type b;
                };
            };
        };

        MatrixToVectorContainer() : matrix() {}
        ~MatrixToVectorContainer() {}
    };

    template <typename Type>
    struct MatrixToVectorContainer <4, 1, Type> { 
        using MatCont = MatrixContainer<4, 1, Type, ((4 * 1) > MAX_MATRIX_SIZE)>;
        union {
            MatCont matrix;
            Type array[4];  /// for compatibility with opengl
            struct {
                union {
                    Type x;
                    Type r;
                };
                union {
                    Type y;
                    Type g;
                };
                union {
                    Type z;
                    Type b;
                };
                union {
                    Type w;
                    Type a;
                };
            };
        };

        MatrixToVectorContainer() : matrix() {}
        ~MatrixToVectorContainer() {}
    };

    template <int row_count, int column_count, typename Type>
    class Matrix : public MatrixToVectorContainer <row_count, column_count, Type> {
    public:
        constexpr static const int rows = row_count;
        constexpr static const int cols = column_count;

        using MatCont = MatrixToVectorContainer <row_count, column_count, Type>;

        ~Matrix() {}

        template <typename T>
        static constexpr const bool is_matrix{is_same_template<T, Matrix<1, 1, float>>};

        template <typename TypeCols>
        static constexpr int get_col_number() {
            if constexpr (is_matrix<TypeCols>) {
                return TypeCols::cols;
            }
            else {
                return 1;
            }
        }

        /// Mathematical Stuff:
        class MatrixEpsilon {
        public:
            Type epsilon;

            MatrixEpsilon() {
                if constexpr (std::is_arithmetic<Type>::value) {
                    epsilon = std::numeric_limits<Type>::epsilon();
                }
                else {
                    epsilon = Type(0.00001f);
                }
            }

            template <typename Abs_T = double(*)(double)>
            bool areEqual (Type arg1, Type arg2, Abs_T abs = std::abs) {
                return (abs(arg1 - arg2) < epsilon);
            }

            template <typename Abs_T = double(*)(double)>
            bool isZero (Type arg1, Abs_T abs = std::abs) {
                return (abs(arg1 - Type(0)) < epsilon);
            }

            MatrixEpsilon (Type& epsilon) : epsilon(epsilon) {}
            MatrixEpsilon (Type&& epsilon) : epsilon(epsilon) {}
        };

        static MatrixEpsilon defaultEpsilon;

        template <typename Sqrt_T = double(*)(double), typename Abs_T = double(*)(double)>
        Type getFrobeniusNorm (Sqrt_T sqrt = std::sqrt, Abs_T abs = std::abs) {
            Type result = Type(0);

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    result += abs(MatCont::matrix[i][j]) * abs(MatCont::matrix[i][j]);

            return sqrt(result);
        }

        template <typename Sqrt_T = double(*)(double), typename Abs_T = double(*)(double)>
        Type vecNorm2 (Sqrt_T sqrt = std::sqrt, Abs_T abs = std::abs) {
            Type result = Type(0);

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    result += MatCont::matrix[i][j] * MatCont::matrix[i][j];

            return sqrt(result);    
        }

        template <typename Abs_T = double(*)(double)>
        Type vecNorm1 (Abs_T abs = std::abs) {
            Type result = Type(0);

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    result += abs(MatCont::matrix[i][j]);

            return result;  
        }

        template <typename Abs_T = double(*)(double)>
        Type vecNormInf (Abs_T abs = std::abs) {
            Type result = abs(MatCont::matrix[0][0]);

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    if (result < abs(MatCont::matrix[i][j]))
                        result = abs(MatCont::matrix[i][j]);

            return result;  
        }

        template <typename TypeArg>
        Matrix<3, 1, decltype(Type() * TypeArg())> cross (Matrix<3, 1, TypeArg>& arg) {
            Matrix <3, 1, decltype(Type() * TypeArg())> result;

            static_assert((3 == rows && cols == 1), "The vectors must be equal to use cross");

            return Matrix <3, 1, decltype(Type() * TypeArg())> (
                MatCont::y * arg.z - MatCont::z * arg.y,
                MatCont::z * arg.x - MatCont::x * arg.z,
                MatCont::x * arg.y - MatCont::y * arg.x
            );
        }

        template <typename TypeArg>
        Matrix<3, 1, decltype(Type() * TypeArg())> cross (Matrix<3, 1, TypeArg>&& arg) {
            return cross(arg);
        }

        template <int rowsArg, typename TypeArg>
        decltype(Type() * TypeArg()) dot (Matrix<rowsArg, 1, TypeArg>& arg) {
            static_assert((rowsArg == rows && cols == 1), "The vectors must be equal to use dot");

            decltype(Type() * TypeArg()) result = 0;
            for (int i = 0; i < rows; i++)
                result += MatCont::matrix[i][0] * arg[i][0];

            return result;
        }

        template <int rowsArg, typename TypeArg>
        decltype(Type() * TypeArg()) dot (Matrix<rowsArg, 1, TypeArg>&& arg) {
            return dot(arg);
        }

        template <typename Sqrt_T = double(*)(double), typename Abs_T = double(*)(double)>
        MatrixEpsilon getSugestedEpsilon(Sqrt_T sqrt = std::sqrt, Abs_T abs = std::abs) {

            Type norm = getFrobeniusNorm(sqrt, abs);
            if (abs(norm) < 1)
                norm = 1;
            if constexpr (std::is_arithmetic<Type>::value) {
                return MatrixEpsilon(Type(std::numeric_limits<Type>::epsilon()) * norm);
            }
            else {
                return MatrixEpsilon(Type(0.00001f) *  norm);
            }
        }

        template <typename Abs_T = double(*)(double)>
        Type det(MatrixEpsilon &epsilon = defaultEpsilon, Abs_T abs = std::abs) {
            static_assert((rows == cols), "need to have a square matrix to use determinant!");

            Type result = Type(1);
            Type sign = Type(1);
            auto temp = *this;
            for (int k = 0; k < rows; k++) {
                Type maxVal = abs(temp[k][k]);
                std::pair <int, int> pivot(k, k);

                for (int i = k; i < rows; i++)
                    for (int j = k; j < cols; j++)
                        if (abs(temp[i][j]) > maxVal)
                            maxVal = abs(temp[i][j]),
                            pivot = std::pair<int, int>(i, j);

                temp.swapLines(pivot.first, k);
                temp.swapColls(pivot.second, k);

                if (pivot.first != k)
                    sign *= -1;

                if (pivot.second != k)
                    sign *= -1;

                if (epsilon.isZero(temp[k][k]))
                    return Type(0);

                for (int i = k + 1; i < cols; i++)
                    temp[k][i] /= temp[k][k];

                result *= temp[k][k];
                temp[k][k] = Type(1);

                for (int i = k + 1; i < rows; i++) {
                    for (int j = k + 1; j < cols; j++) {
                        temp[i][j] -= temp[k][j] * temp[i][k];
                    }
                    temp[i][k] = 0;
                }
            }
            return result * sign;
        }

        void swapLines (int l1, int l2) {
            for (int i = 0; i < cols; i++)
                std::swap(row(l1)[i], row(l2)[i]);
        }

        void swapColls (int c1, int c2) {
            for (int i = 0; i < rows; i++)
                std::swap(coll(c1)[i], coll(c2)[i]);
        }

        /// Line, column indexers:
    // private:
        class LineIndexer {
        public:
            using MatType = Matrix<1, cols, Type>;

            Matrix <rows, cols, Type>& parentMatrix;
            int line;

        // public:
            LineIndexer(Matrix <rows, cols, Type>& parentMatrix, int line) 
            : parentMatrix(parentMatrix), line(line) {}

        // private: 
            LineIndexer (LineIndexer&& lineIndex) : parentMatrix(lineIndex.parentMatrix) {
                for (int i = 0; i < cols; i++)
                    (*this)[i] = lineIndex[i];
            }

            LineIndexer (LineIndexer& lineIndex) : parentMatrix(lineIndex.parentMatrix) {
                for (int i = 0; i < cols; i++)
                    (*this)[i] = lineIndex[i];
            }

            LineIndexer operator = (LineIndexer& lineIndex) {
                for (int i = 0; i < cols; i++)
                    (*this)[i] = lineIndex[i];

                return (*this);
            }

            LineIndexer operator = (LineIndexer&& lineIndex) {
                for (int i = 0; i < cols; i++)
                    (*this)[i] = lineIndex[i];

                return (*this);
            }

        // public:
            MatType getAsMatrix() {
                return MatType(*this);
            }

            operator MatType () {
                MatType result;

                for (int i = 0; i < cols; i++)
                    result[0][i] = parentMatrix[line][i];

                return result;
            }

            Type& operator [] (int index) {
                return parentMatrix[line][index];
            }

            Type& operator () (int index) {
                return parentMatrix[line][index];
            }       
        };

        class CollumnIndexer {
        public:
            using MatType = Matrix<rows, 1, Type>;

            Matrix <rows, cols, Type>& parentMatrix;
            int collumn;

        // public:
            CollumnIndexer (Matrix <rows, cols, Type>& parentMatrix, int collumn) 
            : parentMatrix(parentMatrix), collumn(collumn) {}

        // private:
            CollumnIndexer (CollumnIndexer&& colIndex) : parentMatrix(colIndex.parentMatrix) {
                for (int i = 0; i < rows; i++)
                    (*this)[i] = colIndex[i];
            }

            CollumnIndexer (CollumnIndexer& colIndex) : parentMatrix(colIndex.parentMatrix) {
                for (int i = 0; i < rows; i++)
                    (*this)[i] = colIndex[i];
            }

            CollumnIndexer operator = (CollumnIndexer& colIndex) {
                for (int i = 0; i < rows; i++)
                    (*this)[i] = colIndex[i];

                return (*this);
            }

            CollumnIndexer operator = (CollumnIndexer&& colIndex) {
                for (int i = 0; i < rows; i++)
                    (*this)[i] = colIndex[i];

                return (*this);
            }

        // public:
            MatType getAsMatrix() {
                return MatType(*this);
            }

            operator MatType () {
                MatType result;

                for (int i = 0; i < rows; i++)
                    result[i][0] = parentMatrix[i][collumn];

                return result;
            }

            Type& operator [] (int index) {
                return parentMatrix[index][collumn];
            }

            Type& operator () (int index) {
                return parentMatrix[index][collumn];
            }
        };

    // public:
        friend LineIndexer;
        friend CollumnIndexer;

        LineIndexer row(int line) {
            return LineIndexer(*this, line);
        }

        CollumnIndexer coll(int coll) {
            return CollumnIndexer(*this, coll);
        }

        Type& operator () (int lin, int col) {
            return MatCont::matrix[lin][col];
        }

        Type& operator () (int index) {
            static_assert (cols == 1, "Can use this operator only on vectors");
            return MatCont::matrix[index][0];
        }

        /// main operators
        template <int rowsB, int colsB, typename TypeB>
        Matrix<rows, cols, decltype(Type() + TypeB())>& operator = (Matrix<rowsB, colsB, TypeB>& mat) {
            static_assert((cols == colsB && rows == rowsB), "Cannot equal, sizes don't match");

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    MatCont::matrix[i][j] = mat[i][j];

            return (*this);
        }

        template <int rowsB, int colsB, typename TypeB>
        Matrix<rows, cols, decltype(Type() + TypeB())>& operator = (Matrix<rowsB, colsB, TypeB>&& mat) {
            static_assert((cols == colsB && rows == rowsB), "Cannot equal, sizes don't match");

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    MatCont::matrix[i][j] = mat[i][j];

            return (*this);
        }

        template <int rowsB, int colsB, typename TypeB>
        Matrix<rows, cols, decltype(Type() + TypeB())> operator + (Matrix<rowsB, colsB, TypeB>& mat) {
            Matrix <rows, cols, decltype(Type() + TypeB())> result;

            static_assert((cols == colsB && rows == rowsB), "Cannot add, sizes don't match");

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    result[i][j] = MatCont::matrix[i][j] + mat[i][j];

            return result;
        }

        template <int rowsB, int colsB, typename TypeB>
        Matrix<rows, cols, decltype(Type() - TypeB())> operator - (Matrix<rowsB, colsB, TypeB>& mat) {
            Matrix <rows, cols, decltype(Type() - TypeB())> result;

            static_assert((cols == colsB && rows == rowsB), "Cannot substract, sizes don't match");

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    result[i][j] = MatCont::matrix[i][j] - mat[i][j];

            return result;
        }

        template <int rowsB, int colsB, typename TypeB>
        Matrix<rows, colsB, decltype(Type() * TypeB())> operator * (Matrix<rowsB, colsB, TypeB>& mat) {
            Matrix <rows, colsB, decltype(Type() * TypeB())> result;

            static_assert((cols == rowsB), "Cannot multiply, sizes don't match");

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < colsB; j++)
                    for (int k = 0; k < cols; k++)
                        result[i][j] += MatCont::matrix[i][k] * mat[k][j];

            return result;
        }

        Matrix<rows, cols, Type> operator - () {
            return (*this) * Type(-1);
        }

        template <typename ScalarType>
        Matrix<rows, cols, decltype(Type() * ScalarType())> operator * (ScalarType& scalar) {
            Matrix<rows, cols, decltype(Type() * ScalarType())> result;

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    result[i][j] = MatCont::matrix[i][j] * scalar;

            return result;
        }

        template <typename ScalarType>
        Matrix<rows, cols, decltype(Type() / ScalarType())> operator / (ScalarType& scalar) {
            Matrix<rows, cols, decltype(Type() / ScalarType())> result;

            static_assert((scalar != 0), "Divide by zero");

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    result[i][j] = MatCont::matrix[i][j] / scalar;

            return result;
        }

        /// friends, references and rvalues
        template <typename ScalarType>
        friend Matrix<rows, cols, decltype(Type() * ScalarType())> operator * 
                (ScalarType& scalar, Matrix<rows, cols, decltype(Type() * ScalarType())>& mat) 
        {
            return (mat * scalar);
        }

        template <typename ScalarType>
        friend Matrix<rows, cols, decltype(Type() * ScalarType())> operator * 
                (ScalarType& scalar, Matrix<rows, cols, decltype(Type() * ScalarType())>&& mat) 
        {
            return (mat * scalar);
        }

        template <typename ScalarType>
        friend Matrix<rows, cols, decltype(Type() * ScalarType())> operator * 
                (ScalarType&& scalar, Matrix<rows, cols, decltype(Type() * ScalarType())>& mat) 
        {
            return (mat * scalar);
        }

        template <typename ScalarType>
        friend Matrix<rows, cols, decltype(Type() * ScalarType())> operator * 
                (ScalarType&& scalar, Matrix<rows, cols, decltype(Type() * ScalarType())>&& mat) 
        {
            return (mat * scalar);
        }

        template <typename ScalarType>
        Matrix<rows, cols, decltype(Type() / ScalarType())> operator / (ScalarType&& scalar) {
            return (*this / scalar);
        }

        template <typename ScalarType>
        Matrix<rows, cols, decltype(Type() * ScalarType())> operator * (ScalarType&& scalar) {
            return (*this * scalar);
        }

        template <int rowsB, int colsB, typename TypeB>
        Matrix<rows, colsB, decltype(Type() * TypeB())> operator * (Matrix<rowsB, colsB, TypeB>&& mat) {
            return (*this * mat);
        }

        template <int rowsB, int colsB, typename TypeB>
        Matrix<rows, cols, decltype(Type() - TypeB())> operator - (Matrix<rowsB, colsB, TypeB>&& mat) {
            return (*this - mat);
        }

        template <int rowsB, int colsB, typename TypeB>
        Matrix<rows, cols, decltype(Type() + TypeB())> operator + (Matrix<rowsB, colsB, TypeB>&& mat) {
            return (*this + mat);
        }

        /// operator <something>=
        template <int rowsB, int colsB, typename TypeB>
        Matrix<rows, cols, decltype(Type() + TypeB())> operator += (Matrix<rowsB, colsB, TypeB>& mat) {
            return ((*this) = (*this) + mat);
        }

        template <int rowsB, int colsB, typename TypeB>
        Matrix<rows, cols, decltype(Type() - TypeB())> operator -= (Matrix<rowsB, colsB, TypeB>& mat) {
            return ((*this) = (*this) - mat);
        }

        template <int rowsB, int colsB, typename TypeB>
        Matrix<rows, colsB, decltype(Type() * TypeB())> operator *= (Matrix<rowsB, colsB, TypeB>& mat) {
            return ((*this) = (*this) * mat);
        }

        template <typename ScalarType>
        Matrix<rows, cols, decltype(Type() * ScalarType())> operator *= (ScalarType& scalar) {
            return ((*this) = (*this) * scalar);
        }

        template <typename ScalarType>
        Matrix<rows, cols, decltype(Type() / ScalarType())> operator /= (ScalarType& scalar) {
            return ((*this) = (*this) / scalar);
        }

        Type* operator [] (int index) {
            return MatCont::matrix[index];
        }

        Matrix<cols, rows, Type> tr() {
            Matrix<cols, rows, Type> result;

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    result[j][i] = MatCont::matrix[i][j];

            return result;
        }

        // The matrix will be constructed as follows: 
        // -> we will consider basic types to be a matrix of 1 X 1
        // -> we will take the first matrix and place it at (0, 0)
        // -> if the second matrix has space on the right of the first matrix then we place it on the right
        // -> else we jump on the line under the matrix we just placed at collumn 0 and 
        // continue placing the next matrix from there

        Matrix () {}

        template <typename ArgType, typename ...Args>
        Matrix (ArgType& arg, Args ...args) {
            fill_mat <0, 0, is_matrix<ArgType>, ArgType, Args...> (arg, args...);
        }

        template <typename ArgType, typename ...Args>
        Matrix (ArgType&& arg, Args ...args) {
            fill_mat <0, 0, is_matrix<ArgType>, ArgType, Args...> (arg, args...);
        }

        template <int lin, int col, bool is_matrix_val, typename ArgType, typename NextType, typename... Args>
        typename std::enable_if<!is_matrix_val, void>::type fill_mat(ArgType& arg, NextType& nextArg, Args ...args) {
            MatCont::matrix[lin][col] = arg;

            if constexpr (col + 1 + get_col_number<NextType>() <= cols) // we have space for the next matrix
                fill_mat <lin, col + 1, is_matrix<NextType>, NextType, Args...> (nextArg, args...);
            else {                                                      // we don't have space for the matrix
                fill_mat <lin + 1, 0, is_matrix<NextType>, NextType, Args...> (nextArg, args...);
            }
        }

        template <int lin, int col, bool is_matrix_val, typename ArgType, typename NextType, typename... Args>
        typename std::enable_if<is_matrix_val, void>::type fill_mat(ArgType& arg, NextType& nextArg, Args ...args) {
            for (int i = 0; i < ArgType::rows; i++)
                for (int j = 0; j < ArgType::cols; j++)
                        MatCont::matrix[lin + i][col + j] = arg[i][j];

            if constexpr (col + ArgType::cols + get_col_number<NextType>() <= cols) {
                fill_mat <lin, col + ArgType::cols, is_matrix<NextType>, NextType, Args...> (nextArg, args...);
            }
            else {                                              // we don't have space for the matrix
                fill_mat <lin + ArgType::rows, 0, is_matrix<NextType>, NextType, Args...> (nextArg, args...);
            }
        }

        template <int lin, int col, bool is_matrix_val, typename ArgType, typename NextType>
        typename std::enable_if<!is_matrix_val, void>::type fill_mat(ArgType& arg, NextType& nextArg) {
            MatCont::matrix[lin][col] = arg;

            if constexpr (col + 1 + get_col_number<NextType>() <= cols) {
                fill_mat <lin, col + 1, is_matrix<NextType>, NextType> (nextArg);
            }
            else {                                              // we don't have space for the matrix
                fill_mat <lin + 1, 0, is_matrix<NextType>, NextType> (nextArg);
            }
        }

        template <int lin, int col, bool is_matrix_val, typename ArgType, typename NextType>
        typename std::enable_if<is_matrix_val, void>::type fill_mat(ArgType& arg, NextType& nextArg) {
            for (int i = 0; i < ArgType::rows; i++)
                for (int j = 0; j < ArgType::cols; j++)
                        MatCont::matrix[lin + i][col + j] = arg[i][j];

            if constexpr (col + ArgType::cols + get_col_number<NextType>() <= cols) {
                fill_mat <lin, col + ArgType::cols, is_matrix<NextType>, NextType> (nextArg);
            }
            else {                                              // we don't have space for the matrix
                fill_mat <lin + ArgType::rows, 0, is_matrix<NextType>, NextType> (nextArg);
            }
        }

        template <int lin, int col, bool is_matrix_val, typename ArgType>
        typename std::enable_if<!is_matrix_val, void>::type fill_mat(ArgType& arg) {
            MatCont::matrix[lin][col] = arg;
        }

        template <int lin, int col, bool is_matrix_val, typename ArgType>
        typename std::enable_if<is_matrix_val, void>::type fill_mat(ArgType& arg) {
            for (int i = 0; i < ArgType::rows; i++)
                for (int j = 0; j < ArgType::cols; j++)
                        MatCont::matrix[lin + i][col + j] = arg[i][j];
        }

        /// ostream, istream:
        friend std::ostream& operator << (std::ostream& stream, Matrix<rows, cols, Type>& arg) {
            stream << "rows: " << rows << ", cols:" << cols << std::endl; 
            for (int i = 0; i < rows; i++) {
                for (int j = 0; j < cols; j++) {
                    stream << arg[i][j] << " ";
                }
                if (i < rows - 1)
                    stream << std::endl;
            }

            return stream;
        }

        friend std::ostream& operator << (std::ostream& stream, Matrix<rows, cols, Type>&& arg) {
            stream << "rows: " << rows << ", cols:" << cols << std::endl;
            for (int i = 0; i < rows; i++) {
                for (int j = 0; j < cols; j++) {
                    stream << arg[i][j] << " ";
                }
                if (i < rows - 1)
                    stream << std::endl;
            }

            return stream;
        }

        friend std::istream& operator << (std::istream& stream, Matrix<rows, cols, Type>& arg) {
            for (int i = 0; i < rows; i++) {
                for (int j = 0; j < cols; j++) {
                    stream >> arg[i][j];
                }
            }

            return stream;
        }
    };

    template <int rows, int cols, typename Type>
    typename Matrix<rows, cols, Type>::MatrixEpsilon Matrix<rows, cols, Type>::defaultEpsilon;
}
#endif

MatrixHelper.h:

#ifndef MATRIX_HELPER_H
#define MATRIX_HELPER_H

#include "Matrix.h"

namespace MathLib {
    template <int size, typename Type>
    using Vector = Matrix <size, 1, Type>;

    template <int x, int y>
    using Matd = Matrix <x, y, double>;

    template <int x, int y>
    using Matf = Matrix <x, y, float>;

    template <int x, int y>
    using Mati = Matrix <x, y, int>;

    template <int x>
    using Vecd = Vector <x, double>;

    template <int x>
    using Vecf = Vector <x, float>;

    template <int x>
    using Veci = Vector <x, int>;

    using Vec3d = Vecd<3>;
    using Vec3f = Vecf<3>;
    using Vec3i = Veci<3>;

    using Point3d = Vec3d;
    using Point3f = Vec3f;
    using Point3i = Vec3i;

    using Vec4d = Vecd<4>;
    using Vec4f = Vecf<4>;
    using Vec4i = Veci<4>;

    using Point4d = Vec4d;
    using Point4f = Vec4f;
    using Point4i = Vec4i;

    using Vec2d = Vecd<2>;
    using Vec2f = Vecf<2>;
    using Vec2i = Veci<2>;

    using Point2d = Vec2d;
    using Point2f = Vec2f;
    using Point2i = Vec2i;

    using Matrix2d = Matd<2, 2>;
    using Matrix2f = Matf<2, 2>;
    using Matrix2i = Mati<2, 2>;

    using Matrix3d = Matd<3, 3>;
    using Matrix3f = Matf<3, 3>;
    using Matrix3i = Mati<3, 3>;

    using Matrix4d = Matd<4, 4>;
    using Matrix4f = Matf<4, 4>;
    using Matrix4i = Mati<4, 4>;

    template <int rows, typename Type>
    Matrix <rows, rows, Type> Identity () {
        Matrix <rows, rows, Type> mat;

        for (int i = 0; i < rows; i++)
            mat[i][i] = Type(1);

        return mat;
    }
}

#endif

Here is a very small usage example:

#include <iostream>
#include "MathLib.h"

void testBasic();
// int testBasicOld();

int main(int argc, char const *argv[])
{
    using namespace MathLib;
    using namespace std;

    testBasic();
}

void testBasic() {
    using namespace MathLib;

    Matrix4f mat1(
        1, 0, 0, 0,
        0, 2, 0, 0,
        0, 0, 2, 0,
        0, 0, 0, 1);

    Vec4f vec1(0, 1, 0, 0);
    Vec4f vec2(0, 1, 1, 0);

    Matrix4f mat2 = mat1;

    std::cout << vec1.x << std::endl;
    std::cout << vec1.y << std::endl;
    std::cout << vec1.z << std::endl;

    std::cout << mat1 << std::endl;
    std::cout << mat2 << std::endl;

    std::cout << vec1 * vec2.tr() << std::endl;
    std::cout << vec2.tr() * vec1 << std::endl;
}

Answer 1

Given that there's around a thousand lines of code to review, I trust you're not expecting anything too in depth.

Almost all of the concerns I have with this code stem from what I suspect is the desire to have certain interface characteristics. In other words, I suspect that a lot of your motivation in this code is that you want to be able to write C++ code that looks exactly like mathematical notation, or some other style. Unfortunately, I think you're going to have to give a lot of that up.

MatrixContainer

Let's start with the "large" MatrixContainer class. The core of the class is this:

Type **matrix;

MatrixContainer() {
    matrix = new Type*[row_count]{0};

    for (int i = 0; i < row_count; i++)
        matrix[i] = new Type[column_count] {Type(0)};
}

Now I suspect the intention here was that you want to be able to address elements in your matrix as m[row][col]. Frankly, I think you're going to have to give that dream up. At least for now. That's because, in your efforts to get that notation, this class has become over-complicated, bloated, and inefficient. And also unsafe. And since it serves as the base of your entire matrix hierarchy, the whole thing is built on quicksand.

Let's start with the default constructor. Imagine you're constructing a 500×200 element matrix. The first line of the constructor will allocate an array of 500 pointers. Then the loop will allocate 500 200-element arrays. Now image the loop is chugging along, and it allocates 499 200-element arrays no problem... then the 500th allocation fails. std::bad_alloc is thrown. The stack unwinds. And your constructor has just leaked 500 pointers and 99800 numbers in memory.

Now there are several techniques you could use to fix that leak, but they're all ugly. They also won't help fix all your problems.

So let's say you do manage to safely allocate and construct your 500×200 element matrix. Great! Now you have a valid matrix object m1. Everything's cool. But then you do auto m2 = m1. Unfortunately, due to historical compatibility with C, your class has an implicit copy constructor. The implicit copy constructor just copies all the (non-static) data members. You have only one: matrix, a pointer (to an array of pointers to arrays of numbers). So the pointer gets copied. At first, you may notice no problems. Then you may notice that when you make changes to m1, those changes also appear in m2. Hmm. But then the real disaster happens when you finally destroy one of the objects. Let's say m2 dies first. Your destructor is correctly written, and will properly clean up all the memory you allocated. And then... m1 dies. (Or not. Even if you just try to use it, you'll probably trigger an access violation.) Its destructor also tries to run... but the memory has already been freed. Boom. Double-delete. Program crash.

This too is fixable, but if you thought the code to fix the default constructor was ugly, you ain't seen nothing yet. And you need not only a copy constructor, but also copy assignment, and probably move constructor and assignment, too. All of them will be complicated and ugly.

At this point I'd say stop. It's time to rethink the design.

Let's give up that nice, textbook-pretty notation. Instead of arrays of arrays, let's just allocate a flat array, and calculate multidimensional offsets. In other words, instead of m[row][col], let's have m[row * column_count + col]. We can hide that extra ugliness in a function, so we could write, for example, m(row, col). That ain't bad. And you already have that interface in your matrix class.

There's another reason to prefer a flat array to an array of arrays, and that has to do with efficiency. A flat array uses less memory, and plays better with caching, which means it can be orders of magnitude faster for common operations.

What would the class look like with a flat array? Maybe something like this:

template <int Rows, int Cols, typename Type, bool bigMatrix>
class MatrixContainer {
public: 
    static_assert(Rows >= 1);
    static_assert(Cols >= 1);

    MatrixContainer() :
        matrix_{new Type[Rows * Cols]{0}}
    {}

    // Copy constructor
    MatrixContainer(MatrixContainer const& m) :
        matrix_{new Type[Rows * Cols]{0}}
    {
        std::copy_n(m.matrix_, Rows * Cols, matrix_);
    }

    // Move constructor
    MatrixContainer(MatrixContainer&& m) noexcept
    {
        swap_(m);
    }

    ~MatrixContainer() {
        delete[] matrix_;
    }

    Type& operator()(int row, int col) noexcept {
        return matrix_[row * Cols + col];
    }

    Type const& operator()(int row, int col) const noexcept {
        return matrix_[row * Cols + col];
    }

    // Copy assign
    MatrixContainer& operator=(MatrixContainer const& m)
    {
        // Copy and swap technique.
        auto temp = m;
        swap_(temp);
        return *this;
    }

    // Move constructor
    MatrixContainer& operator=(MatrixContainer&& m) noexcept
    {
        swap_(m);
        return *this;
    }

private:
    void swap_(MatrixContainer& other) noexcept
    {
        std::swap(matrix_, other.matrix_);
    }

    Type* matrix_ = nullptr;
};

That's quite a bit bigger than what you've got, but it's also a lot simpler. Most of the functions are boilerplate - even swap() could be made part of the interface, instead of an implementation detail. All the loops are gone. The free memory function is gone.

Can this be improved? It can. It turns out there is already a class in the standard library for dealing with a dynamically allocated array. It's called vector.

That would give us this:

template <int Rows, int Cols, typename Type, bool bigMatrix>
class MatrixContainer {
public: 
    static_assert(Rows >= 1);
    static_assert(Cols >= 1);

    MatrixContainer() :
        matrix_(Rows * Cols)
    {}

    Type& operator()(int row, int col) noexcept {
        return matrix_[row * Cols + col];
    }

    Type const& operator()(int row, int col) const noexcept {
        return matrix_[row * Cols + col];
    }

private:
    std::vector<Type> matrix_;
};

No joke: You cannot write a more efficient matrix class than this. (Assuming you can't take advantage of small-vector optimizations, which vector can't do.) And this class doesn't leak, or trigger crashes when you try to copy or move it. Even if Type throws, this will still work properly. That function call operator will probably be reduced to a single instruction or two.

If you're not convinced yet to give up the "array of array of numbers" design in favour of the flat array design, then let's put a pin in that. We'll get back to it.

For now, let's move on to the "small" MatrixContainer.

This class is much smaller and simpler, but it still has nested loops. Here's what it might look like with a flat design:

template <int Rows, int Cols, typename Type, bool bigMatrix>
class MatrixContainer<Rows, Cols, Type, false> {
public: 
    static_assert(Rows >= 1);
    static_assert(Cols >= 1);

    MatrixContainer()
    {
        std::fill_n(matrix_, Rows * Cols, T(0));
    }

    Type& operator()(int row, int col) noexcept {
        return matrix_[row * Cols + col];
    }

    Type const& operator()(int row, int col) const noexcept {
        return matrix_[row * Cols + col];
    }

private:
    std::array<Type, Rows * Cols> matrix_;
};

Almost exactly the same as the other version, except that it uses an array rather than a vector. If you really want to use a C-array, you could.

However, I would recommend against using this "small" version for anything that's not really small. The reason is that when you use an array, you give up all the wonderful benefits of move semantics. The "large" class is more expensive to construct and copy... but those are things you avoid anyway when writing performant code. The real cost of your class is probably going to be in moves and swaps. And vector is trivial to move and swap - almost free! (It's like 3 pointer copies/swaps.) array, however, cannot be moved. And it cannot be cheaply swapped - it must always be swapped element by element.

That means that if your workload is almost all moves and swaps - which it often is in hot code segments - vector will probably be faster than array for 4 or more elements.

Of course, you'll have to profile your code to know for sure if this workload is what you're dealing with.

Before I move on from MatrixContainer, there's one last thing. To get this increase in efficiency, simplicity, and reliability, I told you you'd have to give up your pretty m[row][col] syntax. I lied. Even with a flat design - with all its benefits - you can still have that syntax. Unfortunately, getting it is not easy, and it requires proxy objects and other arcane tricks. You already seem to have at least part of the solution with your line and column indexers - though I'm not sure because there's just too much code to review, so I didn't really look into what those do. But the bottom line is that this is a situation where you can have your cake and eat it too.

So I would strongly suggest giving up the array-of-arrays design, and going with a flat design.

MatrixToVectorContainer

As far as I can tell, the entire purpose of this family of classes is to create synonyms to access certain elements in certain cases. Like, when you have a Matrix<3, 1, int>, you want to be able to access m[0][0] (or m(0,0)) as m.x or m.r.

Unfortunately, the method you're using to do it is not legal C++.

Using unions this way, known as type-punning, is legal C, but not legal C++. (Most if not all compilers have extensions that will allow you to get away with it, though.) In C++, the only union member you are allowed to read from is the last one you wrote to. In other words, you could do m.r = 1; cout << m.r; ... but you could not do m.r = 1; cout << m.x; // or m[0][0] or m(0,0).

Once again, if you're willing to give up your dream syntax, you can get what you want.

In this case, instead of m.x, you could have m.x(). Otherwise, everything else is the same.

All you'd need to do is forget about all those unions. Just get rid of them all. (You don't even need them for the OpenGL type cast, if you're using a flat array design as I recommend.) Then just replace them with standard accessors. For example:

Type& x() noexcept {
    return matrix_[0]; // or (*this)(0, 0);
}

Type const& x() const noexcept {
    return matrix_[0];
}

Once you're doing this, the MatrixToVectorContainer classes all become unnecessary. What you'd do instead is use SFINAE/enable_if in your Matrix class to enable x(), g(), or whatever depending on the Rows and Cols template parameters.

Matrix

And now to the monster, which there's no way I can get even halfway deep into it.

I'll skip the is_matrix and get_col_number templates, because I don't see a point to them, but maybe there's a purpose a hundred or three lines later.

I don't see the purpose of the MatrixEpsilon class. It seems to me that all you need is static Type defaultEpsilon = is_arithmetic_v<Type> ? numeric_limits<Type>::epsilon() : Type(0.00001f);. I certainly don't see the point of the copy/move constructors - in fact, I'm surprised they compile.

The main problem in the MatrixEpsilon class is that you do this:

template <typename Abs_T = double(*)(double)>
bool isZero (Type arg1, Abs_T abs = std::abs)

The reason that's a problem is that there's no guarantee that std::abs has the signature double (double). I know that sounds crazy, but its true. The standard works on an "as if" basis; std::abs has to be "as if" it has the signature double (double) when you call it. That could mean it actually has a hidden default parameter, or it might have some bizarre calling convention.

The upshot is that you cannot assume that you can take the address of (the right overload of) std::abs, and put it in a double (*)(double).

What you'd need to do instead is this:

template <typename Abs_T>
bool isZero(Type arg1, Abs_T&& abs) {
    return (abs(arg1 - Type(0)) < epsilon);
}

bool isZero(Type arg1) {
    return isZero(arg1, [](auto&& x){ return std::abs(x); });
}

The same holds true for all cases of taking a standard library function's address.

Next is a bunch of computation functions. These all use nested loops, and could probably be massively simplified with a flat design. Just the first example, getFrobeniusNorm() could become:

template <typename Sqrt_T, typename Abs_T>
Type getFrobeniusNorm (Sqrt_T&& sqrt, Abs_T&& abs) {
    return sqrt(
        std::transform_reduce(
            matrix_.begin(), matrix.end(), T{},
            [&](auto&& x) { return abs(x) * abs(x); }));
}

Type getFrobeniusNorm () {
    return getFrobeniusNorm(
        [](auto&& x) { return std::sqrt(x); },
        [](auto&& x) { return std::abs(x); });
}

As a bonus, when you can use an algorithm (like transform_reduce()), you can also make it parallelized, or even vectorized.

A lot of the calculation functions could easily be moved out of the class. The cross product, for example, seems like a natural candidate for putting outside. You wouldn't need the static assert then. (Also, you declare result unnecessarily.)

template <typename T, typename U>
Matrix<3, 1, decltype(std::declval<T>() * std::declval<U>())>
cross(Matrix<3, 1, T> const& t, Matrix<3, 1, U> const& u)
{
    return Matrix <3, 1, decltype(std::declval<T>() * std::declval<U>())>{
        t.y() * u.z() - t.z() * u.y(),
        t.z() * u.x() - t.x() * u.z(),
        t.x() * u.y() - t.y() * u.x()
    };
}

The row and line indexers are next. I'm going to skip them.

Next is the assignment operators. Both can be defaulted now that the MatrixContainer handles that correctly.

Next is the other arithmetic operators. Most of these can and should be moved out of the class. With a flat design, most can be written using simple algorithms, like std::transform().

Next are some constructors and fill_mat() functions. I didn't look too deeply into these, but it seems over-complex and brittle. For example, if I'm building a 3×3 matrix, and the constructor arguments are a 2×2 matrix {{a,b},{c,d}}, then 5 values e, f, g, h, and i... what do i get in the final matrix? By my reading of the intention, it sounds like I get {{a,b,e},{f,g,h},{i,0,0}} - the values c and d just got nuked.

If there's really a use case for building a matrix hodgepodge like that, it should probably be a factory function, not a constructor. Constructors should be reserved for basic construction of the object - for the obvious and necessary ways of construction. For a N×M matrix, the only constructors I need are default (to zero initialize the elements), copy from another N×M matrix, efficient move from another N×M matrix, and initialize from up to N×M values (setting the unset values to zero). That's it. If I have a bunch of static factory functions for useful things like identity matrices and so on, that's a nice bonus.

And finally stream insertion/extraction.

The extraction looks fine (though, as always, with a flat design it becomes much simpler), but I don't see the purpose for two insertion functions. Especially since nothing is ever moved in the rvalue-reference version. Normally all you need is a single insertion function taking const-ref, and you're done. (I'd also replace all the std::endls with simple newlines. Granted, you probably don't care much about efficiency in your stream ops, but why make them hundreds of times slower unnecessarily?)

Summary

That was a fairly rushed review, but if there's a single takeaway here, I'd say that it's to use a flat architecture, and accept that the usage syntax won't be exactly what you want. (Though if you're really finicky about the syntax, you can get it, mostly. It's just a lot of pain.) That should allow you to vastly simplify everything, and once everything is simplified, it will be easier to spot optimization opportunities.

Answer 2

Data layout

A first impression gives me some concerns about your data representation. Firstly, in modern C++, we can avoid the responsibility for ownership of resources - instead, we can use (without performance overhead) objects such as smart pointers and standard containers to manage that ownership for you, and so avoid having to write destructors and (often) copy constructors and assignment operators.

Secondly, the array-of-arrays storage can give poor locality of your data, as well as taking more time to create and destroy. A more cache-friendly representation is to create a flat array (probably as a std::vector) and provide an operator()(std::size_t row, std::size_t col) to access the elements. (In fact, if you're converting code, you can start by implementing that method; after the code is updated to use it, it's then much easier to change the representation - that's the beauty of abstraction). Using a standard collection such as std::vector will (assuming a decent Standard Library implementation) also take care of the "small object" optimisation you've laboriously created.

Anonymous structs

C++ does not permit anonymous structs, so this code gives errors:

    union {
        MatCont matrix;
        Type array[3];  /// for compatibility with opengl
        struct {
            union {
                Type x;
                Type r;
            };
            union {
                Type y;
                Type g;
            };
            union {
                Type z;
                Type b;
            };
        };

The preceding one is even weirder - a union of a single anonymous struct and nothing else:

    union {
        struct {
            MatCont matrix;
            Type array[2];  /// for compatibility with opengl
            union {
                Type x;
                Type r;
            };
            union {
                Type y;
                Type g;
            };
        };
    };

Unions containing class types are dangerous anyway - we have std::variant for those.

Naming/typos

std::istream& operator <<()

We conventionally call that operator >>, to reduce confusion with the output operator.

Throughout, "columns" is misspelt with two l's.