Simple matrix library in C++

Question

Today I started learning C++ and at the end of the day have made a simple generic matrix class. I'm looking for feedback on techniques and features. It's still not complete, though. But everything works!

I put some functions outside the class and inside a namespace because I was trouble defining template member functions with partial specialization.

Matrix.hpp

#pragma once

#include <iostream>

template <typename T, int m, int n = m>
class mat {

template <typename, int, int>
friend class mat;

private:

    T * data;

public:

    mat(const std::initializer_list<T> & ini) {
        std::copy(ini.begin(), ini.end(), data);
    }

    mat() : data(new T[m * n]) {
    }

    mat(T * values) : data(values) {

    }

    mat(const mat & mat2) : mat(){
        for (int i = 0; i < m; ++i)
            for (int j = 0; j < n; ++j)
                data[i * n + j] = mat2(i, j);
    }

    ~mat() {
        delete data;
    }

    void fill(T val) {
        for (int i = 0; i < m; ++i)
            for (int j = 0; j < n; ++j)
                data[i * n + j] = val;
    }

    void print() {
        for (int i = 0; i < m; ++i) {
            for (int j = 0; j < n; ++j)
                std::cout << data[i * n + j] << " ";
            std::cout << std::endl;
        }
    }

    mat<T, m, 1> col(int j) const {
        mat<T, m, 1> col;
        for (int i = 0; i < m; ++i)
            col.data[i] = data[i * n + j];
        return col;
    }

    mat<T, 1, n> row(int i) const {
        mat<T, 1, n> row;
        std::copy_n(data + i * n, n, row.data);
        return row;
    }

    T operator () (int row, int col) const {
        return data[row * n + col];
    }

    T & operator () (int row, int col) {
        return data[row * n + col];
    }

    mat<T, n, m> transpose() const {
        mat<T, n, m> res;
        for (int i = 0; i < m; ++i)
            for (int j = 0; j < n; ++j)
                res(j, i) = this->data[i * n + j];
        return res;
    }

    template <int p>
    mat<T, m, p> operator * (const mat<T, n, p> & other) const {
        mat<T, m, p> res;
        for (int i = 0; i < m; ++i)
            for (int j = 0; j < p; ++j) {
                T sum = 0;
                for (int k = 0; k < n; ++k)
                    sum += this->data[i * n + k] * other.data[k * n + j];
                res(i, j) = sum;
            }
        return res;
    }

    mat<T, m, n> & operator *= (T val) {
        for (int i = 0; i < m; ++i)
            for (int j = 0; j < n; ++j)
                data[i * n + j] *= val;
        return *this;
    }

    mat<T, m, n> & operator /= (T val) {
        return this *= 1 / val;
    }

    mat<T, m, n> & operator += (const mat & other) {
        for (int i = 0; i < m; ++i)
            for (int j = 0; j < n; ++j)
                data[i * n + j] += other(i, j);
        return *this;
    }

    mat<T, m, n> & operator -= (const mat & other) {
        for (int i = 0; i < m; ++i)
            for (int j = 0; j < n; ++j)
                data[i * n + j] -= other(i, j);
        return *this;
    }

    mat<T, m, n> operator * (T val) const {
        auto res(*this);
        res *= val;
        return res;
    }

    mat<T, m, n> operator / (T val) const {
        return (*this) * (1 / val);
    }

    mat<T, m, n> operator + (const mat & other) const {
        auto res(*this);
        res += other;
        return res;
    }

    mat<T, m, n> operator - (const mat & other) const {
        auto res(*this);
        res -= other;
        return res;
    }

    mat<T, m - 1, n - 1> cut(int row, int col) const {
        mat<T, m - 1, n - 1> res;
        int index = 0;

        for (int i = 0; i < m; ++i)
            for (int j = 0; j < n; ++j) {
                if (i == row || j == col)
                    continue;
                res.data[index++] = data[i * n + j];
            }

        return res;
    }

};

namespace matrix {

    template <typename T>
    T det(const mat<T, 1, 1> & arg) {
        return arg(0, 0);
    }

