7
\$\begingroup\$

I have written a minimal matrix class. I would like to utilize the latest C++ features, follow best practices, and utilize multi-threading in some of the operations. I have also setup a GoogleTest environment with test cases like verifying matrix initialization and calculations but I am not sure what other test cases are meaningful. I am using this project to gain more understanding in C++, OOP, testing, multi-threading, and basically writing good code.

I want some feedback on what I am doing right, what I am doing wrong, and what I can improve on.

I will be updating my code on github

Compiled with clang++ -pedantic-errors -Wall -Weffc++ -Wextra -Wsign-conversion -std=c++2b matrix.hpp.

#include <cassert>
#include <cmath>
#include <cstddef>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <random>
#include <stdexcept>
#include <vector>

using Vector = std::vector<float>;
using Index = Vector::size_type;
class Matrix {
private:
    Index rows_{};
    Index cols_{};
    Vector elements_{};
    
public:
    Matrix(Index rows = 0, Index cols = 0, float f = 0.0f) : rows_{ rows }, cols_{ cols }, elements_{ Vector(rows*cols, f) } {
        assert(rows >= 0);
        assert(cols >= 0);

    }
    Matrix(Index rows, Index cols, Vector elements) : rows_{ rows }, cols_{ cols }, elements_{ elements } {
        assert(rows >= 0);
        assert(cols >= 0);
    }

    float& operator[] (Index index) {
        assert(index >= 0);
        if (index >= elements_.size())
            throw std::out_of_range("operator[]: Out of range.\n");

        return elements_.at(index);
    }

    float& operator[] (Index row, Index col) {
        assert(row >= 0);
        assert(col >= 0);
        if (row*cols_ + col >= elements_.size())
            throw std::out_of_range("operator[]: Out of range.\n");
        return elements_.at(row*cols_ + col);
    }

    friend std::ostream& operator<<(std::ostream& out, const Matrix& matrix) {
        assert(matrix.elements_.size() >= 0);
        out << std::fixed;
        out << std::setprecision(4);

        Index rows{ matrix.rows_ };
        Index cols{ matrix.cols_ };

        out << '[';
        for(Index i{ 0 }; i < rows; ++i)
        {
            if (i != 0)
                out << ' ';
            out << '[';

            for(Index j{ 0 }; j < cols; ++j)
            {
                out << matrix.elements_.at(i*cols+j);
                if (j != (cols - 1))
                    out << ' ';
            }
            out << ']';
            if (i != (rows - 1))
                out << '\n';
        }
        out << ']';
        return out;

    }

    friend Matrix operator+(const Matrix& m, const Matrix& n) {
        assert(m.rows_ == n.rows_);
        assert(m.cols_ == n.cols_);
        Vector u = m.elements_;
        Vector v = n.elements_;
        std::transform(u.cbegin(), u.cend(), v.cbegin(), u.begin(), std::plus<float>());

        return Matrix{ m.rows_, m.cols_, u };
    }

    friend Matrix operator+(const Matrix& m, float value) {
        Vector v = m.elements_;
        std::transform(v.cbegin(), v.cend(), v.begin(), [value](float f){ return f + value; });

        return Matrix{ m.rows_, m.cols_, v };
    }

    friend Matrix operator-(const Matrix& m, const Matrix& n) {
        assert(m.rows_ == n.rows_);
        assert(m.cols_ == n.cols_);
        Vector u = m.elements_;
        Vector v = n.elements_;
        std::transform(u.cbegin(), u.cend(), v.cbegin(), u.begin(), std::minus<float>());

        return Matrix{ m.rows_, m.cols_, u };
    }

    friend Matrix operator-(const Matrix& m, float value) {
        Vector v = m.elements_;
        std::transform(v.cbegin(), v.cend(),v.begin(), [value](float f){return f - value;});

        return Matrix{ m.rows_, m.cols_, v };
    }

    Matrix operator-() const {
        Matrix negativeMat = *this;
        std::transform(negativeMat.elements_.cbegin(), negativeMat.elements_.cend(), negativeMat.elements_.begin(), std::negate<float>());

        return negativeMat;
    }

    void alloc(Vector elements) {
        assert(rows_*cols_ == std::size(elements));
        elements_ = elements;
    }

    auto capacity() {
        return elements_.capacity();
    }

