[eigen] Special Matrices: first get diagonal matrices done

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

Hi,

Before getting started with BandedMatrix and other kinds of special
matrices, I think the perfect starting point is to first nail the
situation with diagonal matrices.

Currently we have:
* class Diagonal nests a matrix xpr and is a vector xpr
* class DiagonalMatrixWrapper nests a vector xpr and is a matrix xpr
* class DiagonalMatrix stores a vector and is a matrix xpr

I think that we agree that we want instead DiagonalMatrix to NOT be a
matrix xpr (ie NOT inherit MatrixBase) and have a coeff() method that
just asserts row==col.

DiagonalMatrixWrapper is still as much needed as before for situations like

mat * vec.asDiagonal();

where asDiagonal() is actually just being used as a thin wrapper
around the vector allowing to select the matrix product strategy.

What is open for discussion is whether DiagonalMatrixWrapper should
still inherit MatrixBase? And if not, should there be an AsDense kind
of class constructing a MatrixBase expression out of a
DiagonalMatrixWrapper without evaluating?

Let me explain why I think that the answer to both questions is _no_.
The solution i'm going to advocate will be to keep the same class
names, only make DiagonalMatrixBase no longer inherit MatrixBase.

First of all we could have a DiagonalMatrixBase::toMatrix() that
returns a plain matrix (see at the end for naming discussion), so that
the user who wants something that behaves like a matrix, can get a
plain matrix if he wants to:

cout << vec.asDiagonal().toMatrix() << endl;

The argument for having instead a asDense() returning a MatrixBase, is
to avoid the cost of evaluating the matrix, but then remember that
this expression would have a very expensive coeff() method as it would
have to do a if(row==col). So as soon as you're accessing a large part
of the matrix, it's not good anyway (performance wise) to have an
asDense expression here.

I also think that there isn't much of a use case for doing directly
"vec.asDiagonal()(i,j)" -- it's OK, i think, to force the user to
replace that by "i==j?vec(i):0" making him realize where there is room
for optimization.

Then let's look at the main use case for asDiagonal():

mat * vec.asDiagonal()

Currently, since DiagonalMatrixWrapper is a MatrixBase, we have to
catch that in Product.h, adding complexity in there.

If instead DiagonalMatrixWrapper were NOT a MatrixBase, we could catch
that in a completely separate operator*. Less code overall and
especially in Product.h, shorter compilations times, "you don't pay
for what you don't use", same result.

Is there anything I didn't think of?

Yes that means that one could no longer do:
mat * mat2.marked<DiagonalBits>()
and hope for that product to be optimized, but frankly !! if you know
that a matrix happens to be diagonal, then you can probably avoid
storing it as a full matrix in the first place! So I don't think
there's a real use case for that; and even if there were one could
still do:
mat * mat2.diagonal().asDiagonal().

To summarize my plan:
*no more changes to Diagonal (the ex-DiagonalCoeffs)
*make DiagonalMatrixBase no longer inherit MatrixBase. It doesn't even
need a coeff() method.
*DiagonalMatrix still inherits DiagonalMatrixBase
*DiagonalMatrixWrapper still inherits DiagonalMatrixBase
*DiagonalMatrixBase is the class against which we implement the
operator* catching diagonal products, and by the way we can also have
operator= etc.
*DiagonalMatrixBase::toMatrix() returns a plain Matrix.

(the last method toMatrix(), i couldn't find a nicer name: eval()
suggests it returns a DiagonalMatrix, toDense() doesn't make it clear
it evaluates, ...)

Notice that if you are really in favor of a diagonal matrix
expression, we can still add this to this scheme using multiple
inheritance on MatrixBase and DiagonalMatrixBase...

Cheers,
Benoit