Linear algebra module

Question

I'm working on a linear algebra module to improve my knowledge with mathematics and, because I'll need a lightweight linear algebra module for my future work with Vulkan!

I tried to keep a blas-like routine style, for exemple, the related functions are called like this:

template <typename T, std::size_t M, std::size_t N>
Matrix<T, N, M> transposed(const DenseBase<T, M, N>& mat)

The whole module is templated and works with compile-time size matrices. I use an std::vector as a storage container because I want to keep it allocated on the heap (for large matrices)

The 2nd goal is to provide an access to row/column element by reference, for example:

Matrix<float, 10, 10> matrix{};
matrix.col(1)[1] = 5; // 2nd element of 2nd matrix's column is now set to 5

To achieve this, I have a hierarchy with a base class:

template <typename T, std::size_t M, std::size_t N>
class DenseBase
{
 public:
 using Matrixbase = Matrix<T, M, N>;
 using Row = Rowview<DenseBase, T, N>;
 using Col = Colview<DenseBase, T, M>;

 using ConstRow = const Rowview<const DenseBase, const T, N>;
 using ConstCol = const Colview<const DenseBase, const T, M>;

  Row row(std::size_t index) { return Row{*this, index}; }
  Col col(std::size_t index) { return Col{*this, index}; }

  ConstRow row(std::size_t index) const { return ConstRow{*this, index}; }
  ConstCol col(std::size_t index) const { return ConstCol{*this, index}; }
  T& operator()(std::size_t row, std::size_t col) { return coeff(row, col); }
  T& operator[](std::size_t index) { return coeff(index); }

  const T& operator()(std::size_t row, std::size_t col) const { return coeff(row, col); }
  const T& operator[](std::size_t index) const { return coeff(index); }

  virtual T& coeff(std::size_t row, std::size_t col) = 0;
  virtual const T& coeff(std::size_t row, std::size_t col) const = 0;

  virtual T& coeff(std::size_t index) = 0;
  virtual const T& coeff(std::size_t index) const = 0;
};

Then, I have a Rowview, Colview class which takes a Densebase as an argument and implements the coeff functions:

template <typename Parent, typename T, std::size_t N>
class Rowview : public DenseBase<T, 1, N>
{
 public:
  Rowview(Parent& parent, std::size_t row) : parent_{parent}, row_{row} {}

  T& coeff(std::size_t, std::size_t col) { return parent_(row_, col); }
  T& coeff(std::size_t index) { return parent_(row_, index); }

  const T& coeff(std::size_t, std::size_t col) const { return parent_(row_, col); }
  const T& coeff(std::size_t index) const { return parent_(row_, index); }

 private:
  Parent& parent_;
  std::size_t row_;
};

And finally, the matrix class itself, with the std::vector and the implementation of coeff methods to access its elements.

I really need some advice, maybe some ideas about how to improve the architecture. I don't like the idea of an abstract template class, but I have no alternatives yet.

Note that you can find the whole code on GitHub.

Answer 1

Since you are using matrices, I don't think that you need the runtime polymorphic behaviour offered by virtual but that you created virtual functions so that you could call methods of the derived class from methods of the base class. If this is the case, you could use static polymorphism instead, thanks to the Curiously Recurring Template Pattern idiom (CRTP).

In other words, make DenseBase take a Derived template parameter which corresponds to the type of the derived class:

template <typename T, std::size_t M, std::size_t N, typename Derived>
class DenseBase
{
    // ...
}

Then make the derived classes feed their own type to this template parameter:

template <typename Parent, typename T, std::size_t N>
class Rowview : public DenseBase<T, 1, N, Rowview<Parent, T, N>>
{
    // ...
}

Now, add a couple of methods to DenseBase to make a downcast when needed. Since you know the derived type at compile time thanks to the template parameter Derived, you can make a safe downcast with a static_cast:

Derived& derived() { return static_cast<Derived&>(*this); }
const Derived& derived() const { return static_cast<Derived&>(*this); }

Now you can rewrite the coeff methods in DenseBase so that they don't need to use virtual anymore to call the coeff methods from the derived class:

T& coeff(std::size_t row, std::size_t col) { return derived().coeff(row, col); }
const T& coeff(std::size_t row, std::size_t col) const { return derived().coeff(row, col); }

T& coeff(std::size_t index) { return derived().coeff(index); }
const T& coeff(std::size_t index) const { return derived().coeff(index); }

That way, you have some of the benefits of polymorphism without having to pay the cost of runtime polymorphism which I feel you won't need for such a class. That's a little bit tricky to understand and use at first, but it has clear benefits.