    Matrix dot(Matrix m) {
        Matrix dst(rows_, m.cols_);
        assert(cols_ == m.rows_);
        Index size{ cols_ };
        assert(dst.rows_ == rows_);
        assert(dst.cols_ == m.cols_);

        for (Index i{ 0 }; i < dst.rows_; ++i) {
            for (Index j{ 0 }; j < dst.cols_; ++j) {
                dst.elements_.at(i*dst.cols_+j) = 0.0f;
                for (Index k{ 0 }; k < size; ++k) {
                    dst.elements_.at(i*dst.cols_+j) += elements_.at(i*cols_+k) * m.elements_.at(k*m.cols_+j);
                }
            }
        }
        return dst;
    }

    Matrix col(Index c) {
        Matrix dst(rows_, 1);
        for (Index i{ 0 }; i < rows_; ++i)
            dst.elements_.at(i) = elements_.at(c + i*cols_);

        return dst;
    }

    Index getCols() {
        return cols_;
    }

    Index getRows() {
        return rows_;
    }

    void rand(float low, float high) {
        std::mt19937 mt{ std::random_device{}() };
        std::uniform_real_distribution rand(0.0, 1.0);

        for (Index i{ 0 }; i < rows_*cols_; ++i)
            elements_[i] = (rand(mt)*(high-low)+low);
    }

    Matrix row(Index r) {
        Matrix dst(1, cols_);
        Index begin_offset{ r*cols_ };
        Index end_offset{ (r+1)*cols_ };
        Vector slice(elements_.begin()+static_cast<long>(begin_offset), elements_.begin()+static_cast<long>(end_offset));
        dst.alloc(slice);

        return dst;
    }

    void sigmoid() {
        std::for_each(elements_.begin(), elements_.end(), [](float &f){ f = 1.0f / (1.0f + std::exp(-f)); });
    }

    auto size() {
        return elements_.size();
    }
};
\$\endgroup\$

4 Answers 4

8
\$\begingroup\$

Use std::mdspan

Since you tagged the post , consider using std::mdspan in your code. Not only can you use that internally to simplify accessing a given element using its row and column index, you can also have row() and col() return a std::mdspan, which avoids making copies.

Keep it simple

There are several functions that can be removed from the class and just be completely defined outside it. For example, both the binary and unary operator-() can be defined outside of the class. They don't need to be friends; they can read and modify the elements of a Matrix object just fine using other public functions.

This especially goes for more complex functions like rand() and sigmoid(). Where does it stop? What if someone wants to do a matrix transpose or inverse, are you going to add those as member functions as well? What about all of linear algebra?

Make it more generic

What if you need a matrix of doubles? What if you need a three-dimensional matrix? Consider making your class more generic. Making the type of elements you store generic is easy:

template<typename T>
class Matrix {
    using value_type = T;
    using container_type = std::vector<value_type>;
    using size_type = container_type::size_type;

    container_type elements_{};
public:
    Matrix(size_type rows, size_type cols, value_type value = {}): … {}
    …
};

Making the number of dimensions generic is possible as well, using std::dextents for example. Also have a look at the proposed std::mdarray.

