Barebones C++ Matrix class

Question

I have implemented a C++ Matrix class using std::vector and a number for rows/cols. The implementation works decently from the QA I've done. I have implemented a Vector class as a derived class but I'm unsure whether this implementation is optimal. Would it make more sense to make Vector the base class and Matrix the derived class with Vector being just m_data and Matrix adding m_nrow and m_ncol? Any general feedback more than welcome.

matrix.h:

#pragma once
#include <vector>

class Vector;

class Matrix
{
protected:
    size_t m_nrow{};
    size_t m_ncol{};
    std::vector<double> m_data{};

public:
    Matrix(const size_t nrow = 0, const size_t ncol = 0, const std::vector<double>& data = std::vector<double>{});
    void clear();
    void setElement(const double value, const size_t row_idx, const size_t col_idx);
    void setElement(const double value, const size_t idx);
    void removeRow(const size_t row_idx);
    void removeCol(const size_t col_idx);
    double& at(const size_t row_idx, const size_t col_idx);
    double& at(const size_t idx);
    double getElement(const size_t row_idx, const size_t col_idx) const;
    double getElement(const size_t idx) const;
    const std::vector<double>& getDataAsVector() const;
    size_t nRow() const;
    size_t nCol() const;
    void print() const;
    bool isSquare() const;
    void operator+=(const Matrix& other_matrix);
    void operator-=(const Matrix& other_matrix);
    void operator*=(const double multiplier);
    const Matrix operator+(const Matrix& other_matrix) const;
    const Matrix operator-(const Matrix& other_matrix) const;
    const Matrix operator*(const double multiplier) const;
    const Matrix operator*(const Matrix& other_matrix) const;
    const Vector operator*(const Vector& other_vector) const;
    double operator[](const size_t idx) const;
};

class Vector : public Matrix
{
private:
    using Matrix::isSquare; //hide this from the Vector class

public:
    Vector(const std::vector<double>& data = std::vector<double>{});
    double length() const;
    double dotProduct(const Vector& other_vector) const;
    double distanceTo(const Vector& other_vector) const;
    const Vector operator+(const Vector& other_vector) const;
    const Vector operator-(const Vector& other_vector) const;
    const Vector operator*(const double multiplier) const;
};

matrix.cpp:

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

Matrix::Matrix(const size_t nrow, const size_t ncol, const std::vector<double>& data) //constructor, called when an object is created, don't include default vars here
    : m_nrow{ nrow }, m_ncol{ ncol }, m_data{ data }
{
    if ((nrow * ncol) != data.size())
        throw std::invalid_argument("Data size does not match the dimensions!");
    if ((nrow == 0 && ncol > 0) || (nrow > 0 && ncol == 0))
        throw std::invalid_argument("One dimension is zero while the other one is not!");
}

void Matrix::clear()
{
    m_nrow = 0;
    m_ncol = 0;
    m_data.clear();
}

void Matrix::setElement(const double value, const size_t row_idx, const size_t col_idx)
{
    if (row_idx >= m_nrow || col_idx >= m_ncol)
        throw std::invalid_argument("Index exceeds dimensions!");
    m_data.at(row_idx * m_ncol + col_idx) = value;
}

void Matrix::setElement(const double value, const size_t idx)
{
    if (idx >= m_nrow * m_ncol)
        throw std::invalid_argument("Index exceeds dimensions!");
    m_data.at(idx) = value;
}

void Matrix::removeRow(const size_t row_idx)
{
    if (row_idx >= m_nrow)
        throw std::invalid_argument("Index exceeds dimensions!");
    else if (m_nrow == 1)
    {
        m_nrow = 0;
        m_ncol = 0;
        m_data.clear();
    }
    else
    {
        size_t idx_start{ row_idx * m_ncol };
        m_data.erase(std::next(m_data.begin(), idx_start), std::next(m_data.begin(), idx_start + m_ncol));
        m_nrow--;

        if ((m_nrow * m_ncol) != m_data.size())
            throw std::invalid_argument("Data size does not match the dimensions!");
    }
}

