Template Matrix Class: implemented some basic functionalities

Question

Basically I have written a matrix class for addition, multiplication and scalar multiplication.

I need your review of the class implementation below in terms of efficiency, memory consumption and new C++11/14/17 features that can be used.

Here goes the code:

Matrix.hpp

// Created by prajwal.sapare on 2/18/2020.
//

#ifndef UNTITLED2_MATRIX_HPP
#define UNTITLED2_MATRIX_HPP

#include <vector>
#include <exception>
#include <iostream>
#include <algorithm>
#include "helper.h"

template<class T, uint m, uint n>
class Matrix;

template <class T, uint n>
using rowVector = Matrix<T, 1, n>;

template <class T, uint m>
using colVector = Matrix<T, m, 1>;

template<class T, uint m, uint n>
class Matrix {
public:
    explicit Matrix();
    explicit Matrix(const std::vector<T> matrixValue);

    int getRows()const;
    int getColoumns() const;

    T& operator()(const uint row, const uint col);

    Matrix<T,m,n> operator=(std::vector<T> matrixIntializationValue);

    template<uint N>
    Matrix<T,m,N> operator*(Matrix<T,n,N>& other);

    template<uint a, uint b>
    bool operator==(Matrix<T,a,b>& other);

    Matrix<T,m,n> operator+(Matrix<T,m,n>& other);
    Matrix<T,m,n> operator*(T scalar);

    template<class T, uint m, uint n> friend std::ostream& operator<<(std::ostream& os, const Matrix<T,m,n>& matrix);

private:
    uint rows;
    uint cols;
    std::vector<std::vector<T>> data;
};

template<class T, uint m, uint n>
Matrix<T,m,n>::Matrix(): rows(m), cols(n)
{
    if(rows == 0 || cols == 0)
        throw std::invalid_argument( "received zero as argument" );

    data.resize(rows);
    for (auto& colData : data)
        colData.resize(cols);
}

template<class T, uint m, uint n>
Matrix<T,m,n>::Matrix(const std::vector<T> matrixValue): rows(m), cols(n)
{
    if(rows == 0 || cols == 0)
        throw std::invalid_argument( "received zero as argument" );

    if(matrixValue.empty())
        throw std::invalid_argument( "Empty vector" );

    if(rows * cols != matrixValue.size())
        throw std::runtime_error( "Total number of matrix values does not match with rows and coloumns provided" );

    data.resize(rows);
    for (auto& colData : data)
        colData.resize(cols);

    for (auto i = 0; i<rows; i++)
        data[i] = {matrixValue.begin() + (i * cols), matrixValue.begin() + ((i+1) * cols)};
}

template<class T, uint m, uint n>
int Matrix<T,m,n>::getRows() const
{
    return this->rows;
}

template<class T, uint m, uint n>
int Matrix<T,m,n>::getColoumns() const
{
    return this->cols;
}

template<class T, uint m, uint n>
T& Matrix<T,m,n>::operator()(const uint row, const uint col)
{
    return this->data[row][col];
}

template<class T, uint m, uint n>
Matrix<T,m,n> Matrix<T,m,n>::operator+(Matrix& other)
{
    if ((this->rows == other.getRows()) && (this->cols == other.getColoumns()))
    {
        Matrix<T,m,n> resultantMatrix;

        for(auto i = 0; i< this->rows; i++)
        {
            for(auto j = 0; j < this->cols; j++)
            {
                auto& valueFirst = this->data[i][j];
                auto& valueSecond = other(i, j);

                if((additionOverflow(valueFirst, valueSecond)) || (additionUnderflow(valueFirst, valueSecond)))
                    throw std::out_of_range("Resultant value of matrix is out of range");
                else
                    resultantMatrix(i, j) = valueFirst + valueSecond;
            }
        }
        return resultantMatrix;
    }
    else
        throw std::runtime_error("Matrices cannot be added, sizes do not match");
}

