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.
- 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?
- 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.
- You don’t support different character or traits types.
- … 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:
- 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.)
- 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.)
- You don’t need to worry about synchronization in the case of multi-threading.
- 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.