void Matrix::removeCol(const size_t col_idx)
{
    if (col_idx >= m_ncol)
        throw std::invalid_argument("Index exceeds dimensions!");
    else if (m_ncol == 1)
    {
        m_nrow = 0;
        m_ncol = 0;
        m_data.clear();
    }
    else
    {
        size_t idx_from_back{ (m_nrow - 1) * m_ncol + col_idx };

        for (size_t i{ 0 }; i < m_nrow; i++)
        {
            m_data.erase(m_data.begin() + idx_from_back);
            idx_from_back -= m_ncol;
        }
        m_ncol--;

        if ((m_nrow * m_ncol) != m_data.size())
            throw std::invalid_argument("Data size does not match the dimensions!");
    }
}

double& Matrix::at(const size_t row_idx, const size_t col_idx)
{
    if (row_idx >= m_nrow || col_idx >= m_ncol)
        throw std::invalid_argument("Index exceeds dimensions!");
    return m_data.at(row_idx * m_ncol + col_idx);
}

double& Matrix::at(const size_t idx)
{
    if (idx >= m_nrow * m_ncol)
        throw std::invalid_argument("Index exceeds dimensions!");
    return m_data.at(idx);
}

double Matrix::getElement(const size_t row_idx, const size_t col_idx) const
{
    if (row_idx >= m_nrow || col_idx >= m_ncol)
        throw std::invalid_argument("Index exceeds dimensions!");
    return m_data.at(row_idx * m_ncol + col_idx);
}

double Matrix::getElement(const size_t idx) const
{
    if (idx >= m_nrow * m_ncol)
        throw std::invalid_argument("Index exceeds dimensions!");
    return m_data.at(idx);
}

const std::vector<double>& Matrix::getDataAsVector() const
{
    return m_data;
}

size_t Matrix::nRow() const
{
    return m_nrow;
}

size_t Matrix::nCol() const
{
    return m_ncol;
}

void Matrix::print() const
{
    std::cout << "Nrow: " << m_nrow << " Ncol: " << m_ncol << '\n';
    std::cout << "Size: " << m_data.size() << " Capacity: " << m_data.capacity() << '\n';
    size_t idx{ 0 };
    for (size_t i{ 0 }; i < m_nrow; i++)
    {
        for (size_t j{ 0 }; j < m_ncol; j++)
        {
            std::cout << m_data.at(idx) << ' ';
            idx++;
        }
        std::cout << '\n';
    }
    std::cout << '\n';
}

bool Matrix::isSquare() const
{
    return (m_nrow == m_ncol);
}

void Matrix::operator+=(const Matrix& other_matrix)
{
    if (m_ncol != other_matrix.m_ncol || m_nrow != other_matrix.m_nrow)
        throw std::invalid_argument("Matrices have different dimensions!");

    for (size_t i{ 0 }; i < m_data.size(); i++)
    {
        m_data.at(i) += other_matrix.m_data.at(i);
    }
}

void Matrix::operator-=(const Matrix& other_matrix)
{
    if (m_ncol != other_matrix.m_ncol || m_nrow != other_matrix.m_nrow)
        throw std::invalid_argument("Matrices have different dimensions!");

    for (size_t i{ 0 }; i < m_data.size(); i++)
    {
        m_data.at(i) -= other_matrix.m_data.at(i);
    }
}

void Matrix::operator*=(const double multiplier)
{
    for (size_t i{ 0 }; i < m_data.size(); i++)
    {
        m_data.at(i) *= multiplier;
    }
}

const Matrix Matrix::operator+(const Matrix& other_matrix) const
{
    if (m_ncol != other_matrix.m_ncol || m_nrow != other_matrix.m_nrow)
        throw std::invalid_argument("Matrices have different dimensions!");
    Matrix result{ *this };
    result += other_matrix;
    return result;
}

const Matrix Matrix::operator-(const Matrix& other_matrix) const
{
    if (m_ncol != other_matrix.m_ncol || m_nrow != other_matrix.m_nrow)
        throw std::invalid_argument("Matrices have different dimensions!");
    Matrix result{ *this };
    result -= other_matrix;
    return result;
}

