Re: [eigen] patch for tridiagonal matrices in eigen

["Gael Guennebaud" <gael.guennebaud@xxxxxxxxx>]

Hi,

wow that sounds really great. I have not much to say.

About the shortcuts, yes Upper must remain unchanged (Upper == UpperTriangularBit), and I'm also OK to rename it to UpperTriangular.

About the vectorization of DiagonalProduct, actually some cases was already vectorized for a long time. Also note that it does not have the EvaluateBeforeAssign bit, so it can already works in-place. What is missing is the vectorization of the other cases with all the mess we discussed last week... Now we have a TridiagonalProduct as well, I also think that's better to go for a C-style function with in-place evaluation.

Marco: in ei_traits<Product<TridiagonalProduct> >::Flags you should remove the VectorizableBit flag, we will vectorize it later using an auxiliary C-style function.

Also, I think we need two different implementations of TridiagonalProduct: one when we take the tridiagonal part of a complete matrix (or when we mark a complete matrix as tridiagonal), and another one when the tridiagonal matrix relies on a sparse storage (set of 3 vectors). Currently it relies on dynamic branching in TridiagonalProduct::coeff(). That's ok for fixed-size because our meta-unrollers will allow the compiler to remove them, but that's pretty bad for large sizes. Actually, an easy and generic way to implement it would be to use a sum of 3 diagonal products, but that would not be cache friendly, so...

I also think the TridiagonalMatrix _expression_ should overload the diagonal() function (to return a vector) and provide the respective sub-diagonal access functions.

We probably also need to add a specialization of assign to optimize simple operations on tridiagonal matrices (sum, etc.), more like we did for symmetric matrices.

Then there is the story of tridiagonal symmetric matrices....

We should also come up with a mechanism to define TriagonalMatrix with their own storage. For instance with DiagonalMatrix, one can currently do:

typedef DiagonalMatrix<NestByValue<VectorXf> > DiagonalMatrixXf;

to almost get a "real" diagonal matrix of float with dynamic size. I say "almost", because we cannot set its size, resize it, etc... Well one can add take a vector, and add a ".asDiagonal()" all the time, but, 1) that's not super convenient, and 2) this would be even less convenient for TridiagonalMatrix, so we really need to come up with something here.

To conclude, I'm obviously 100% OK that you commit your stuff whenever you can, but I'm not in favor to add such new beta features in Core. So I would suggest to put the support for tridiagonal matrix in the Sparse module or into a new module that we could call Experimental to give us the time we need to refine it a bit. What do you think ?

About the QR module, I think that MatrixType is already square, isn't it ? so I don't really follow you... About the "if", don't bother too much since this need to be redesigned anyway (see the other thread), but your use case (the matrix is already tridiagonal) will help to make the right choices...


cheers,
Gael.

On Mon, Dec 15, 2008 at 2:59 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

>>> -const unsigned int Upper = UpperTriangularBit;
>>> -const unsigned int StrictlyUpper = UpperTriangularBit | ZeroDiagBit;
>>> -const unsigned int Lower = LowerTriangularBit;
>>> -const unsigned int StrictlyLower = LowerTriangularBit | ZeroDiagBit;
>>> +const unsigned int Upper = UpperTriangularBit | UpperHessenbergBit;
>>> +const unsigned int StrictlyUpper = Upper | ZeroDiagBit;
>>> +const unsigned int Lower = LowerTriangularBit | LowerHessenbergBit;
>>> +const unsigned int StrictlyLower = Lower | ZeroDiagBit;
>>>
>>> here I don't really follow. These are parameters for part(). So unless
>>> you want to add Hessenberg support in part (which would be a good
>>> idea) no need to change that. If you want to add Hessenberg support,
>>> you need to make the corresponding changes in Part.h.
>>
>> I saw that the product helper checked for diagonals using the flags,
>> so I did the same for tridiagonals.
>
> OK right, sorry. These are always convenient shortcuts. What needs to
> be updated then is the comment saying that these are the parameters
> for part(), actually they're not used only for part().

Wait, no, I answered too fast. We really can't do such a thing as

+const unsigned int Upper = UpperTriangularBit | UpperHessenbergBit;

for example how would we then in part() distinguish between upper
triangular and upper hessenberg...

the real problem is that in the first place we allowed ourselves this
silly shortcut Upper for upper triangular. Now that there is upper
hessenberg too, this is getting ambiguous. I think we just rename
Upper to UpperTriangular and then we add UpperHessenberg alongside.

Yes this is an incompatible API change, there's a reason why we
haven't released yet. KDE / KOffice is not affected, as far as I can
remember.

Cheers,
Benoit

---