Bug / RFC: Expected behaviour of BCRSMatrix::axpy()
You can currently put all kinds of matrices into a BCRSMatrix
. When you're assembling a lumped mass matrix, maybe you'll want to use a ScaledIdentityMatrix
instead of a FieldMatrix
-- saves space and matrix-vector multiplication at the block level will be cheaper.
The BCRSMatrix
class also has a function axpy
which does m1 += c * m2
, i.e. it adds the c
-fold multiple of m2
to m1
. Patterns need to be compatible, naturally, and the way it is implemented, it hands over the work to the axpy
method of the blocks inside the loop.
Now here's the problem: Some matrix classes do not have an axpy
function. These include DiagonalMatrix
and ScaledIdentityMatrix
(IdentityMatrix
is already deprecated so I won't bother with pointing out issues that class had). So when you try m1.axpy(c, m2)
with blocks that are instances of the ScaledIdentityMatrix
class, you'll get an error at compile-time because there is no ScaledIdentityMatrix::axpy
while m1.axpy(c, m2)
will compile and work as expected when the blocks are of type FieldMatrix
instead.
My question is: Is all of that intentional? Should it be that way? I would expect BCRSMatrix<ScaledIdentityMatrix>::axpy(scalar, matrix)
to work when matrix is itself of type BCRSMatrix<ScaledIdentityMatrix>
. Should BCRSMatrix::axpy
fail differently? Should there be a ScaledIdentityMatrix::axpy()
?
Here's the code that I used for testing
#include "config.h"
#include <dune/common/diagonalmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/dynmatrix.hh>
#include <dune/common/identitymatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/scaledidmatrix.hh>
using namespace Dune;
template <class M>
struct MatrixProvider {
M static get();
};
template <>
struct MatrixProvider<Dune::DynamicMatrix<double>> {
Dune::DynamicMatrix<double> static get() { return {{1, 2}, {3, 4}}; }
};
template <>
struct MatrixProvider<Dune::FieldMatrix<double, 2, 2>> {
Dune::FieldMatrix<double, 2, 2> static get() { return {{1, 5}, {9, 16}}; }
};
template <>
struct MatrixProvider<Dune::IdentityMatrix<double, 2>> {
Dune::IdentityMatrix<double, 2> static get() { return {}; }
};
template <>
struct MatrixProvider<Dune::DiagonalMatrix<double, 2>> {
Dune::DiagonalMatrix<double, 2> static get() { return {25}; }
};
template <>
struct MatrixProvider<Dune::ScaledIdentityMatrix<double, 2>> {
Dune::ScaledIdentityMatrix<double, 2> static get() { return {11}; }
};
template <class T>
void test() {
using M = BCRSMatrix<T>;
T const a = MatrixProvider<T>::get();
M m(2, 2, BCRSMatrix<T>::random);
m.setrowsize(0, 2);
m.setrowsize(1, 2);
m.endrowsizes();
m.addindex(0, 0);
m.addindex(0, 1);
m.addindex(1, 0);
m.addindex(1, 1);
m.endindices();
for (auto &row : m)
for (auto &col : row)
col = a;
M m2 = m;
m2.axpy(2.5, m);
}
int main() {
test<FieldMatrix<double, 2, 2>>();
test<DynamicMatrix<double>>();
//test<DiagonalMatrix<double, 2>>(); // has no axpy
//test<IdentityMatrix<double, 2>>(); // has no blocklevel
//test<ScaledIdentityMatrix<double, 2>>(); // has no axpy
}