    template <typename T>
    T det(const mat<T, 2, 2> & arg) {
        return arg(0, 0) * arg(1, 1) - arg(1, 0) * arg(0, 1);
    }

    template <typename T, int n>
    T det(const mat<T, n, n> & arg) {
        T res = 0, coef = 1;
        for (int i = 0; i < n; ++i, coef *= -1) {
            res += coef * arg(0, i) * matrix::det(arg.cut(0, i));
        }
        return res;
    }

    template <typename T>
    mat<T, 2, 2> inv(const mat<T, 2, 2> & arg) {
        mat<T, 2, 2> helper;
        helper(1, 1) = arg(0, 0);
        helper(0, 0) = arg(1, 1);
        helper(0, 1) = -arg(0, 1);
        helper(1, 0) = -arg(1, 0);
        return helper / det(arg);
    }

    template <typename T, int m>
    mat<T, m, m> id() {
        mat<T, m, m> res;
        for (int i = 0; i < m; ++i)
            res(i, i) = 1;
        return res;
    }

};

Main.cpp

#include <iostream>
#include <exception>

#include "Matrix.hpp"

using namespace std;


int main() {

    constexpr int size = 4;

    long * data = new long[size * size]{ 0 };
    for (int i = 0; i < size; ++i)
        data[size * i + i] = 2;

    mat<long, size> big(data);
    cout << matrix::det(big) << endl;

    return 0;
}

Answer 1

mat (const std::initializer_list & ini)

You forgot to allocate any memory.

Rule of 0/5

Your class should follow the rule of 5 which states you should define copy-constructor, copy-assignment operator, move-constructor, move-assignment operator and destructor if you define one of them. You are missing copy-assignment operator, move-constructor, and move-assignment operator. You code cannot be efficiently moved. Implicitly defined copy-assignment operator is incorrect because you are dealing with memory ownership and the default copy-assignment operator can leads to double freeing the memory.

You could also follow the rule of zero and delegate the responsibility of memory management to something like std::unique_ptr.

Since the size of your matrix are template parameters you could just allocate the data as part of your class rather then dynamically allocating the data.

T data[n][m];

mat (T * values) : data (values)

In this constructor you are stealing the ownership of the allocated memory. This is bad style because you could accidentally free the memory from outside of the class or pass the same memory to multiple mat object, each will free the same memory.

You should allocate new memory and copy the values over.

mat (const mat & mat2)

Could call fixed version of mat (T * values) instead to reduce code duplication. Could also be done using std::copy_n

~mat ()

Should be delete[] data; since you are using new[].

void fill (T val)

Can be done using std::fill_n

Row and column iterator

You could provide row and column iterators and implement some of your functions using those iterator. That would make your code easier to understand and easier to spot error.

operator *=, operator * (T val)

Can be done using std::for_each_n

operator += (const mat & other), operator -= (const mat & other)

Can be done using std::transform

operator += (T value)

You are missing the version that add constant to matrix. Also what about matrix multiplication?

print()

Print should take in std::ostream & that way you can print to files as well as std::cout.

Answer 2

The C++ side seems covered, so in this answer I take more of a linear algebra angle.

The algorithms used for the determinant is not practical for matrixes that are not small. Cofactor expansion works fine with small matrixes and it is a useful mathematical definition, but its computational expense scales as the factorial of the matrix size, which quickly grows to unreasonable levels.

There are some other options. You could implement triangular factorization (eg LUP), after which the determinant can be found by multiplying the diagonal elements. The cost of this approach scales only as the cube of the matrix size, though it involves divisions which makes it not so efficient for small matrixes where such overhead still dominates. The divisions also mean that it only works with T being float or double, there are some techniques that do work for integers but they are probably out of scope for a beginner question. To counter these disadvantages you could keep cofactor expansion as well, and apply the best algorithm for the case depending on T and the matrix dimensions (this seems like a nice place to apply some template magic).

Triangular factorization can also be used to efficiently solve equations of the form Ax=b, and therefore also to find the explicit inverse of A by solving for all the basis vectors, though most commonly this wouldn't be needed since the purpose of the inverse would typically be to solve such equations in the first place.