Re: [eigen] patch for tridiagonal matrices in eigen

["Mauro Iazzi" <mauro.iazzi@xxxxxxxxx>]

Hi Benoit,
thanks for the comments. I'll apply the changes as soon as Gael
comments on them and I adapt them to the correct coding style (do you
have guidelines written somewhere?).

>>  - in SelfAdjointEigenSolver: use the _MatrixType::SquareMatrixType
>> for storing the eigenvalues, so that it can be instantiated with e.g.
>> _MatrixType == TridiagonalMatrix
>>  - skip tridiagonalization if the matrix is already tridiagonal (not
>> very happy with this, since it's a runtime check if compiler does not
>> optimize).
>
> I comment below, this should really use a meta-selector instead of a 'if'.

or maybe there's a better way.

> Do you mean adding new storage types, or just new expression / product
> types as you just did?
> If  you mean new product types, sure it's always good to make Eigen
> more extensible, but looking at Product.h it already looks quite
> extensible to me. Maybe you can elaborate after your experience?

I meant to define a new matrix class and specializing expressions,
like product. It was very easy to do, so it's already very extensible,
but I had to change the headers for that. Ideally I could just define
without having to put new code in the library.


> Actually I don't think that it should be difficult at all to add,
> e.g., a storage version of TridiagonalMatrix.
> The best solution, if possible, would be to just say that the storage
> version is just a partial specialization of TridiagonalMatrix where
> the three vector types are actual storage vector types (not abstract
> expressions). I haven't looked closely enough at that to be sure that
> it makes sense. That would probably require a flag forcing the coeffs
> vectors to be stored by value. Then some neat typedefs like
> TridiagonalMatrixXd so it's mostly transparent to the user. Or default
> to storage by value, I don't know...

I expressed myself badly, sorry. However I think that should work.

> Now some comments on your patch:

> +  protected:
> +    const typename CoeffsVectorType::Nested m_coeffs;
> +    const typename UCoeffsVectorType::Nested m_ucoeffs;
> +    const typename DCoeffsVectorType::Nested m_dcoeffs;
>
> Hm... apparently you use the letter D for "down" ? Very slippery as it
> could mean diagonal.
> So I suggest more explicit names: Diagonal,Upper,Lower...

very wise. I'll do that.

> -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.

> -  TridiagonalizationType::decomposeInPlace(m_eivec, diag, subdiag,
> computeEigenvectors);
> +  if ( (Flags & Tridiagonal) == Tridiagonal ) {
> +         for (int i=0;i<n-1;i++) {
> +                 diag[i] = matrix(i, i);
> +                 subdiag[i] = matrix(i, i+1);
> +         }
> +         diag[n-1] = matrix(n-1, n-1);
> +         m_eivec = EigenVectorMatrixType::Identity(n, n);
> +  } else {
> +         TridiagonalizationType::decomposeInPlace(m_eivec, diag, subdiag,
> computeEigenvectors);
> +  }
>
> Here, the compiler will have to go through the whole logic of BOTH
> cases, which means uselessly long compilation times. Yes the compiler
> will always optimize the "if" itself (even in debug mode) so it's not
> a runtime speed issue, nor a code size issue, it's really a
> compilation-time issue, and we're extremely careful about that.

understood. Maybe switching to MatrixType::SquareMatrixType also
inside Tridiagonalization<> would allow to specialize it for the
tridiagonal case. This would avoid storing and checking the flags and
also the manual copy of the vectors, since we could use the ones that
are already there (m_diagonalCoeffs, etc...).
> And, the line
>  if ( (Flags & Tridiagonal) == Tridiagonal ) {
> can be written as just
>  if(Flags & Tridiagonal) {

I think that would incorrectly get Upper/Lower Hessenberg as
Tridiagonals... I've seen (Flags & Diagonal) == Diagonal in other
places for the same reason.

cheers,

mauro

---