3
\$\begingroup\$

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.

\$\endgroup\$
7
  • \$\begingroup\$ @user2296177 Sure. But how would I implement the double subscript operator? \$\endgroup\$
    – coderodde
    Sep 29, 2016 at 16:35
  • \$\begingroup\$ @user2296177 Why not make an answer with the two points? \$\endgroup\$
    – coderodde
    Sep 29, 2016 at 16:40
  • 1
    \$\begingroup\$ I've just posted very similar code as an answer to a different question - I made the dimensions part of the type, and use modern decltype() and Concepts, but you may find it interesting anyway. \$\endgroup\$ Sep 29, 2016 at 17:09
  • \$\begingroup\$ My opinion is to have dimensions as part of the type. On top of that, operator()(std::size_t will need to return std::array<> in order to avoid users trying to resize it. Additionally, multiplication will issue compile time error if used with wrong dimensions. I think just storing raw pointer to std::array<std::array<>> is ok. \$\endgroup\$ Sep 29, 2016 at 17:19
  • 1
    \$\begingroup\$ I don't have time for a full answer but your matrix multiplication is slow for large matrices. There are more efficient methods with lower memory bandwidth, see Wikipedia for example: en.m.wikipedia.org/wiki/Matrix_multiplication_algorithm in general I advise against building your own math primitives for anything else than a learning experience. It's not worth the time and you will likely end up with something suboptimal any way. \$\endgroup\$
    – Emily L.
    Sep 29, 2016 at 20:56

1 Answer 1

7
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.