template<class T, uint m, uint n>
template<uint N>
Matrix<T,m,N> Matrix<T,m,n>::operator*(Matrix<T,n,N>& other)
{
    if ((this->cols == other.getRows()))
    {
        Matrix<T,m,N> resultantMatrix;
        T temp = 0;

        for (auto i = 0; i < this->rows; i++) {
            for (auto j = 0; j < other.getColoumns(); j++) {
                temp = 0.0;
                for (auto k = 0; k < this->cols; k++) {
                    if(other(k, j) != 0)
                    {
                       if(multiplicationOverflow(this->data[i][k], other(k, j)) || multiplicationUnderflow(this->data[i][k], other(k, j)))
                            throw std::out_of_range("Resultant value of matrix is out of range");
                    }
                    auto tempMul = this->data[i][k] * other(k, j);

                    if((additionOverflow(temp, tempMul)) || (additionUnderflow(temp, tempMul)))
                        throw std::out_of_range("Resultant value of matrix is out of range");

                    temp = temp + tempMul;
                }
                resultantMatrix(i, j) = temp;
            }
        }
        return resultantMatrix;
    }
    else
        throw std::runtime_error("Matrices cannot be multiplied, invalid sizes");
}

template<class T, uint m, uint n>
Matrix<T,m,n> Matrix<T,m,n>::operator*(T scalar){
    Matrix<T,m,n> resultantMatrix;

    for (auto i = 0; i < this->rows; i++)
    {
        for (auto j = 0; j < this->cols; j++)
        {
            auto valueFirst = this->data[i][j];
            if (scalar == 0)
                resultantMatrix(i,j) = 0;
            else if((multiplicationOverflow(valueFirst, scalar)) || (multiplicationUnderflow(valueFirst, scalar)))
                throw std::out_of_range("Resultant value of matrix is out of range");
            else
                resultantMatrix(i,j) = this->data[i][j] * scalar;
        }
    }
    return resultantMatrix;
}

template<class T, uint m, uint n>
Matrix<T,m,n> Matrix<T,m,n>::operator=(std::vector<T> matrixIntializationValue)
{
    if (!matrixIntializationValue.empty())
    {
        if(rows * cols != matrixIntializationValue.size())
            throw std::invalid_argument( "Total number of matrix values does not match with rows and coloumns provided" );

        for (auto i = 0; i<rows; i++)
            data[i] = {matrixIntializationValue.begin() + (i * cols), matrixIntializationValue.begin() + ((i+1) * cols)};

        return *this;
    }
    else
        throw std::runtime_error("Matrices cannot be multiplied, invalid sizes");
}

template<class T, uint m, uint n>
std::ostream& operator<<(std::ostream &os, const Matrix<T,m,n> &matrix)
{
    if(!matrix.data.empty())
    {
        for (const auto& rowVal : matrix.data)
        {
            for(const auto& colVal : rowVal)
                os << colVal << " ";
            os << std::endl;
        }
    }
    return os;
}

template<class T, uint m, uint n>
template<uint a, uint b>
bool Matrix<T,m,n>::operator==(Matrix<T,a,b>& other)
{
    if( (m!=a) || (n!=b))
        return false;

    for(auto i =0; i< m ; i++)
    {
        for(auto j=0;j <n;j++)
        {
            if(this->data[i][j]!= other(i,j))
                return false;
        }
    }

    return true;
}

#endif //UNTITLED2_MATRIX_HPP

Helper.h

//
// Created by prajwal.sapare on 2/22/2020.
//

#ifndef UNTITLED2_HELPER_H
#define UNTITLED2_HELPER_H

using uint = unsigned int;

#include <limits>

template <typename T>
constexpr bool additionOverflow(const T& x, const T& y) {
    return (y >= 0) && (x > std::numeric_limits<T>::max() - y);
}

template <typename T>
constexpr bool additionUnderflow(const T& x, const T& y) {
    return (y < 0) && (x < std::numeric_limits<T>::min() - y);
}

template <typename T>
constexpr bool multiplicationOverflow(const T& x, const T& y) {
    return ((y >= 0) && (x >= 0) && (x > std::numeric_limits<T>::max() / y))
           || ((y < 0) && (x < 0) && (x < std::numeric_limits<T>::max() / y));
}

template <typename T>
constexpr bool multiplicationUnderflow(const T& x, const T& y) {
    return ((y >= 0) && (x < 0) && (x < std::numeric_limits<T>::min() / y))
           || ((y < 0) && (x >= 0) && (x > std::numeric_limits<T>::min() / y));
}

