mat.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_MATRIX_H
#define THEORETICA_MATRIX_H
 
#ifndef THEORETICA_NO_PRINT
#include <sstream>
#include <ostream>
#endif
 
#include <array>
#include <vector>
 
#include "../core/error.h"
#include "../core/constants.h"
#include "../core/real_analysis.h"
#include "./algebra.h"
#include "./transform.h"
#include "./vec.h"
 
 
namespace theoretica {
 
 
    template<typename Matrix, typename ReturnType = matrix_element_t<Matrix>&>
    class mat_iterator {
 
        private:
 
            Matrix& matrix;
 
            size_t row;
 
            size_t col;
 
        public:
 
            using iterator_category = std::forward_iterator_tag;
            using value_type = matrix_element_t<Matrix>;
            using pointer = value_type*;
            using reference = value_type&;
 
            mat_iterator(
                Matrix& matrix,
                size_t row = 0,
                size_t col = 0)
            : matrix(matrix), row(row), col(col) {}
 
 
            ReturnType operator*() {
                return matrix(row, col);
            }
 
 
            mat_iterator& operator++() {
 
                col++;
 
                if(col == matrix.cols()) {
                    col = 0;
                    row++;
                }
 
                return *this;
            }
 
 
            size_t row_index() {
                return row;
            }
 
 
            size_t col_index() {
                return col;
            }
 
 
            // mat_iterator& operator--() {
            //  if(row == 0) {
            //      TH_MATH_ERROR(
            //          "iterator::operator--",
            //          row, IMPOSSIBLE_OPERATION);
            //      return *this;
            //  }
            //  if(col == 0) {
            //      col = max_col;
            //      row--;
            //  }
            //  col--;
            //  return *this;
            // }
 
 
            // Comparison operators
 
            bool operator==(const mat_iterator& other) const {
                return (row == other.row) &&
                    (col == other.col);
            }
 