const Matrix Matrix::operator*(const double multiplier) const
{
    Matrix result{ *this };
    result *= multiplier;
    return result;
}

const Matrix Matrix::operator*(const Matrix& other_matrix) const
{
    if (m_ncol != other_matrix.m_nrow)
        throw std::invalid_argument("Matrices are not compatible for multiplication!");

    size_t sum_over{ m_ncol };
    size_t nrow{ m_nrow };
    size_t ncol{ other_matrix.m_ncol };

    Matrix result{ nrow, ncol, std::vector<double>(nrow * ncol) };

    for (size_t i{ 0 }; i < nrow; i++)
    {
        for (size_t j{ 0 }; j < ncol; j++)
        {
            double element_value{ 0.0 };

            for (size_t k{ 0 }; k < sum_over; k++)
            {
                element_value += this->getElement(i, k) * other_matrix.getElement(k, j);
            }

            result.setElement(element_value, i, j);
        }
    }

    return result;
}

const Vector Matrix::operator*(const Vector& other_vector) const
{
    if (m_ncol != other_vector.m_nrow)
        throw std::invalid_argument("Matrices are not compatible for multiplication!");

    size_t sum_over{ m_ncol };
    size_t nrow{ m_nrow };

    std::vector<double> result(nrow);

    for (size_t i{ 0 }; i < nrow; i++)
    {
        double element_value{ 0.0 };

        for (size_t k{ 0 }; k < sum_over; k++)
        {
            element_value += this->getElement(i, k) * other_vector.getElement(k);
        }

        result.at(i) = element_value;
    }

    return Vector{ result };
}

double Matrix::operator[](const size_t idx) const
{
    if (idx >= m_nrow * m_ncol)
        throw std::invalid_argument("Index exceeds dimensions!");
    return m_data.at(idx);
}

Vector::Vector(const std::vector<double>& data)
    : Matrix{ data.size(), (data.size() == 0) ? 0 : 1, data }
{
}

double Vector::length() const
{
    if (m_data.size() == 0)
        throw std::invalid_argument("Vector is empty!");
    double sum_of_squares{ 0.0 };
    for (size_t i{ 0 }; i < m_data.size(); i++)
    {
        sum_of_squares += m_data.at(i) * m_data.at(i);
    }
    return std::sqrt(sum_of_squares);
}

double Vector::dotProduct(const Vector& other_vector) const
{
    if (m_data.size() != other_vector.m_data.size())
        throw std::invalid_argument("Vector don't have the same dimension!");
    else if (m_data.size() == 0)
        throw std::invalid_argument("Vectors are empty!");

    double sum{ 0.0 };
    for (size_t i{ 0 }; i < m_data.size(); i++)
    {
        sum += m_data.at(i) * other_vector.m_data.at(i);
    }
    return sum;
}

double Vector::distanceTo(const Vector& other_vector) const //const functions can be used by const class members
{
    if (m_data.size() != other_vector.m_data.size())
        throw std::invalid_argument("Vector don't have the same dimension!");
    else if (m_data.size() == 0)
        throw std::invalid_argument("Vectors are empty!");

    double sum_of_squares{ 0.0 };
    for (size_t i{ 0 }; i < m_data.size(); i++)
    {
        sum_of_squares += (m_data.at(i) - other_vector.m_data.at(i)) * (m_data.at(i) - other_vector.m_data.at(i));
    }
    return std::sqrt(sum_of_squares);
}

const Vector Vector::operator+(const Vector& other_vector) const
{
    if (m_nrow != other_vector.m_nrow)
        throw std::invalid_argument("Vectors have different dimensions!");
    Vector result{ *this };
    result += other_vector;
    return result;
}

const Vector Vector::operator-(const Vector& other_vector) const
{
    if (m_nrow != other_vector.m_nrow)
        throw std::invalid_argument("Vectors have different dimensions!");
    Vector result{ *this };
    result -= other_vector;
    return result;
}

const Vector Vector::operator*(const double multiplier) const
{
    Vector result{ *this };
    result *= multiplier;
    return result;
}