#endif //UNTITLED2_HELPER_H

Answer 1

Overview

Potentially not very efficient.

I have seen other Matrix libraries not do the actual computation until the value is needed. So what is passed around is a Matrix wrapper that contains references to the original matrixes and the operations on them. You can then delay the actual work till it is needed. Then if you multiply by the null matrix you can can optimize out all the numeric operations.

I will assume you know how to do the basic maths of matrix manipulation. I will look at the interfaces and idioms.

Code Review

Nice:

template <class T, uint n>
using rowVector = Matrix<T, 1, n>;

template <class T, uint m>
using colVector = Matrix<T, m, 1>;

---

So the only way to initialize the matrix is via a vector?

    explicit Matrix(const std::vector<T> matrixValue);
                    ^^^^^^^^^^^^^^^^^^^^
                    Note: passing by value.
                          You are making a copy. Probably want to pass
                          by const reference (see below)             

    explicit Matrix(std::vector<T> const& matrixValue);
    explicit Matrix(std::vector<T>&&      matrixValue); // don't forget move semantics.

I would have liked to be able to initialize it via iterators

    template<typename I>
    explicit Matrix(I begin, I end);

That way I can simply read the data from a file:

    std::ifstream    file("Matrix.file");
    Matrix<int, 5, 6>  data(std::istream_iterator<int>{file}, std::istream_iterator<int>{});

---

I like this:

    int getRows()const;
    int getColoumns() const;

But some people will complain that you should use unsigned int (i.e. std::size_t).

---

Yes you need this:

    T& operator()(const uint row, const uint col);

But you also need the const version. A lot of time you may pass a const reference of the matrix to a function. You should still be able to read data from the matrix.

    T const& operator()(const uint row, const uint col) const;

If you want to fo whole hog you can add the [] operators. A tiny bit more work but doable.

    // Define inside Matrix so the template stuff
    // Is already available.
    struct Row
    {
        Matrix&   parent;
        int       row;
        Row(Matrix& parent, int row) : parent(parent), row(row) {}
        T&       operator[](int col)       {return parent(row, col);}
    }
    Row operator[](int row){return Row(*this, row);}

    // Very similar for the const version.

---

Fair:

    Matrix<T,m,n> operator=(std::vector<T> matrixIntializationValue);

But if you define this. Then you also need to define the standard assignment operator. Otherwise you can't do simple Matrix assignment.

---

You don't need all this template stuff if you define this inside the class declaration.

    template<class T, uint m, uint n> friend std::ostream& operator<<(std::ostream& os, const Matrix<T,m,n>& matrix);

By doing it inside the class you make everything much simpler to read.

---

Why are you storing the data as a vector of vectors?

std::vector<std::vector<T>> data;

You have this wrapper class for managing all the interactions. You can use a more effecient storage method internally because you have wrapped the object and limit access to the implementation.

Simply use a single vector of data.

std::vector<T> data;

It makes initialization easy. Access will be quicker (less look-ups and data locality will mean that data will already be in the processor cache).

---

Since m and n are defined at compile time. This runtime check should not be required. You simply make the zero sized versions invalid and you get a compile time error.

