Access to off-diagonal entries of a DiagonalMatrix
I feel like this has been discussed before, but I cannot find anything in here or in Flyspray, so I bring this up again.
The current behaviour of a Dune::DiagonalMatrix
is such that the following snippet happily claims that all off-diagonal entries of A are 1. In a debug build it this is asserted.
DiagonalMatrix<double, 2> A(1);
for (int i=0; i<A.rows; ++i)
for (int j=0; i<A.cols; ++i)
if (i != j)
std::cout << "Offdiagonal entry: " << A[i][j] << std::endl;
I consider this a bug that should be fixed, even if it involves introducing the if
that checks for the entry being off-diagonal. Otherwise DiagonalMatrix
can hardly claim to be an implementation of the mathematical concept of a diagonal matrix and does not deserve its name.
For those interested, here is a scenario, where this is biting us heavily: We are generating finite element assembly code with different backends one of them being generic, in the sense that it may be used with all sorts of grids. However, it currently fails for YaspGrid, because the generated code does not (cannot!) make assumptions on the absence of off diagonal terms and computes garbage. Changing the code generator to avoid direct matrix accesses in favor of mv
methods etc. would heavily violate design principles of the DSL -> IR -> Code Toolchain.