Answer 1

The code looks correct at first glance.

History

These days it is very hard to compete with any serious implementation of matrix operations and general BLAS. Blaze, Eigen, uBLAS, armadillo and the rest already put years into making their libraries create fast code for modern hardware. I will not mention any performance related aspects, as the submitted code doesn't seem to implement any optimizations.

I believe that every piece of code should be written with a purpose in mind. It is very hard to justify writing a BLAS library today without having any aspects that would be favorable compared to the competition.

Interface

The interface review will be from my computer vision background.

Not so compatible with OpenCV

Yes, I can write &m.at(0, 0), but I would much prefer to write m.data(). The latter follows the general convention of returning value_type*.

Only double

I cannot use this matrix to pass into any C related functions to read images, as they take char* or the like. The usual choice would be float rather than double for further processing.

Unusual function names

FYI, multidimensional subscript operator was voted into standard, so it might be worth it to provide the new subscript operator circa 2023. The old code uses function call operator to provide subscript support. There is also quite a bit of duplication in terms of functionality with no benefit. Not being able to use .at() on a const object would also upset quite a few users.

Unnecessary checks

It would be better to write assert rather than if. Nobody will use a BLAS library that does bounds checking on every access.

Unidirectional vector

I would like different names for column vector and row vector. Those are really matrices with one of the dimensions being set to 1.

No fixed size support

Providing a fixed size support would fix the unidirectional vector problem, as one dimension could be fixed and the rest be something like constexpr inline std::size_t dynamic_extent = std::numeric_limits<std::size_t>::max(); or straight up use the one from std.

---

If I would start a BLAS library today, I would probably write with SYCL or CUDA in mind. Those two are designed for this kind of data parallel computations in mind.

Answer 2

Design

Like @Incomputable mentioned you are not going to beat an existing matrix library out there. They are highly optimized and tuned for performance.

So lets look at this as an exercise.

I don't like that you have get and set methods for the members in the matrix. (ie getElement() and setElement()). This looks like it should all be done via a standard accesses operator.

  // These are just like your `at()` methods.
  // You should have normal operator accesses
  doubel&       operator()(int x, int y);
  double const& operator()(int x, int y) const;

I am using operator() just to make it simpler. Personally I would use operator[] but that takes slightly more work.

Note the difference here is that at() normally checks. While operator() does not check the bounds.

Think of the situaiton:

 Matrix    data(x,y, getData());
 for (int row = 0; row < data.nRow(); ++row) {
      for (int col = 0; col < data.nCol(); ++col) {
          data.at(row, col) += 12;
      }
 }