template<class T, uint m, uint n>
Matrix<T,m,n>::Matrix(): rows(m), cols(n)
{
    // Compile time check.
    static_assert(m > 0 && n > 0);

---

No need for this loop:

    data.resize(rows);
    for (auto& colData : data)
        colData.resize(cols);

    ----

Simply do the declaration in the initializer list:

template<class T, uint m, uint n>
Matrix<T,m,n>::Matrix()
     : rows(m)
     , cols(n)
     , data(rows, std::vector<T>(cols))
{}

---

Stop using: this->

    return this->rows;

Its not very C++. It also causes errors. Not for the compiler but for developers. The only reason you need this-> is when you have shadowed a variable and need to distinguish between the local and the object variable. Unfortunately when you miss this-> off the compiler can't tell you have made a mistake.

But by using distinct names you know which variable you are using and the compiler can warn you about shadowed variables (turn on the warning and make it an error). This will reduce the number of programmer mistakes in your code.

Well named and distinct variables will produce better code and less mistakes.

---

You don't need this check:

template<class T, uint m, uint n>
Matrix<T,m,n> Matrix<T,m,n>::operator+(Matrix& other)
{

Did you not make sure that this addition can only happen with very specific types in the class definition? This means that Matrixes of the wrong size would generate a compile time error and thus you don't need a runtime error.

    if ((this->rows == other.getRows()) && (this->cols == other.getColoumns()))
    {

---

template<class T, uint m, uint n>
std::ostream& operator<<(std::ostream &os, const Matrix<T,m,n> &matrix)
{

Can the data ever be empty?

    if(!matrix.data.empty())

Don't you guarantee in the constructor that this never happens?

---

Don't use std::endl when "\n" will do.

            os << std::endl;

The difference is that endl will force a stream flush. Manually flushing the stream is almost always (I say almost because somebody will point out an obscure situation were it would be nice) the wrong thing to do. You should never be using it until you can show with testing that it is worth the flush.

When people complain about the speed of the C++ streams this is one of the main culprits. Excessive flushing of the stream is one of the main slowdown people see. When you stop doing it the C++ streams are about the same speed as C streams.

Answer 2

Correctness:

You are dividing by 0 here

((y >= 0) && (x >= 0) && (x > std::numeric_limits<T>::max() / y)

and here

((y >= 0) && (x < 0) && (x < std::numeric_limits<T>::min() / y))

Note that division can also underflow and overflow.

---

Efficiency:

Your code is about 1000 times slower than MKL for multiplying large float matrices, which is state-of-the-art, but there are other implementations that have similar performance to MKL, such as OpenBLAS.

The reasons your code is hugely inefficient include

1. The layout of your data requires extra indirection (one chunk of data is better)

2. The loops are cache-unfriendly (The triple loop is inefficient for dot products)

3. Division is much more expensive than multiplication, typically, and you are doing it twice to check that multiplying is "OK" (which this fails to do)

---

Some minor things I noticed:

constexpr seems pointless here, since you don't know the arguments at compile time.

this-> is unnecessary

Matrix dimensions are uint, but the methods returning them return int

There is no need to store rows and cols in the class, since they are known to the compiler (template parameters)

---

Odd choice:

Your class can't handle arbitrary matrix sizes, unless they are known at compile time.

But since you are doing that, you also should avoid throwing a runtime error, if matrix sizes do not match -- you should know this before the program runs.

---

Compilation errors:

Shadowing template parameters:

template<class T, uint m, uint n>
class Matrix {
// ...
template<class T, uint m, uint n> friend std::ostream& operator<<(std::ostream& os, const Matrix<T,m,n>& matrix);
// ...
};

Answer 3

I am not going to be nice. That's how bad it is. The design is terrible. You are greatly confused on so many levels.

Internally, matrix stores runtime rows, cols, and dynamically the data in std::vector<std:vector> - while the rows and cols are defined in the template arguments and thus compile time constants. It makes no sense at all.

Generally there are two types of matrix implementations:

(1) struct with fixed rows and cols and its data is stored in, say, double mData[rows][cols]; - a static format.

(2) dynamically allocated matrix with run-time chosen rows and cols with dynamically allocated data in, say, std::vector<double> data;.

(Note: eigen linear algebra library also has more intermediate types for greater compatiability and efficiency.)

Method (1) has the advantage that everything is known at compile time and data is allocated statically. Con: It requires compile time knowledge of the size of matrices and has to compile a new class for each size. This is the preferred method for small sized matrices.

Method (2) can serve all matrices of all sizes but it requires dynamic allocation of data and lack of compile time knowledge of sizes might infer with optimizations and some bugs will not be discoverable at compile time. This is the preferred method for large sized matrices and those when one lacks knowledge of matrix size at compile time.

Your version has all the cons and none of the advantages. It requires compile time knowledge of the matrix's size and it requires run-time allocations and for some mysterious reasons it stores the size (compile-time defined) inside the matrix instance. Furthermore, you store data inside std::vector<std::vector<double>> which is a horrible idea in as it requires rows+1 allocations - it is slow on both initialization and usage and in addition results in memory fragmentation problems. This is how bad it is.

Please, just use an open source matrix library. There are plenty.