            bool operator!=(const mat_iterator& other) const {
                return !(*this == other);
            }
    };
 
 
    template<typename Type = real, unsigned int N = 0, unsigned int K = 0>
    class mat {
        public:
 
#ifdef THEORETICA_ROW_FIRST
        Type data[N][K];
#else
        Type data[K][N];
#endif
 
 
        mat() {
            algebra::mat_zeroes(*this);
        }
 
 
        template<typename Matrix>
        mat(const Matrix& m) {
            algebra::mat_copy(*this, m);
        }
 
 
        template<typename T = Type>
        inline mat(const std::initializer_list<std::initializer_list<T>>& rows) {
 
            if(rows.size() != N) {
                TH_MATH_ERROR("mat::mat", rows.size(), INVALID_ARGUMENT);
                algebra::mat_error(*this);
                return;
            }
 
            unsigned int i = 0;
            unsigned int j = 0;
 
            for (const auto& row : rows) {
 
                if (row.size() != K) {
                    TH_MATH_ERROR("mat::mat", rows.size(), INVALID_ARGUMENT);
                    algebra::mat_error(*this);
                    return;
                }
 
                for (const auto& x : row) {
                    
                    at(i, j) = x;
                    j = (j + 1) % K;
                }
 
                i++;
            }
        }
 
 
        template<typename Matrix>
        inline mat<Type, N, K>& operator=(const Matrix& other) {
            return algebra::mat_copy(*this, other);
        }
 
 
        inline void make_zeroes() {
            algebra::mat_zeroes(*this);
        }
 
 
        inline static mat<Type, N, K> zeroes() {
            mat<Type, N, K> res;
            algebra::mat_zeroes(res);
            return res;
        }
 
 
        template<typename Matrix>
        inline mat<Type, N, K> operator+(const Matrix& other) const {
            mat<Type, N, K> res;
            return algebra::mat_sum(res, *this, other);
        }
 
 
        template<typename Matrix>
        inline mat<Type, N, K> operator-(const Matrix& other) const {
            mat<Type, N, K> res;
            return algebra::mat_diff(res, *this, other);
        }
 
 
        inline mat<Type, N, K> operator*(Type scalar) const {
            mat<Type, N, K> res;
            return algebra::mat_scalmul(res, scalar, *this);
        }
 
 
        inline friend mat<Type, N, K> operator*(Type a, const mat<Type, N, K>& B) {
            return B * a;
        }
 
 
        template<typename VecType, unsigned int M>
        inline friend vec<VecType, K> operator*(
            const vec<VecType, M>& a, const mat<Type, N, K>& B) {
            return B * a;
        }
 
 
        inline mat<Type, N, K> operator/(Type scalar) const {
 
            mat<Type, N, K> res;
 
            if(abs(scalar) < MACH_EPSILON) {
                TH_MATH_ERROR("mat::operator/", scalar, DIV_BY_ZERO);
                return algebra::mat_error(res);
            }
            
            return algebra::mat_scalmul(res, 1.0 / scalar, *this);
        }
 
 
        template<typename Vector>
        inline Vector transform(const Vector& v) const {
 
            if(v.size() != cols()) {
                TH_MATH_ERROR("mat::transform", v.size(), INVALID_ARGUMENT);
                Vector res;
                res.resize(cols());
                algebra::vec_error(res);
                return res;
            }
 
            return algebra::transform(*this, v);
        }
 
 
        inline vec<Type, N> transform(const vec<Type, K>& v) const {
            return algebra::transform(*this, v);
        }
 
 
        inline vec<Type, N> operator*(const vec<Type, K>& v) const {
            return transform(v);
        }
 
 
        template<unsigned int M>
        inline mat<Type, N, M> mul(const mat<Type, K, M>& B) const {
            mat<Type, N, M> res;
            return algebra::mat_mul(res, *this, B);
        }
 
 
        template<typename Matrix>
        inline Matrix mul(const Matrix& B) const {
 
            Matrix res;
            res.resize(N, B.cols());
 
            if(B.rows() != K) {
                TH_MATH_ERROR("mat::transform", B.rows(), INVALID_ARGUMENT);
                algebra::mat_error(res);
                return res;
            }
 
            return algebra::mat_mul(res, *this, B);
        }
 
 
        template<typename Matrix>
        inline auto operator*(const Matrix& B) const {
 
            Matrix res;
            res.resize(N, B.cols());
 
            if(B.rows() != K) {
                TH_MATH_ERROR("mat::transform", B.rows(), INVALID_ARGUMENT);
                algebra::mat_error(res);
                return res;
            }
 
            return algebra::mat_mul(res, *this, B);
        }
 
 
        template<typename Matrix>
        inline mat<Type, N, K>& operator+=(const Matrix& other) {
            return algebra::mat_sum(*this, other);
        }
 
 
        template<typename Matrix>
        inline mat<Type, N, K>& operator-=(const Matrix& other) {
            return algebra::mat_diff(*this, other);
        }
 
 
        inline mat<Type, N, K>& operator*=(Type scalar) {
            return algebra::mat_scalmul(scalar, *this);
        }
 
 
        inline mat<Type, N, K>& operator/=(Type scalar) {
 
            if(abs(scalar) < MACH_EPSILON) {
                TH_MATH_ERROR("mat::operator/", scalar, DIV_BY_ZERO);
                return algebra::mat_error(*this);
            }
 
            return algebra::mat_scalmul(1.0 / scalar, *this);
        }
 
 
        template<typename Matrix>
        inline mat<Type, N, K>& operator*=(const Matrix& B) {
            return (*this = this->operator*(B));
        }
 
 
        inline mat<Type, N, K>& transpose() {
            static_assert(
                N == K, "The matrix must be square to be transposed in place.");
            return algebra::make_transposed(*this);
        }
 
 
        inline mat<Type, K, N> transposed() const {
            return algebra::transpose<mat<Type, N, K>, mat<Type, K, N>>(*this);
        }
 
 
        inline Type& at(unsigned int i, unsigned int j) {
 
#ifdef THEORETICA_ROW_FIRST
            return data[i][j];
#else
            return data[j][i];
#endif
        }
 
 
        inline const Type& at(unsigned int i, unsigned int j) const {
 
#ifdef THEORETICA_ROW_FIRST
            return data[i][j];
#else
            return data[j][i];
#endif
        }
 
 
        inline Type& operator()(unsigned int i, unsigned int j) {
            return at(i, j);
        }
 
 
        inline const Type& operator()(unsigned int i, unsigned int j) const {
            return at(i, j);
        }
 
 
        inline Type get(unsigned int i, unsigned int j) const {
 
#ifdef THEORETICA_ROW_FIRST
            return data[i][j];
#else
            return data[j][i];
#endif
        }
 
 
        using iterator = mat_iterator<mat<Type, N, K>, Type&>;
 
 
        using const_iterator = mat_iterator<const mat<Type, N, K>, const Type&>;
 
 
        inline auto begin() {
            return iterator(*this, 0, 0);
        }
 
 
        inline auto end() {
            return iterator(*this, rows(), 0);
        }
 
 
        inline auto begin() const {
            return const_iterator(*this, 0, 0);
        }
 
 
        inline auto end() const {
            return const_iterator(*this, rows(), 0);
        }
 
 
        TH_CONSTEXPR inline unsigned int rows() const {
            return N;
        }
 
 
        TH_CONSTEXPR inline unsigned int cols() const {
            return K;
        }
 
 
        inline unsigned int size() const {
            return N * K;
        }
 
 
        template<typename Matrix>
        inline bool operator==(const Matrix& other) const {
            return algebra::mat_equals(*this, other);
        }
 
 
        template<typename Matrix>
        inline bool operator!=(const Matrix& other) const {
            return !algebra::mat_equals(*this, other);
        }
 
 
        inline bool is_square() const {
            return algebra::is_square(*this);
        }
 
 
        inline bool is_diagonal() const {
            return algebra::is_diagonal(*this);
        }
 
 
        inline bool is_symmetric() const {
            return algebra::is_symmetric(*this);
        }
 
 
        inline real density(real tolerance = 1E-12) {
            return algebra::density(*this, tolerance);
        }
 
 
        inline real sparsity(real tolerance = 1E-12) {
            return algebra::sparsity(*this, tolerance);
        }
 
 
        inline Type trace() {
            return algebra::trace(*this);
        }
 
 
        inline Type diagonal_product() {
            return algebra::diagonal_product(*this);
        }
 
 
        inline Type det() const {
            static_assert(N == K, "The matrix must be square to compute the determinant.");
            return algebra::det(*this);
        }
 
 
        inline mat<Type, N, K> inverse() const {
            static_assert(N == K, "The matrix must be square to be invertible.");
            return algebra::inverse(*this);
        }
 
 
        inline mat<Type, N, K>& invert() {
            static_assert(N == K, "The matrix must be square to be invertible.");
            return algebra::invert(*this);
        }
 
 
        inline mat<Type, N, K> resize(unsigned int n, unsigned int k) const {
 
            if(rows() != n) {
                TH_MATH_ERROR("mat::resize", n, INVALID_ARGUMENT);
            } else if(cols() != k) {
                TH_MATH_ERROR("mat::resize", k, INVALID_ARGUMENT);
            }
 
            return *this;
        }
 
 
        // Transformation matrices
 
 
        inline static mat<Type, N, K> identity() {
            return algebra::identity<mat<Type, N, K>>();
        }
 
 
        inline static mat<Type, N, K> diagonal(Type diag) {
            return mat<Type, N, K>(diag);
        }
 
 
        template<typename Vector = vec<real, N - 1>>
        inline static mat<Type, N, K> translation(Vector&& t) {
            return algebra::translation<mat<Type, N, K>>(t);
        }
 
 
        inline static mat<Type, N, K> rotation_2d(real theta) {
            static_assert(N >= 2 && K >= 2, "The matrix must be 2x2 or bigger");
            return algebra::rotation_2d<mat<Type, N, K>>(theta);
        }
 
 
        inline static mat<Type, N, K> rotation_3d_xaxis(real theta) {
            static_assert(N >= 3 && K >= 3, "The matrix must be 3x3 or bigger");
            return algebra::rotation_3d_xaxis<mat<Type, N, K>>(theta);
        }
 
 
        inline static mat<Type, N, K> rotation_3d_yaxis(real theta) {
            static_assert(N >= 3 && K >= 3, "The matrix must be 3x3 or bigger");
            return algebra::rotation_3d_yaxis<mat<Type, N, K>>(theta);
        }
 
 
        inline static mat<Type, N, K> rotation_3d_zaxis(real theta) {
            static_assert(N >= 3 && K >= 3, "The matrix must be 3x3 or bigger");
            return algebra::rotation_3d_zaxis<mat<Type, N, K>>(theta);
        }
 
 
        template<typename Vector = vec<real, 3>>
        inline static mat<Type, N, K> rotation_3d(real theta, Vector&& axis) {
            static_assert(N >= 3 && K >= 3, "The matrix must be 3x3 or bigger");
            return algebra::rotation_3d<mat<Type, N, K>>(theta, axis);
        }
 
 
        inline static mat<Type, N, K> perspective(
            real left, real right, real bottom,
            real top, real near, real far) {
 
            static_assert(N >= 4 && K >= 4, "The matrix must be 4x4 or bigger");
            return algebra::perspective<mat<Type, N, K>>(
                left, right, bottom, top, near, far);
        }
 
 
        inline static mat<Type, N, K> perspective_fov(
            real fov, real aspect, real near, real far) {
 
            static_assert(N >= 4 && K >= 4, "The matrix must be 4x4 or bigger");
            return algebra::perspective_fov<mat<Type, N, K>>(fov, aspect, near, far);
        }
 
 
        inline static mat<Type, N, K> ortho(
            real left, real right, real bottom, real top, real near, real far) {
            static_assert(N >= 4 && K >= 4, "The matrix must be 4x4 or bigger");
            return algebra::ortho<mat<Type, N, K>>(left, right, bottom, top, near, far);
        }
 
 
        template<typename Vector1, typename Vector2, typename Vector3>
        inline static mat<Type, 4, 4> look_at(
            const Vector1& camera, const Vector2& target, const Vector3& up) {
            return algebra::look_at<mat<Type, 4, 4>>(camera, target, up);
        }
 
 
        inline static mat<Type, N, K> symplectic(unsigned int n = 0, unsigned int k = 0) {
            static_assert(N == K && (N % 2 == 0),
                "N must equal K and they should be a multiple of 2");
            return algebra::symplectic<mat<Type, N, K>>(n, k);
        }
 
 
 
#ifndef THEORETICA_NO_PRINT
 