Here we are already guranteeing that the values (row/col) are in range (we guranteed by checking against nRow() and nCol() so there is no need to call at() and validate the row/col are good (we already guranteed this) so we should call a non checked interface.

Why do you need a separate Vector class?. Is this not just a Matrix with one dimension is only 1.

Review

If you are passing data that is going to become part of the object (like your constructor). Then the new consensus is to pass by value. This allows you to pass by l-Value and r-Value and achieve move where it is normal and copy where move is not an option.

So looking at your constructor:

Matrix::Matrix(const size_t nrow, const size_t ncol, const std::vector<double>& data)

Yes you can pass data here that is a temp object or a reference to some existing data and that's fine. But you are Always going to have a copy as the data is copied from the parameter into your internal structure.

It would be nice to rewrite this so you get a copy when you have no other choice but a move when it is available. The simple way is to make data a value parameter.

Matrix::Matrix(const size_t nrow, const size_t ncol, const std::vector<double> data)

But you are thinking hey now I am always paying the price of a copy. But think one second you are paying the price of a copy here at the interface but you can move it into the object (so you have move the location of the copy not added one) the advantage is that when sombody moves the data to your interface it can be moved to the parameter then moved it into the object so in some situations you have zero copies.

Matrix::Matrix(const size_t nrow, const size_t ncol, const std::vector<double> data)
    : m_nrow( nrow )
    , m_ncol( ncol )
    , m_data( std::move(data) )
{ /* STUFF */ }

usage:

 std::vector<double>  data{1,2,3,4}
 Matrix               m1(data);            // One Copy before and after change.
                                           // The location of the copy has moved.
 Matrix               m2(std::move(data)); // Now a move every step of the way.
 Matrix               m3(std::vector<double>{1,2,3}); // Was one copy before.
                                                      // But now is a guranteed move.

---

You have to have data in the array?

    if ((nrow * ncol) != data.size())
        throw std::invalid_argument("Data size does not match the dimensions!");

Can we not default initialize an array with all zero values?

---

I think you can simplify this test:

    if ((nrow == 0 && ncol > 0) || (nrow > 0 && ncol == 0))
        throw std::invalid_argument("One dimension is zero while the other one is not!");

Why not:

    if (nrow == 0 || ncol == 0)
        throw std::invalid_argument("At least one dimension is zero sized!");

---

But now you can reset the matrix to zero size.

void Matrix::clear()
{
    m_nrow = 0;
    m_ncol = 0;
    m_data.clear();
}

If you can not construct a Matrix of zero size why should you be allowed to reset the size to zero?

---

Why are you using at() here.

void Matrix::setElement(const double value, const size_t row_idx, const size_t col_idx)
{
    if (row_idx >= m_nrow || col_idx >= m_ncol)
        throw std::invalid_argument("Index exceeds dimensions!");
    m_data.at(row_idx * m_ncol + col_idx) = value;
}

You just validated that the values are in range and payed the cost of doing that. Why are you revalidating the value is in range! Do the validation once only use operator[].

---

Don't lile this.

void Matrix::setElement(const double value, const size_t idx)
{
    if (idx >= m_nrow * m_ncol)
        throw std::invalid_argument("Index exceeds dimensions!");
    m_data.at(idx) = value;
}

There is an implied need to understand the underlying implementation. You need to know that the data is represented as a single dimenstion array and that you have pre-multiplied the values together.

---

I would expect these to throw std::range_error

double& Matrix::at(const size_t row_idx, const size_t col_idx)
{
    if (row_idx >= m_nrow || col_idx >= m_ncol)
        throw std::invalid_argument("Index exceeds dimensions!");
    return m_data.at(row_idx * m_ncol + col_idx);
}

double& Matrix::at(const size_t idx)
{
    if (idx >= m_nrow * m_ncol)
        throw std::invalid_argument("Index exceeds dimensions!");
    return m_data.at(idx);
}

---

void Matrix::print() const

Sure this is fine. But it is nice to be able to send this to any stream (not just std::cout).

So I would add a parameter (that defaults to `std::cout)

void Matrix::print(std::ostream& stream = std::cout) const

Now that we have a stream based print() we can also add the operator<< that uses print. Thus allows us to print the Matrix using the standard stream operators.

 friend std::ostream& operator<<(std::ostream& stream, Matrix const& m)
 {
     m.print(stream);
     return stream;
 }

Now I can use the standard C++ way of outputing the Matrix:

Matrix  x(1, 15, getData());
std::cout << x << "\n\n";

---

The return value should NOT be const here:

const Matrix Matrix::operator+(const Matrix& other_matrix) const
const Matrix Matrix::operator-(const Matrix& other_matrix) const

I would write this like:

Matrix Matrix::operator-(const Matrix& other_matrix) const
{
    if (m_ncol != other_matrix.m_ncol || m_nrow != other_matrix.m_nrow)
        throw std::invalid_argument("Matrices have different dimensions!");
    Matrix result{ *this };
    return result -= other_matrix;  // compute and return
}

Answer 3

1. Firstly, Matrix and Vectors should not be parent or child classes. An ISA relationship is not enough to classify a parent and child relationship. If you really want to reuse code. It is better to use the vector code first and use it for composition.
2. Secondly using something like int a = 0 is more common than using something like int a = {0} and hence more readable.
3. Why have const return by values? I think replacing them with non-const by values is better.
4. length(): That function name is confusing, better change it to something like magnitude to avoid confusing with the size of the vector.