\$\endgroup\$
4
  • \$\begingroup\$ Depending on what you mean, either do not support "3D matrices" (they're tensors), or use Matrix<Matrix<double>> could also be an approach. \$\endgroup\$ Dec 7, 2023 at 13:46
  • \$\begingroup\$ Thanks for your insight! Eventually, I want to do linear algebra operations and create a neural network framework based on this. If that is the case, should I create a separate header file dedicated for linear algebra operations? \$\endgroup\$ Dec 7, 2023 at 22:40
  • 1
    \$\begingroup\$ Yes, I would indeed separate that into its own header file. \$\endgroup\$
    – G. Sliepen
    Dec 7, 2023 at 23:13
  • \$\begingroup\$ A matrix is 2D by definition. There are no 3D matrices. My upvote is for where you encourage free functions over making everything be members. \$\endgroup\$ Dec 8, 2023 at 2:18
8
\$\begingroup\$

Missing all includes and the C++ version, so I had to make some guesses to get this to compile.

    assert(rows >= 0);
    assert(cols >= 0);

OK, but Index comes from std::vector<T>::size_type which is required to be an unsigned integer type. So these asserts do nothing, but they shouldn't be harmful.

Allowing a dimension to be zero does not seem useful to me.

rows*cols could wrap, resulting in a vector elements_ that cannot hold as many elements as expected, and no exception.

The constructor that takes a vector of elements as input does nothing to ensure that the number of elements is correct.

    if (row*cols_ + col >= elements_.size())
        throw std::out_of_range("operator[]: Out of range.\n");
    return elements_.at(row*cols_ + col);

This does do something, but maybe not quite what you wanted. Checking only the combined index means that col may be out of range and the index may still be accepted. Additionally, row*cols_, or row*cols_ + col, may wrap if any of the values are sufficiently large. The resulting index may be less than elements_.size(), that way certain out-of-bounds indexing may succeed. Using at effectively repeats the same range test that you just performed manually so I don't know what the added value is.

The last point also applies to the other indexer, the one with a one-dimensional index. Also, I don't think you should have that indexer, it seems to me that it invites the simple error of forgetting the specify either the row or column. It would be safer to offer that style of indexing via a plain old method (if you offer it at all), since it could not be confused with the indexing operator.

dot

Probably the most important operation, so let's look at it in more detail:

Matrix dot(Matrix m) {
    Matrix dst(rows_, m.cols_);
    assert(cols_ == m.rows_);
    Index size{ cols_ };
    assert(dst.rows_ == rows_);
    assert(dst.cols_ == m.cols_);

    for (Index i{ 0 }; i < dst.rows_; ++i) {
        for (Index j{ 0 }; j < dst.cols_; ++j) {
            dst.elements_.at(i*dst.cols_+j) = 0.0f;
            for (Index k{ 0 }; k < size; ++k) {
                dst.elements_.at(i*dst.cols_+j) += elements_.at(i*cols_+k) * m.elements_.at(k*m.cols_+j);
            }
        }
    }
    return dst;
}

assert(dst.rows_ == rows_) and assert(dst.cols_ == m.cols_) are very paranoid. The dst matrix has just been created with those dimensions. If it turns out to not have those dimensions now, you're dealing with a problem against which there is no defense.

Using at here should be unnecessary. You know how big the matrices are and your loops by design stay in the bounds. Using at unnecessarily doesn't make the code safer, but does make it harder for the compiler to generate good code.

Accumulating the result in the matrix element (instead of in a local variable which is stored in the matrix in the end) makes it harder for the compiler to generate good code. Even if the compiler can understand what you're doing, expect at least a redundant store in every iteration.

Making matrix multiplication anywhere near efficient would take me a dozen or so paragraphs to explain, I can do it if you want but since you didn't tag this question I haven't bothered with it yet. Even keeping it single-threaded, you can make this two orders of magnitude faster (exact ratio depends on the size of the matrix and machine specifics). That's not even an exaggeration, the impact is really that large.

You've mentioned threading though, so I'll note that adding threading to fundamentally inefficient code won't be effective. The main problem with inefficient matrix multiplication is that it doesn't use the caches well and thereby makes cores sit around idle waiting for data from main memory. Doing that on more threads mostly makes more cores fight over the same memory bandwidth, though it may help for some sizes of matrix. By contrast, efficient matrix multiplication uses the caches very well and therefore can scale well with multiple threads since the cores wouldn't be fighting over a shared resource.

\$\endgroup\$
1
  • \$\begingroup\$ Thank you so much for your input. I changed my code according to your suggestions. I have included all the missing includes and the command I used to compile the code. \$\endgroup\$ Dec 7, 2023 at 7:25
2
\$\begingroup\$

Overview

If your Matrix sizes are always known at compile time then you could potentially have written them to do compile time checking of dimensions. Of course if you use any dynamic input this constraint will probably fall apart (but there are a lot of situations were static sizes are perfectly valid).


Don't like this:

float& operator[] (Index index);

This is a bug waiting to happen. I would like some time of compile time check to make sure that that this index is not being used incorrectly. It's nice that C++23 allows multiple parameters to operator[] that has solved the problem for you.


I think these are baddly defined.

    Matrix col(Index c) {
    Matrix row(Index c) {

There usage is not intuitive enough that they would always be used without error. I personally would make the return a view into the matrix.

Codereview

Float?

using Vector = std::vector<float>;

That could be very limiting on your use cases. Personally I would have made the Matrix class templateable with the type (and potentially defaulted to float).


Why are you making this hard to read?

public:
    Matrix(Index rows = 0, Index cols = 0, float f = 0.0f) : rows_{ rows }, cols_{ cols }, elements_{ Vector(rows*cols, f) } {

Same principle for member variable initialization as local variable initialization. All coding standards (I have read) say one variable per line (please).

public:
    Matrix(Index rows = 0, Index cols = 0, float f = 0.0f)
        : rows_{ rows }
        , cols_{ cols }
        , elements_{ Vector(rows*cols, f) }
    {

I saw somebody else also point this out.
But if either of these is zero the Matrix collapses down to empty.

        assert(rows >= 0);
        assert(cols >= 0);

Also note that if neither of these can be initialized with a negative value. If you pass a negative value then integer conversion will be done by the compiler and you will probably end up with an excessively large poitive value (that is not caught by the assert).


You should check the size of elements.

    Matrix(Index rows, Index cols, Vector elements) : rows_{ rows }, cols_{ cols }, elements_{ elements } {

Also that is a copy operation. You may want to move array content.

    Matrix(Index rows, Index cols, Vector elements)
        : rows_{ rows }
        , cols_{ cols }
        , elements_{ std::move(elements) }
    {
        assert(rows_ * cols_ == elements_.size());
    }

That's not how I would expect that to work!!!!

    float& operator[] (Index index) {
        assert(index >= 0);
        if (index >= elements_.size())
            throw std::out_of_range("operator[]: Out of range.\n");

        return elements_.at(index);
    }

I would expect that to return a row of elements.

I would change the input type so that you could only get the type by combining two Index values together. So that a user explicitly knows that they are retrieving a pre-computed value.

This is a bug waiting to happen.


Comments in-line:

    float& operator[] (Index row, Index col) {

        // These asserts do nothing.
        assert(row >= 0);
        assert(col >= 0);

        // This check in inadequate
        // There are still lots of invalid ways to accesses the
        // elements that pass this test.
        if (row*cols_ + col >= elements_.size())
            throw std::out_of_range("operator[]: Out of range.\n");

        // This is a superfluous use of `at()`.
        // You have already done the validation of rows/cols
        // No need to check it again!
        return elements_.at(row*cols_ + col);
    }

    friend std::ostream& operator<<(std::ostream& out, const Matrix& matrix) {
        // Why.
        // Did you not validate this in the constructor.
        // This can never be true.
        assert(matrix.elements_.size() >= 0);

        // You have modified the state.
        // Do you not want to reset it back after this function.
        // otherwise you are polluting the stream and it will
        // affect other users.
        out << std::fixed;
        out << std::setprecision(4);

    friend Matrix operator+(const Matrix& m, const Matrix& n) {
        assert(m.rows_ == n.rows_);
        assert(m.cols_ == n.cols_);

        // Sure you can make copy of this.
        // Assuming you also fix the move in the constructor.
        // But is making a copy the most efficient technique?
        //     I don't know the answer.
        //     But you could make an empty array and reserve the required
        //     size then use `std::transform` with `std::back_inserter` (maybe)
        Vector u = m.elements_;


        // Don't think you want to actually copy this.
        // I would rather take a reference to the original.
        // Or simply use the iterator directly.
        Vector v = n.elements_;



        std::transform(u.cbegin(), u.cend(), v.cbegin(), u.begin(), std::plus<float>());

        return Matrix{ m.rows_, m.cols_, u };
    }

I would also note: That operator+ is usually implemented in terms of operator+=. Which basically gives you both operations very simply.


The minus operators seem to be identical to the plus operators. With one function different. Might be worth implementing these in a common function and passing the function as a parameter.

    friend Matrix operator-(const Matrix& m, const Matrix& n);
    friend Matrix operator-(const Matrix& m, float value);

I think a move operation is in order here:

    void alloc(Vector elements) {
        assert(rows_*cols_ == std::size(elements));
        elements_ = elements;
    }

You are making a copy of the column here.

    Matrix col(Index c) {
        Matrix dst(rows_, 1);
        for (Index i{ 0 }; i < rows_; ++i)
            dst.elements_.at(i) = elements_.at(c + i*cols_);

        return dst;
    }

Not sure that is obvious. If I use it to get a column I may want to be able to use it to manipulate the original Matirx.

Same comment for:

    Matrix row(Index r) {

Personally I would make both these return a view into the original Matrix. If the user wants them become a separate Matrix simply allow them to assign the view to a matrix to make a copy.


The mt19937 is relatively expensive to create. I would create this once and re-use.

    void rand(float low, float high) {

        static std::mt19937 mt{ std::random_device{}() };
        static std::uniform_real_distribution rand(0.0, 1.0);
     // ^^^^^^ Add static so they are only created once.
     //        The first time this function is called they are
     //        created. All other times they will be re-used.

        for (Index i{ 0 }; i < rows_*cols_; ++i)
            elements_[i] = (rand(mt)*(high-low)+low);
    }

These should be const methods:

    Index getCols() {
        return cols_;
    }

    Index getRows() {
        return rows_;
    }


    auto size() {
        return elements_.size();
    }
\$\endgroup\$
2
\$\begingroup\$

I’m not going to do a full review, because the reviews you already have are excellent. I’ll just a few extra notes since you mentioned you are interested in next-gen, C++23-and-beyond best practices.

Type aliases

using Vector = std::vector<float>;
using Index = Vector::size_type;

Both of these aliases should be defined in the class, since they only make sense in context of the class. However, I think you could even go a few steps further.

Let’s start with the index type. Rather than using a single type for both row and column indices, it might be a better idea to create strong aliases to distinguish them. This will prevent this kind of error:

// I want a 2x3 matrix:
auto mat = Matrix{3, 2}; // oops!

With strong aliases, not only could you detect mixing up the parameters at compile time, you could even make that a feature:

auto mat1 = Matrix{Matrix::column_index_type{3}, Matrix::row_index_type{2}}; // works! (makes a 2x3 matrix)
auto mat2 = Matrix{Matrix::row_index_type{2}, Matrix::column_index_type{3}}; // also makes a 2x3 matrix

This looks like a lot more typing, but in practice it could be minimized in a dozen ways. Like using UDLs, for example:

auto mat = Matrix{2_row, 3_col};

There is also a standard type for dimensional extents: std::extents. In addition to the constructors that take row and column extents, you could have a constructor that takes a std::extents<IndexType, I1, I2> object, like so:

template <std::integral IndexType, std::size_t I1, std::size_t I2>
constexpr explicit Matrix(std::extents<IndexType, I1, I2> e)
    : Matrix{e.extent(0), e.extent(1)}
{}

(While you’re at it, you may also want to actually store the extents as a std::extents, and have a extents() function, just like mdspan.)

As for the vector type alias, I’d suggest getting rid of that entirely. It is only used 2 places in the public interface, and both are unwise. More on that in a bit.

Allocator awareness

Because your matrix uses dynamic allocation, it should support the standard allocator model. Unfortunately, this isn’t easy… but it doesn’t need to be extremely hard.

For starters, you could add an allocator template parameter. This is the “classic” way of adding allocator support. But a better, more modern way is to use polymorphic allocators. Then you don’t need to make the matrix class a template (though, as @G.Sliepen notes, you probably want to anyway).

Either way, you’d need an allocator_type nested type, and each constructor needs a duplicate that takes an allocator… even the special constructors (default, copy, move). Plus you need a couple other things (like a member swap()). There’s a fair bit of boilerplate, but it’s not too bad:

class Matrix
{
    using vector_type = std::pmr::vector<float>;

    vector_type::size_type rows_ = {};
    vector_type::size_type cols_ = {};
    vector_type elements_ = {};

public:
    using index_type     = vector_type::size_type;
    using allocator_type = vector_type::allocator_type;

    // Default construct:
    constexpr Matrix() noexcept = default;
    constexpr explicit Matrix(allocator_type const& a) noexcept
        : elements_{a}
    {}

    // Copy construct:
    constexpr Matrix(Matrix const&) = default;
    constexpr Matrix(Matrix const& m, allocator_type const& a)
        : elements_{m.elements_, a}
    {}

    // Move construct:
    constexpr Matrix(Matrix&&) noexcept = default;
    constexpr Matrix(Matrix&& m, allocator_type const& a)
        : elements_{std::move(m.elements_), a}
    {}

    // Custom constructors (basically just delegate each one or use default arguments):

    // Example of delegating:
    constexpr Matrix(index_type r, index_type c)
        : Matrix{r, c, allocator_type{}}
    {}
    constexpr Matrix(index_type r, index_type c, allocator_type const& a)
        : elements_(r * c, a}
    {}

    // Example with default argument:
    constexpr Matrix(index_type r, index_type c, float v, allocator_type const& a = {})
        : elements_(r * c, v, a}
    {}

    // Need this:
    constexpr auto get_allocator() const noexcept { return elements_.get_allocator(); }

    // Also need this (must be member function):
    constexpr auto swap(Matrix& m)
        noexcept(std::allocator_traits<allocator_type>::propagate_on_container_swap::value or std::allocator_traits<allocator_type>::is_always_equal::value)
    {
        return elements_.swap(m.elements_);
    }

    // Don’t forget assignment ops:
    constexpr auto operator=(Matrix const&) -> Matrix& = default;
    constexpr auto operator=(Matrix&&) noexcept(std::allocator_traits<allocator_type>::propagate_on_container_move_assignment::value or std::allocator_traits<allocator_type>::is_always_equal::value) -> Matrix& = default;

    // ... and the rest is basically what you have...
};

Better constructors

I am not a fan of default arguments, and this constructor pretty ably illustrates why:

    Matrix(Index rows = 0, Index cols = 0, float f = 0.0f) : rows_{ rows }, cols_{ cols }, elements_{ Vector(rows*cols, f) } {
        assert(rows >= 0);
        assert(cols >= 0);

    }

You may not have intended it, but you have made this possible:

Matrix m = 2; // whut?

Part of the problem is that the constructor isn’t explicit, but even if were you could write auto m = Matrix{2};… and that’s still pretty silly, because you can’t describe a matrix with a single number.

Default arguments save you some repetition and boilerplate, but they also have hidden costs, and create extra complexity (and, sometimes, as in the case above, absurdity).

You don’t want a constructor that can take 0, 1, 2, or 3 numbers. You want one constructor that takes 0 arguments, and one – a different one, that does something logically distinct – that takes 2 or 3. The zero argument constructor is not the same as calling the 2 (or 3) argument constructor with all zeroes. The zero argument default-initializes the matrix. The 2/3 argument constructor does a specific construction task (creating a matrix with the given size and default value for the elements). Yes, if you pass all zeroes to the 2/3 argument constructor you get effectively the same result… but by a very different mechanism; and this is easy to prove by pointing out that vector<float>{} is noexcept, while vector<float>(size, val) is not, even if you call it with zeroes.

So the above constructor should be broken into:

Matrix() noexcept = default;
// Note: No need for the asserts in the above. And it's noexcept.

// You could leave the last argument with a default, because
// Matrix{a, b} and Matrix{a, b, c} do the same logical operation
// (the former just uses a default value for c).
//
// Personally, though I prefer not to use default args at all, though.
Matrix(Index rows, Index cols, float f = 0.0f)
    : rows_{rows}
    , cols_{cols}
    , elements_(rows * cols, f) // SEE NOTE!
{
    assert(rows >= 0);
    assert(cols >= 0);
}

// NOTE:
//
// In case you didn't notice, I changed your vector initialization from:
//  Vector{rows * cols, f}
// to:
//  Vector(rows * cols, f)
//
// Notice that? I changed the braces to parentheses.
//
// The reason for this is a quirk in std::vector that causes ambiguity
// when giving one or two numeric constructor arguments. It's too
// complex to explain here; if you don't know about it, see the note
// on the cppreference page about vector's constructors:
//  https://en.cppreference.com/w/cpp/container/vector/vector
//
// Right *NOW*, your code is "safe", because the vector's size_type
// will *always* be distinct from float.
//
// However, if you follow the advice to make the matrix element type
// a template parameter, watch out. You could end up with some very
// nasty bugs.

The second constructor is for initializing a matrix with existing data… but right now it requires that data to be in a vector.

But… why? Why can’t the data be in an array? Or a span? Or… anything else?

This is a job for a template:

template <std::ranges::input_range R>
    requires std::convertible_to<ranges::range_reference_t<R>, float>
Matrix(Index rows, Index cols, R&& elements)
    : rows_{rows}
    , cols_{cols}
    , elements_(rows * cols, {})
{
    assert(rows >= 0);
    assert(cols >= 0);

    std::ranges::copy(
        std::forward<R>(elements) | std::views::take(elements_.size()),
        elements_.begin()
    );
}

In the above, if you are copying to a \$M\times N\$ matrix, but the source range has fewer than \$M\times N\$ elements, the remaining elements are value-initialized (which basically means zeroes). If the source range has more than \$M\times N\$ elements, the excess elements are ignored.

You may want to change this behaviour, as it is not really a great idea. It is always a bad idea to hide or ignore strange things; they are usually a sign of bigger problems. It makes more sense to only accept \$M\times N\$ elements. If there are more or less, you should throw an exception. (If someone really wants to allow source ranges with an arbitrary numbers of elements, and pad or trim as necessary, it is trivial to write a wrapper.)

IOStreams sucks

There are a lot of problems with your inserter. Like… a lot. They’re not your fault. They’re because IOStreams sucks.

  1. You don’t take into account the stream’s state. For example, what if the stream’s width has been set to 100? What will that do to your output? What if the locale is set to a locale where there are spaces in the number? For example, the French would write “1 234 567,89”. Won’t that muck up your matrix layout?
  2. You change the stream state, and leave it changed. This can be very bad, and can completely screw up the entire rest of a program’s runtime, in many if not most programs.
  3. You don’t support different character or traits types.
  4. … and so on….

There is one simple form for an inserter that solves most, if not all, of these problems, at the cost of a little extra inefficiency—which usually doesn’t matter much for output—and slightly less flexibility—such as not supporting locale-specific formatting, but you’re not using that anyway. Here is the basic structure:

class T
{
    // ... everything else in the class...

public:
    template <typename Char, typename Traits>
    friend auto operator<<(std::basic_ostream<Char, Traits>& out, T const& t) -> std::basic_ostream<Char, Traits>&
    {
        if (out)
        {
            auto oss = std::ostringstream{};
            oss.imbue(std::locale::classic());

            // DO ALL OUTPUT TO OSS HERE!

            out << std::move(oss).str().c_str();
        }

        return out;
    }
};

So in your case:

class Matrix
{
    // ... everything else in the class...

    template <typename Char, typename Traits>
    friend auto operator<<(std::basic_ostream<Char, Traits>& out, Matrix const& m) -> std::basic_ostream<Char, Traits>&
    {
        if (out)
        {
            auto oss = std::ostringstream{};
            oss.imbue(std::locale::classic());

            oss << std::fixed;
            oss << std::setprecision(4);

            Index rows{ matrix.rows_ };
            Index cols{ matrix.cols_ };

            oss << '[';
            for(Index i{ 0 }; i < rows; ++i)
            {
                if (i != 0)
                    oss << ' ';
                oss << '[';

                for(Index j{ 0 }; j < cols; ++j)
                {
                    oss << matrix.elements_.at(i*cols+j);
                    if (j != (cols - 1))
                        oss << ' ';
                }
                oss << ']';
                if (i != (rows - 1))
                    oss << '\n';
            }
            oss << ']';

            out << std::move(oss).str().c_str();
        }

        return out;
    }
};

Doing all output to a temporary string stream means:

  1. The initial state of the stream is always predictable. (It is usually smart to set the locale to the classic POSIX/C locale, because that’s what you’re probably already assuming you’re using.)
  2. It doesn’t matter if you change the stream state. (Who cares if you’ve changed the floating point output mode? Won’t affect anyone else.)
  3. You don’t need to worry about synchronization in the case of multi-threading.
  4. The whole output operation is seen by the final output stream as a single output operation. So all the formatting options work as expected (assuming they make any sense at all).

And there’s a special trick with IOStreams where you can always stream a NUL-terminated char string to any stream type (wchar_t, different traits types, whatever), and it will automatically be converted. So you just write your output to a plain, ordinary string string, then convert that to a NUL-terminated char string (via a C++ string, via .str()), and then stream that to your final output stream, and poof, it works with everything.

If you want to get fancy, though, you might consider writing a custom formatter for your matrix type, that allows customizing things like the row prefix/suffix, the element delimiter, the format for elements, and so on.

assert() is (almost) the past

Contracts didn’t make C++23, but there is a very good chance they will make C++26. Most notably, it looks like the syntax has finally been settled on.

Most, if not all, of the asserts in your code are not technically assertions… they are preconditions. They are things that should be true before any of your matrix code even starts.

For example:

    friend Matrix operator+(const Matrix& m, const Matrix& n) {
        assert(m.rows_ == n.rows_);
        assert(m.cols_ == n.cols_);
        Vector u = m.elements_;
        Vector v = n.elements_;
        std::transform(u.cbegin(), u.cend(), v.cbegin(), u.begin(), std::plus<float>());

        return Matrix{ m.rows_, m.cols_, u };
    }

Should be:

    friend auto operator+(Matrix const& m, Matrix const& n) -> Matrix
        // You could do this:
        //  pre (m.rows_ == n.rows_)
        //  pre (m.cols_ == n.cols_)
        // But it's better to define preconditions in terms of the
        // public interface, because it is the caller's responsibility
        // to satisfy them, and the caller shouldn't know the class
        // internals:
        pre (m.getRows() == n.getRows())
        pre (m.getCols() == n.getCols())
        // You could also add postconditions:
        //  post (r: r.getRows() == m.getRows())
        //  post (r: r.getCols() == m.getCols())
        // You can even have "expensive" pre/post conditions, that are
        // only checked in certain modes. (Exactly *how* to specify
        // a check is "audit-only" has not yet been finalized.)
        //  post (r: std::ranges::equal(std::views::zip_transform(std::plus<>{}, m.elements_, n.elements_), r.elements_))
    {
        // It would be wiser to define operator+=, then write this as:
        //  auto u = m;
        //  m += n;
        //  return u;
        // But I've left it more or less as is here, except simplified.

        auto u = m;

        std::ranges::transform(u.elements_, n.elements_, u.elements_.begin(), std::plus<>{});

        return u;
    }

Now, as I said, contracts aren’t in C++23… but they are coming. So you might as well plan for them. That means pulling out all assertions that are actually preconditions, moving them to the beginning of the function, and refactoring them in terms of the public interface. That will make transitioning to contract preconditions much easier.

You may also want to at least document postconditons, even though it makes little sense to use assert() to verify them. (assert() doesn’t really work well for postconditions, because it doesn’t really work with static analysis. It doesn’t hurt, it just wastes time.)

A BLAS from the future

Finally, as a last hint about the future: Right now it looks very likely that C++ will have a linear algebra library in C++26. This will not make your matrix class pointless; on the contrary, the standard linear algebra library will be very low level, and will not have a matrix type at all. It will just have matrix functions, that take mdspan. (You mentioned you were considering making a linear algebra library yourself. This is still a good idea, because you could make a very nice high-level, feature-rich interface, and use the low-level standard functions to do the actual number crunching.)

What that means is that @G.Sliepen’s suggestion to better incorporate mdspan into your matrix class is golden advice. If your matrix class can easily convert to a mdspan view of its elements… then you get all the benefits of the upcoming linear algebra library automatically. Which means, for example, you can ditch your dot() function, and instead write:

class Matrix
{
    // ... rest of class ...
    //
    // ... except no need for dot() member function.

public:
    constexpr auto view()      & noexcept { return std::mdspan{elements_, rows_, cols_}; }
    constexpr auto view() const& noexcept { return std::mdspan{elements_, rows_, cols_}; }
};

// And you get dot() for free!
auto m1 = Matrix{2, 3, {11, 12, 13, 21, 22, 23}};
auto m2 = Matrix{3, 2, {11, 12, 21, 22, 31, 32}};

auto result = std::linalg::dot(m1.view(), m2.view());

Or if you prefer:

class Matrix
{
    // ... rest of class ...
    //
    // ... except no need for dot() member function.

public:
    friend auto dot(Matrix const& m1, Matrix const& m2) noexcept
    {
        assert(m1.getRows() == m2.getCols());
        assert(m1.getCols() == m2.getRows());

        auto const v1 = std::mdspan{m1.elements_, m1.rows_, m1.cols_};
        auto const v2 = std::mdspan{m2.elements_, m2.rows_, m2.cols_};

        return std::linalg::dot(v1, v2);
    }
};

// Doing dot simpler:
auto m1 = Matrix{2, 3, {11, 12, 13, 21, 22, 23}};
auto m2 = Matrix{3, 2, {11, 12, 21, 22, 31, 32}};

auto result = dot(m1, m2);

Though you probably still want to have some kind of conversion functions to convert your matrices to and from mdspans easily, so you get all of the linear algebra functions, present and future. You can still add wrappers like the one above, and still add ergonomics like being able to add matrices like “r = a + b”, rather than “std::linalg::add(a, b, r);” (though you could implement that operator+ in terms of std::linalg::add()).

The bottom line is this: make your matrix class easily convertible to and from std::mdspan<T, std::dextents<std::size_t, 2>>, or similar (and std::mdspan<T const, std::dextents<std::size_t, 2>> for const correctness). It will be worth it.

\$\endgroup\$
2
  • \$\begingroup\$ I have a question... I think all constructors should require at least a row and a column parameter. Is this bad practice? And if it is fine, do I = delete the default constructor? \$\endgroup\$ Jan 3 at 23:34
  • \$\begingroup\$ It’s not bad practice at all. And if a default constructor doesn’t make sense, then by all means, don’t include one. (You probably don’t actually need to delete it, though. The default constructor is automatically deleted whenever you declare any other (non-default) constructor.) Bear in mind, though, that a class without a default constructor is generally a lot harder to use, which is why it is recommended that if you can make one, you should. But if you can’t, that’s fine. \$\endgroup\$
    – indi
    Jan 6 at 0:57

Your Answer

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

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