Generic matrix type in C++

Question

This snippet is about a generic matrix type. The nice part is that given two matrices with respective entry types, say, of short and float, the result matrix entry type after the multiplication will be float.

matrix.h:

#ifndef MATRIX_H
#define MATRIX_H

#include <algorithm>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <vector>

template<typename T>
class matrix {
    std::vector<std::vector<T>> content;
    size_t width;
    size_t height;

public:
    matrix(size_t width_, size_t height_) : width{width_}, height{height_}
    {
        content.resize(height);

        for (size_t i = 0; i < height; ++i)
        {
            content[i].resize(width);
        }
    }

    const std::vector<T>& operator[](size_t row_index) const {
        return content[row_index];
    }

    std::vector<T>& operator[](size_t row_index) {
        return content[row_index];
    }

    size_t get_width()  const { return width; }
    size_t get_height() const { return height; }
};

template<typename T1, typename T2>
auto operator*(const matrix<T1>& a, const matrix<T2>& b) -> matrix<decltype(T1{} * T2{})>
{
    if (a.get_height() != b.get_width())
    {
        std::stringstream ss;
        ss << "Matrix dimenstion mismatch: ";
        ss << a.get_height();
        ss << " x ";
        ss << a.get_width();
        ss << " times ";
        ss << b.get_height();
        ss << " x ";
        ss << b.get_width();
        ss << ".";
        throw std::runtime_error(ss.str());
    }

    using value_type = decltype(T1{} + T2{});
    matrix<decltype(T1{} * T2{})> result(a.get_height(), b.get_width());

    for (size_t rowa = 0; rowa != a.get_height(); ++rowa)
    {
        for (size_t colb = 0; colb != b.get_width(); ++colb)
        {
            value_type sum = 0;

            for (size_t i = 0; i != a.get_width(); ++i)
            {
                sum += a[rowa][i] * b[i][colb];
            }

            result[rowa][colb] = sum;
        }
    }

    return result;
}

template<typename T>
std::ostream& operator<<(std::ostream& os, matrix<T> m)
{
    size_t maximum_entry_length = 0;

    for (size_t row = 0; row < m.get_height(); ++row)
    {
        for (size_t col = 0; col < m.get_width(); ++col)
        {
            std::stringstream ss;
            ss << m[row][col];
            std::string entry_text;
            ss >> entry_text;
            maximum_entry_length = std::max(maximum_entry_length,
                                            entry_text.length());
        }
    }

    for (size_t row = 0; row < m.get_height(); ++row)
    {
        for (size_t col = 0; col < m.get_width(); ++col)
        {
            os << std::setw((int) maximum_entry_length) << m[row][col];

            if (col < m.get_width() - 1)
            {
                os << ' ';
            }
        }

        if (row < m.get_height() - 1)
        {
            os << '\n';
        }
    }

    return os;
}

#endif // MATRIX_H

main.cpp:

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

using std::cerr;
using std::cout;
using std::endl;

int main(int argc, const char * argv[]) {
    matrix<int> A(3, 2);

    A[0][0] = 1;
    A[0][1] = 2;
    A[0][2] = 3;
    A[1][0] = 4;
    A[1][1] = 5;
    A[1][2] = 6;

    cout << "A: " << endl;
    cout << A << endl << endl;

    matrix<float> B(2, 3);

    B[0][0] = 5.0f;
    B[0][1] = 4.0f;
    B[1][0] = 3.0f;
    B[1][1] = 2.0f;
    B[2][0] = 1.0f;
    B[2][1] = 0.0f;

    cout << "B: " << endl;
    cout << B << endl << endl;

    auto result = A * B;
    cout << "AB:" << endl;
    cout << result << endl << endl;

    result = B * A;
    cout << "BA:" << endl;
    cout << result << endl;

    matrix<short> F(4, 4);
    try
    {
        F * A;
    }
    catch (std::runtime_error& err)
    {
        cerr << err.what() << endl;
    }
}

As always, please tell me anything that comes to mind.

Answer 1

Use one vector instead of two vectors

You can improve performance by using one vector instead of two vectors; you can calculate the 2nd dimension using a little math. This is actually how primitive 2D arrays work in C++.

Advantages:

• Elements stay in contiguous storage, as opposed to linked-list-like structure where every row leads to an entirely different memory (in terms of adjacency). This is beneficial for the cache.
• Easy to implement

Disadvantages:

• You can no longer use the idiomatic matrix[i][j] syntax, thus you must provide your own access operator.

Possible implementation (please read the comments, they contain additional review items):

#include <cstddef> // added missing header for std::size_t
#include <vector>

// this macro can be replaced by a private (nearly guaranteed to be inlined) function
#define matrix_index(i, j) i * columns + j

template<class T>
class matrix
{
public:
    // use type aliases so that users can correctly refer to values
    using value_type = typename std::vector<T>::value_type;
    using reference = typename std::vector<T>::reference;
    using size_type = typename std::vector<T>::size_type;

    matrix(size_type const r, size_type const c)
        : rows{ r }
        , columns{ c }
    {
        // rows * colums = total memory
        data.resize(r * c);
    }

    reference operator()(size_type const i, size_type const j) noexcept
    {
        return data[matrix_index(i, j)];
    }

    /* ... */        

    // you can make these public and const, they never change
    size_type const rows;
    size_type const columns;

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

Sample usage:

for (matrix<int>::size_type i{ 0 }; i != m.rows; ++i)
{
    for (matrix<int>::size_type j{ 0 }; j != m.columns; ++j)
    {
        m(i, j) = i * j;
    }
}

Provide a constructor that directly initializes the matrix

Currently, you have to call resize(), which simply allocates and default constructs the number of specified elements. This is inefficient unless you want a default-initialized matrix.

Provide a constructor that allows you initialize the matrix directly:

matrix(std::initializer_list<std::initializer_list<value_type>> row_list)
    : rows{ row_list.size() }
    , columns{ row_list.begin()->size() }
{
    data.reserve(rows * columns);
    for (auto const& row : row_list)
    {
        data.insert(data.cend(), row);
    }
}

This means that you can now construct matrices directly like so:

matrix<int> m
{
    { 1, 2, 3 }, // row 1
    { 4, 5, 6 }  // row 2
};

Please note that this is a trivial implementation.

You should verify that the internal std::initializer_list<> instances all have the same size (or fill in the missing values with a default value such as 0). You can turn this into a "no-overhead" check by simply using assert() and having the NDEBUG macro defined in release mode.