make access to off diagonal entries of DiagonalMatrix return zero

There is probably some other idea here that I'm missing but having the DiagonalMatrix return the diagonal value for all entries in a row makes generic programming difficult. Sometimes using mv etc is not an option and I needed access to jacobian[i][j]. Is the idea then to always copy jacobian into a DenseMatrix, which in almost all cases will be an overhead except if jacobian is a DiagonalMatrix or is there some other concept I'm missing? PS: perhaps bugfix should be feature...

I like the idea in principle -- it is surprising that you cannot simply read the off-diagonal entries.

zero_ should probably be static const -- otherwise the lifetime may be shorter than expected:

const DiagonalMatrix<...> mat = ...;
const auto &entry = mat[0][1];
std::cout << entry << std::endl; // error: the life-time of the object referred-to by entry has ended

In your current implementation, it only works for const matrices, which will certainly be surprising. C++17 as_const() can help, but only to make the workaround that users have to apply a little less cumbersome, it would still be surprising. An alternative would be to return a non-const reference to a static object, and declare that modifying that object is undefined behavior. std::string and its terminating '\0' would be a precedent for that.

Hey, I'm a c++98 guy (I just about understand the issue with the static :-) Seriously, your right, the implementation was just a quick hack that demonstrated my problem and solved it for me. Before adding anything more fancy I wanted to see if some change like this was wanted. Removing the operator[] from (at least the non const version) was an approach I was also thinking of perhaps adding a diag(i) method instead for the cases where one knows it's a DiagonalMatrx. At least then the compiler complains

@andreas.dedner: DiagonalMatrix is simply a sparse matrix (with a compile time pattern). Thus you can write generic code using iterators and their index()-method. When I originally implemented ScaledIdentityMatrix and DiagonalMatrix I also tried adding a fake-zero. But the compiler did not optimize this away in general and using the iterator was faster (as expected).

There's several thinks I don't like about this proposal:

• Being allowed to access, but not to modify entries in the non-const version
• Having non-persistent entries
• 'Existing' entries that are not visited by the iterator
• Being able to access the zero entries gives the impression that it's OK to loop over them, while using the iterators gives better performance in general.

To sum up, I think that the current situation is less confusing to users. But if we really add those fake-zeros, we should do the same for BCRSMatrix for consistency.

Good points and some of those were the reason why I didn't really make a full implementation yet. But why are the operator[] methods there at all? Using the mv methods etc or the iterator is the efficient way of doing things anyway - having the methods that then fails at runtime (and only if the assert is enabled) just doesn't seems a good idea. What is the application where one would need dm[5][5] but not dm[5][6]?

It's the same reason as for BCRSMatrix: Writing the entries. E.g., if one assembles a PDE where the local matrices are diagonal, then the operator[] allows simply write the non-zero entries the same way you do it for a FieldMatrix without doing any specialization.

Another reason is consistency. If other sparse matrices (BCRSMatrix) have operator[], then this one should have it, too.

Just to be clear: I consider removing the operator[] and replacing it by e.g. diag() an improvement as well. I'm undecided which of the two possibilities is better.

Maybe I should have highlighted this earlier: There is already DiagonalMatrix::diagonal(size_type). If you know that you're dealing with DiagonalMatrix that's IMO the preferred way for mutable access. But still the other way should exist. I'd object to removing operator[] here, because there's really no difference to BCRSMatrix.

I don't think the connection with BCRSMatrix is so obvious when this is used as the return value of jacobian etc on a grid geometry. I think it is surprising for a user that this would be sparse. I at least wouldn't expect it to behave like a BCRSMatrix but like a dense matrix in that setting (at least if that operator exists).

It's perhaps two different use cases we are thinking about here and those should perhaps have two different implementations?

mentioned in issue #145

mentioned in merge request !953 (closed)

mentioned in issue #253

Can we close this? There's many arguments against this change. On top of the above given ones one can add that is wastes memory. If one really needs a dense implementation this should be done separately.

mentioned in issue #333 (closed)