Re: [eigen] patch for tridiagonal matrices in eigen

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

Hi Mauro,

Thank you for your interest in Eigen, and your patch. It's very nice
to see how you managed to understand and appropriate yourself  the
code -- reminds me of when Gael joined.

Your patch is actually very timely because we were just a few days ago
discussing fixes and improvements in diagonal product and Gael just
committed some fixes. Basically, the performance of diagonal product
could be much improved as it is not yet vectorized; we were still
wondering whether we should use a hardcoded function like Gael's
CacheFriendlyProduct, or vectorize the diagonal Product expression.

[[ Gael: i just saw that you put packet() methods in DiagonalProduct,
so does this mean that at least some cases of diagonal product are
already getting vectorized? ]]

Now your patch looks very good so I suggest that we quickly commit it
and hereafter we take tridiagonal product together with diagonal
product in consideration when doing future improvements.

> To get this I had to apply the following patch. What it does:
>  - define bits for upper/lower hessenberg matrices, and a Tridiagonal flag
>  - define a TridiagonalMatrix class, copied from DiagonalMatrix and
> define TridiagonalProduct accordingly

sounds good

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

> With this changes it was faster than lapack for that task. No actual
> number since I did not benchmark properly. In any case even if it was
> on par it would still be a very good result at this stage.

Impressive, as this part is not yet vectorized.

> Maybe some of this changes could be applied to the tree.

Sure, i want all of these improvements in the tree, and sooner rather
than later! As soon as you can apply the minor improvements that I
request below, and Gael says OK (he is more a specialist than me of
the "products" files in Eigen).

> I think that
> the possibility to define new matrix types (block diagonal comes to
> mind) could be very useful, so making this easy would be very nice.

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?

If you mean new storage types:

> One could use tridiagonals to test how difficult it is to add a new
> storage type to 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...

Now some comments on your patch:

1) please use 2 spaces for indentation
2) please try stay close to our "coding style" except if you really
hate it, e.g. we don't put spaces around parentheses.
  (though there is sometimes some inconsistency between Gael and me already)
3) please note that my current e-mail address is jacob.benoit.1@xxxxxxxxx
4) give yourself a copyright line on any new file you create.

+    /* FIXME: I don't know what it should do
+    template<int LoadMode>
+    const PacketScalar packet(int row, int col) const
+    {
+      if (RhsIsTridiagonal)
+      {
+        ei_assert((_LhsNested::Flags&RowMajorBit)==0);
+        return ei_pmul(m_lhs.template packet<LoadMode>(row, col),
ei_pset1(m_rhs.coeff(col, col)));
+      }
+      else
+      {
+        ei_assert(_RhsNested::Flags&RowMajorBit);
+        return ei_pmul(ei_pset1(m_lhs.coeff(row, row)),
m_rhs.template packet<LoadMode>(row, col));
+      }
+    }
+    */

You did the right thing: let that commented out for now.
As I said, we have yet to decide how to best vectorize these diagonal
/ tridiagonal products.

The packet(row,col) method should return the 128-bit packet (so 4
floats or 2 doubles) starting at the given (row,col) (so it depends on
the storage order).

Gael: for diagonal we were hesitating between a C-style function or
expressions for vectorization; now for tridiagonal product it seems
like the expressions way would be a nightmare, no? So that's one more
argument in favor of C-style functions. Plus, as I noted, for the
diagonal product, a C-style function could work in-place.... not sure
but i think that the tridiagonal case could still partly work in place
in the sense that the amount of temporary buffer memory needed would
be small compared to the memory used by whole matrices...

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

Then m_coeffs could become m_diagonalCoeffs (and then we'd have
m_upperCoeffs and m_lowerCoeffs).

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

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

The solution is to have a meta-selector, you can see tons of examples
throughout Eigen, just grep for '_selector" or "_impl". Basically it's
just using template specialization as a "early compile-time if" that
really allows the compiler to only spend time in the path that's
actually chosen.

Don't worry if this is unclear, we can handle that for you.

And, the line
  if ( (Flags & Tridiagonal) == Tridiagonal ) {
can be written as just
  if(Flags & Tridiagonal) {
it's actually less error-prone in case of future changes (more local,
only 1 place to change)...

Cheers,
Benoit

---