            inline std::string to_string(
                std::string separator = ", ", bool parenthesis = true) const {
 
                std::stringstream res;
 
                for (unsigned int i = 0; i < rows(); ++i) {
                        
                    if(parenthesis)
                        res << "(";
 
                    for (unsigned int j = 0; j < cols(); ++j) {
                        
                        if(j)
                            res << separator;
 
                        if(abs(get(i, j)) < MACH_EPSILON)
                            res << "0";
                        else
                            res << get(i, j);
                    }
        
                    if(parenthesis)
                        res << ")" << std::endl;
                }
 
                return res.str();
            }
 
 
            inline operator std::string() {
                return to_string();
            }
 
 
            inline friend std::ostream& operator<<(
                std::ostream& out, const mat<Type, N, K>& obj) {
                return out << obj.to_string();
            }
 
#endif
 
    };
 
 
 
    template<typename Type>
    class mat<Type, 0, 0> {
        public:
 
        // Container type for storage (alias for std::vector)
        template<typename T>
        using Container = std::vector<T>;
 
        Container<Container<Type>> data;
 
        size_t row_sz {0};
 
        size_t col_sz {0};
 
 
        mat() : row_sz(0), col_sz(0) {}
 
 
        mat(unsigned int n, unsigned int k) {
            resize(n, k);
            algebra::mat_zeroes(*this);
        }
 
 
        template <
            typename Matrix, enable_matrix<Matrix>
        >
        mat(const Matrix& m) {
            resize(m.rows(), m.cols());
            algebra::mat_copy(*this, m);
        }
 
 
        template<typename T = Type>
        inline mat(const std::initializer_list<std::initializer_list<T>>& rows) {
 
            unsigned int i = 0;
            unsigned int j = 0;
 
            for (const auto& row : rows) {
 
                for (const auto& x : row) {
                
                    if (!i && !j) {
                        this->resize(rows.size(), row.size());
                    }
                    else if (row.size() != col_sz) {
                        TH_MATH_ERROR("mat::mat", row.size(), INVALID_ARGUMENT);
                        algebra::mat_error(*this);
                        return;
                    }
 
                    at(i, j) = x;
                    j = (j + 1) % col_sz;
                }
 
                i++;
            }
        }
 
 
        template<typename Matrix>
        inline mat<Type>& operator=(const Matrix& other) {
            resize(other.rows(), other.cols());
            return algebra::mat_copy(*this, other);
        }
 
 
        mat(Type diagonal, unsigned int n, unsigned int k) {
            
            resize(n, k);
            algebra::mat_zeroes(*this);
            const unsigned int m = min(n, k);
 
            for (unsigned int i = 0; i < m; ++i)
                at(i, i) = diagonal;
        }
 
 
        ~mat() {
            row_sz = 0;
            col_sz = 0;
        }
 
 
        inline void make_zeroes() {
            algebra::mat_zeroes(*this);
        }
 
 
        inline static mat<Type> zeroes(unsigned int rows, unsigned int cols) {
            mat<Type> res;
            res.resize(rows, cols);
            algebra::mat_zeroes(res);
            return res;
        }
 
 
        template<typename Matrix>
        inline mat<Type> operator+(const Matrix& other) const {
            mat<Type> res;
            res.resize(rows(), cols());
            return algebra::mat_sum(res, *this, other);
        }
 
 
        template<typename Matrix>
        inline mat<Type> operator-(const Matrix& other) const {
            mat<Type> res;
            res.resize(rows(), cols());
            return algebra::mat_diff(res, *this, other);
        }
 
 
        inline mat<Type> operator*(Type scalar) const {
            mat<Type> res;
            res.resize(rows(), cols());
            return algebra::mat_scalmul(res, scalar, *this);
        }
 
 
        inline friend mat<Type> operator*(Type a, const mat<Type>& B) {
            return B * a;
        }
 
 
        template<typename VecType, unsigned int M>
        inline friend vec<VecType, 0> operator*(
            const vec<VecType, M>& a, const mat<Type, 0, 0>& B) {
            return B * a;
        }
 
 
        inline mat<Type> operator/(Type scalar) const {
 
            mat<Type> res;
            res.resize(rows(), cols());
 
            if(abs(scalar) < MACH_EPSILON) {
                TH_MATH_ERROR("mat::operator/", scalar, DIV_BY_ZERO);
                return algebra::mat_error(res);
            }
            
            return algebra::mat_scalmul(res, 1.0 / scalar, *this);
        }
 
 
        template<typename Vector>
        inline Vector transform(const Vector& v) const {
 
            if(v.size() != rows()) {
                TH_MATH_ERROR("mat::transform", v.size(), INVALID_ARGUMENT);
                Vector res;
                res.resize(rows());
                algebra::vec_error(res);
                return res;
            }
 
            return algebra::transform(*this, v);
        }
 
 
        template<unsigned int N = 0, unsigned int K = 0>
        inline vec<Type, N> transform(const vec<Type, K>& v) const {
            return algebra::transform(*this, v);
        }
 
 
        template<unsigned int N = 0, unsigned int K = 0>
        inline vec<Type, N> operator*(const vec<Type, K>& v) const {
            return transform(v);
        }
 
 
        // inline vec<Type> operator*(const vec<Type>& v) const {
        //  return transform(v);
        // }
 
 
        inline mat<Type> mul(const mat<Type>& B) const {
 
            mat<Type> res;
            res.resize(rows(), B.cols());
 
            if(B.rows() != cols()) {
                TH_MATH_ERROR("mat::mul", B.rows(), INVALID_ARGUMENT);
                algebra::mat_error(res);
                return res;
            }
 
            return algebra::mat_mul(res, *this, B);
        }
 
 
        template<typename Matrix>
        inline Matrix mul(const Matrix& B) const {
 
            Matrix res;
            res.resize(rows(), B.cols());
 
            if(B.rows() != cols()) {
                TH_MATH_ERROR("mat::mul", B.rows(), INVALID_ARGUMENT);
                algebra::mat_error(res);
                return res;
            }
 
            return algebra::mat_mul(res, *this, B);
        }
 
 
        template<typename Matrix>
        inline auto operator*(const Matrix& B) const {
            return mul(B);
        }
 
 
        template<typename Matrix>
        inline mat<Type>& operator+=(const Matrix& other) {
            return algebra::mat_sum(*this, other);
        }
 
 
        template<typename Matrix>
        inline mat<Type>& operator-=(const Matrix& other) {
            return algebra::mat_diff(*this, other);
        }
 
 
        inline mat<Type>& operator*=(Type scalar) {
            return algebra::mat_scalmul(scalar, *this);
        }
 
 
        inline mat<Type>& operator/=(Type scalar) {
 
            if(abs(scalar) < MACH_EPSILON) {
                TH_MATH_ERROR("mat::operator/", scalar, DIV_BY_ZERO);
                return algebra::mat_error(*this);
            }
 
            return algebra::mat_scalmul(1.0 / scalar, *this);
        }
 
 
        template<typename Matrix>
        inline mat<Type>& operator*=(const Matrix& B) {
            return (*this = this->operator*(B));
        }
 
 
        inline mat<Type>& transpose() {
            return algebra::make_transposed(*this);
        }
 
 
        inline mat<Type> transposed() const {
            return algebra::transpose<mat<Type>, mat<Type>>(*this);
        }
 
 
        inline Type& at(unsigned int i, unsigned int j) {
 
#ifdef THEORETICA_ROW_FIRST
            return data[i][j];
#else
            return data[j][i];
#endif
        }
 
 
        inline const Type& at(unsigned int i, unsigned int j) const {
 
#ifdef THEORETICA_ROW_FIRST
            return data[i][j];
#else
            return data[j][i];
#endif
        }
 
 
        inline Type& operator()(unsigned int i, unsigned int j) {
            return at(i, j);
        }
 
 
        inline const Type& operator()(unsigned int i, unsigned int j) const {
            return at(i, j);
        }
 
 
        inline Type get(unsigned int i, unsigned int j) const {
 
#ifdef THEORETICA_ROW_FIRST
            return data[i][j];
#else
            return data[j][i];
#endif
        }
 
 
        using iterator = mat_iterator<mat<Type, 0, 0>, Type&>;
 
 
        using const_iterator = mat_iterator<const mat<Type, 0, 0>, const Type&>;
 
 
        inline auto begin() {
            return iterator(*this, 0, 0);
        }
 
 
        inline auto end() {
            return iterator(*this, rows(), 0);
        }
 
 
        inline auto begin() const {
            return const_iterator(*this, 0, 0);
        }
 
 
        inline auto end() const {
            return const_iterator(*this, rows(), 0);
        }
 
 
        TH_CONSTEXPR inline unsigned int rows() const {
            return row_sz;
        }
 
 
        TH_CONSTEXPR inline unsigned int cols() const {
            return col_sz;
        }
 
 
        inline unsigned int size() const {
            return rows() * cols();
        }
 
 
        template<typename Matrix>
        inline bool operator==(const Matrix& other) const {
            return algebra::mat_equals(*this, other);
        }
 
 
        template<typename Matrix>
        inline bool operator!=(const Matrix& other) const {
            return !algebra::mat_equals(*this, other);
        }
 
 
        inline bool is_square() const {
            return algebra::is_square(*this);
        }
 
 
        inline bool is_diagonal() const {
            return algebra::is_diagonal(*this);
        }
 
 
        inline bool is_symmetric() const {
            return algebra::is_symmetric(*this);
        }
 
 
        inline real density(real tolerance = 1E-12) {
            return algebra::density(*this, tolerance);
        }
 
 
        inline real sparsity(real tolerance = 1E-12) {
            return algebra::sparsity(*this, tolerance);
        }
 
 
        inline Type trace() {
            return algebra::trace(*this);
        }
 
 
        inline Type diagonal_product() {
            return algebra::diagonal_product(*this);
        }
 
 
        inline Type det() const {
            return algebra::det(*this);
        }
 
 
        inline mat<Type> inverse() const {
            return algebra::inverse(*this);
        }
 
 
        inline mat<Type>& invert() {
            return algebra::invert(*this);
        }
 
 
        inline mat<Type>& resize(unsigned int n, unsigned int k) {
 
            // Do nothing if the size is already correct
            if(rows() == n && cols() == k)
                return *this;
 
            // Distinguish between row-first and column-first
            // allocation
#ifdef THEORETICA_ROW_FIRST
            size_t size1 = n, size2 = k;
#else
            size_t size1 = k, size2 = n;
#endif
 
            // The matrix must be allocated anew
            if(!data.size()) {
 
                data.resize(size1);
                for (unsigned int i = 0; i < size1; ++i)
                    data[i].resize(size2);
                
                row_sz = n;
                col_sz = k;
                algebra::mat_zeroes(*this);
 
            // The matrix must be reallocated
            } else {
 
                std::vector<std::vector<Type>> new_data;
                new_data.resize(size1);
                for (unsigned int i = 0; i < size1; ++i)
                    new_data[i].resize(size2);
 
                // Copy data to new memory location
 
                // Bounds on the region to copy
                const size_t row_bound = min(rows(), n);
                const size_t col_bound = min(cols(), n);
 
                for (unsigned int i = 0; i < row_bound; ++i) {
                    for (unsigned int j = 0; j < col_bound; ++j) {
#ifdef THEORETICA_ROW_FIRST
                        new_data[i][j] = get(i, j);
#else
                        new_data[j][i] = get(i, j);
#endif
                    }
                }
 
                // If the new matrix size is bigger,
                // zero out all remaining entries
                for (unsigned int i = 0; i < n; ++i) {
                    for (unsigned int j = col_bound; j < k; ++j) {
 
 
#ifdef THEORETICA_ROW_FIRST
                        new_data[i][j] = (Type) 0;
#else
                        new_data[j][i] = (Type) 0;
#endif
                    }
                }
 
 
                // If the new matrix size is bigger,
                // zero out all remaining entries
                for (unsigned int i = row_bound; i < n; ++i) {
                    for (unsigned int j = 0; j < k; ++j) {
 
 
#ifdef THEORETICA_ROW_FIRST
                        new_data[i][j] = (Type) 0;
#else
                        new_data[j][i] = (Type) 0;
#endif
                    }
                }
 
                data.clear();
                data = new_data;
 
                row_sz = n;
                col_sz = k;
            }
 
            return *this;
        }
 
 
        // Transformation matrices
 
 
        inline static mat<Type> identity(
            unsigned int row_sz, unsigned int col_sz) {
 
            return algebra::identity<mat<Type>>(row_sz, col_sz);
        }
 
 
        inline static mat<Type> diagonal(
            Type diag, unsigned int row_sz, unsigned int col_sz) {
            return mat<Type>(diag, row_sz, col_sz);
        }
 
 
        template<typename Vector>
        inline static mat<Type> translation(Vector&& t) {
            return algebra::translation<mat<Type>>(t);
        }
 
 
        inline static mat<Type> rotation_2d(real theta) {
            return algebra::rotation_2d<mat<Type>>(theta);
        }
 
 
        inline static mat<Type> rotation_3d_xaxis(real theta) {
            return algebra::rotation_3d_xaxis<mat<Type>>(theta);
        }
 
 
        inline static mat<Type> rotation_3d_yaxis(real theta) {
            return algebra::rotation_3d_yaxis<mat<Type>>(theta);
        }
 
 
        inline static mat<Type> rotation_3d_zaxis(real theta) {
            return algebra::rotation_3d_zaxis<mat<Type>>(theta);
        }
 
 
        template<typename Vector = vec<real, 3>>
        inline static mat<Type> rotation_3d(real theta, Vector&& axis) {
            return algebra::rotation_3d<mat<Type>>(theta, axis);
        }
 
 
        inline static mat<Type> perspective(
            real left, real right, real bottom,
            real top, real near, real far) {
 
            return algebra::perspective<mat<Type>>(
                left, right, bottom, top, near, far);
        }
 
 
        inline static mat<Type> perspective_fov(
            real fov, real aspect, real near, real far) {
 
            return algebra::perspective_fov<mat<Type>>(fov, aspect, near, far);
        }
 
 
        inline static mat<Type> ortho(
            real left, real right, real bottom, real top, real near, real far) {
            return algebra::ortho<mat<Type>>(left, right, bottom, top, near, far);
        }
 
 
        template<typename Vector1, typename Vector2, typename Vector3>
        inline static mat<Type> look_at(
            const Vector1& camera, const Vector2& target, const Vector3& up) {
            return algebra::look_at<mat<Type>>(camera, target, up);
        }
 
 
        inline static mat<Type> symplectic(unsigned int rows, unsigned int cols) {
            return algebra::symplectic<mat<Type>>(rows, cols);
        }
 
 
 
 
#ifndef THEORETICA_NO_PRINT
 
            inline std::string to_string(
                std::string separator = ", ", bool parenthesis = true) const {
 
                std::stringstream res;
 
                for (unsigned int i = 0; i < rows(); ++i) {
                        
                    if(parenthesis)
                        res << "(";
 
                    for (unsigned int j = 0; j < cols(); ++j) {
                        
                        if(j)
                            res << separator;
 
                        if(abs(get(i, j)) < MACH_EPSILON)
                            res << "0";
                        else
                            res << get(i, j);
                    }
        
                    if(parenthesis)
                        res << ")" << std::endl;
                }
 
                return res.str();
            }
 
 
            inline operator std::string() {
                return to_string();
            }
 
 
            inline friend std::ostream& operator<<(
                std::ostream& out, const mat<Type>& obj) {
                return out << obj.to_string();
            }
 
#endif
 
    };
 
}
 
#endif