Re: [eigen] patch for tridiagonal matrices in eigen

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

Hi Mauro and Gael,

Finally I started actually thinking about this discussion (until now
my attitude was "noooo resist entering into yet another technical
discussion you don't have time")

First let's repeat one comment about the DiagonalMatrix class: Its
coeff() method does an if(). So it's slow (except in the fixed-size
cases where the if() can go away when loops are unrolled, but that's
not the general case). In 99% of practical use cases, one doesn't want
and one doesn't need to pay that price. For example, the diagonal
product is rewritten in such a way that this coeff() method is
bypassed.

That makes me think that the DiagonalMatrix expression was a bad idea
(my bad, it dates back to early eigen2 days) in the first place. Right
now I think it's completely useless, read on to see what I propose to
have instead.

First let's look at what happened since this DiagonalMatrix class was written:
1) Gael wrote the Diagonal product code that actually completely
bypasses the DiagonalMatrix class (i.e. the coeff() method). In Gael's
DiagonalProduct, the DiagonalMatrix is just used as a trivial wrapper
around the diagonal vector (which is directly used).
2) Mauro is now puzzled as to how to implement TridiagonalMatrix
because the wrong design of DiagonalMatrix is even worse when applied
to tridiagonal matrices.

Before continuing, let's state a rule to determine when an expression
class (here, DiagonalMatrix) is a bad idea (of course it's a
sufficient condition, not a necessary condition):

An expression class MyExpression is a bad idea whenever the
MyExpression::coeff() method is slow and in practice it can be
bypassed when implementing actual algorithms.

For exemple, as we said above, in the case of DiagonalMatrix, the slow
coeff() method is bypassed when implementing the diagonal matrix
product.

What do I propose to do then?

If MyExpression represents a special class of matrices for which the
storage can be optimized and certain algorithms can be optimized, then
instead of an expression class MyExpression, I think that what we want
is a special matrix class MySpecialMatrix that does NOT inherit
MatrixBase, so it's not at all an expression. MySpecialMatrix should
have only the storage that it needs, and should be where we implement
the special optimized algorithms.

For example the DiagonalMatrix class, rethought as a "special matrix",
would look like:

template<typename DiagonalVector>
class DiagonalMatrix
{
  public:
    template<typename OtherDerived>
    DiagonalProductReturnType<blabla> operator*(const
MatrixBase<OtherDerived>& other);

  protected:
    DiagonalVector m_diagonal;
};

Here it remains to be seen what DiagonalProductReturnType<blabla> is.
I think that we can make it a plain expression (i.e. inheriting
MatrixBase) because there the coeff() method can be implemented
without an if(). Not sure though, how that plays with vectorization
(depending on the matrix storage order). For TridiagonalProduct
though, it's not possible to make it return an expression (again
because the coeff() would do an if()). Probably we want some meta
logic here, like we have in Product.h, to distinguish cases where we
want an expression and cases where we don't.

You probably are going to ask: am I talking about a solution for the
diagonal-matrix-expression or for the diagonal-matrix-storage? The
answer is that the DiagonalMatrix above works for both! If you take
DiagonalMatrix<VectorXd> you have a diagonal-matrix-storage. On the
other hand, the diagonal-matrix-expression returned by
MatrixBase<Derived>::asDiagonal() would be a DiagonalMatrix<const
Derived&>. The user is not expected to construct that himself, he just
calls asDiagonal().

Another aspect of my proposal is that Product.h does no longer have to
dispatch all the special products, they are now implemented in the
special matrix classes. So it's more extensible and potentially
lighter on the compiler.

What do you think?

Cheers,
Benoit

---