[eigen] patch for tridiagonal matrices in eigen

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

Hi all,
  first of all thanks to the authors for the library. It's good,
useful and simple.

I am a phd student in physics and I will have to work with electronic
structure simulations, so a lot of linear algebra will be involved. I
briefly tested eigen and it was faster than lapack, which is the main
(and strongly suggested) alternative.

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

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.

Maybe some of this changes could be applied to the tree. 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.
One could use tridiagonals to test how difficult it is to add a new
storage type to the library.

Cheers and thanks again for eigen,

mauro

diff -urN Eigen-svn/Core Eigen-new/Core
--- Eigen-svn/Core	2008-12-14 00:04:04.000000000 +0100
+++ Eigen-new/Core	2008-12-13 23:45:48.000000000 +0100
@@ -40,8 +40,6 @@
 #include <complex>
 #include <cassert>
 #include <functional>
-#include <iostream>
-#include <iostream>
 #include <cstring>
 #include <string>
 
@@ -105,6 +103,7 @@
 #include "src/Core/Transpose.h"
 #include "src/Core/DiagonalMatrix.h"
 #include "src/Core/DiagonalCoeffs.h"
+#include "src/Core/TridiagonalMatrix.h"
 #include "src/Core/Sum.h"
 #include "src/Core/Redux.h"
 #include "src/Core/Visitor.h"
diff -urN Eigen-svn/src/Core/Product.h Eigen-new/src/Core/Product.h
--- Eigen-svn/src/Core/Product.h	2008-12-14 00:03:04.000000000 +0100
+++ Eigen-new/src/Core/Product.h	2008-12-13 23:43:30.000000000 +0100
@@ -79,6 +79,7 @@
  *    - NormalProduct
  *    - CacheFriendlyProduct
  *    - DiagonalProduct
+ *    - TridiagonalProduct
  *    - SparseProduct
  */
 template<typename Lhs, typename Rhs> struct ei_product_mode
@@ -87,6 +88,8 @@
 
     value = ((Rhs::Flags&Diagonal)==Diagonal) || ((Lhs::Flags&Diagonal)==Diagonal)
           ? DiagonalProduct
+          : ((Rhs::Flags&Tridiagonal)==Tridiagonal) || ((Lhs::Flags&Tridiagonal)==Tridiagonal)
+          ? TridiagonalProduct
           : (Rhs::Flags & Lhs::Flags & SparseBit)
           ? SparseProduct
           : Lhs::MaxColsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
diff -urN Eigen-svn/src/Core/TridiagonalMatrix.h Eigen-new/src/Core/TridiagonalMatrix.h
--- Eigen-svn/src/Core/TridiagonalMatrix.h	1970-01-01 01:00:00.000000000 +0100
+++ Eigen-new/src/Core/TridiagonalMatrix.h	2008-12-13 23:48:16.000000000 +0100
@@ -0,0 +1,96 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2006-2008 Benoit Jacob <jacob@xxxxxxxxxxxxxxx>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_TRIDIAGONALMATRIX_H
+#define EIGEN_TRIDIAGONALMATRIX_H
+
+/** \class TridiagonalMatrix
+  *
+  * \brief Expression of a diagonal matrix
+  *
+  * \param CoeffsVectorType the type of the vector of diagonal coefficients
+  *
+  * This class is an expression of a diagonal matrix with given vector of diagonal
+  * coefficients. It is the return
+  * type of MatrixBase::diagonal(const OtherDerived&) and most of the time this is
+  * the only way it is used.
+  *
+  * \sa MatrixBase::diagonal(const OtherDerived&)
+  */
+template<typename CoeffsVectorType, typename UCoeffsVectorType, typename DCoeffsVectorType>
+struct ei_traits<TridiagonalMatrix<CoeffsVectorType, UCoeffsVectorType, DCoeffsVectorType> >
+{
+  // TODO: add assertions for coefficient consistency
+  typedef typename CoeffsVectorType::Scalar Scalar;
+  typedef typename ei_nested<CoeffsVectorType>::type CoeffsVectorTypeNested;
+  typedef typename ei_unref<CoeffsVectorTypeNested>::type _CoeffsVectorTypeNested;
+  enum {
+    RowsAtCompileTime = CoeffsVectorType::SizeAtCompileTime,
+    ColsAtCompileTime = CoeffsVectorType::SizeAtCompileTime,
+    MaxRowsAtCompileTime = CoeffsVectorType::MaxSizeAtCompileTime,
+    MaxColsAtCompileTime = CoeffsVectorType::MaxSizeAtCompileTime,
+    Flags = (_CoeffsVectorTypeNested::Flags & HereditaryBits) | Tridiagonal,
+    CoeffReadCost = _CoeffsVectorTypeNested::CoeffReadCost
+  };
+};
+
+template<typename CoeffsVectorType, typename UCoeffsVectorType, typename DCoeffsVectorType>
+class TridiagonalMatrix : ei_no_assignment_operator,
+   public MatrixBase<TridiagonalMatrix<CoeffsVectorType, UCoeffsVectorType, DCoeffsVectorType> >
+{
+  public:
+
+    EIGEN_GENERIC_PUBLIC_INTERFACE(TridiagonalMatrix)
+
+    inline TridiagonalMatrix(const CoeffsVectorType& coeffs, const UCoeffsVectorType& ucoeffs, const DCoeffsVectorType& dcoeffs) : m_coeffs(coeffs), m_ucoeffs(ucoeffs), m_dcoeffs(dcoeffs)
+    {
+      ei_assert(CoeffsVectorType::IsVectorAtCompileTime
+          && coeffs.size() > 0);
+      ei_assert(DCoeffsVectorType::IsVectorAtCompileTime
+          && dcoeffs.size() > 0 && (coeffs.size()==ucoeffs.size()+1));
+      ei_assert(UCoeffsVectorType::IsVectorAtCompileTime
+          && ucoeffs.size() > 0 && (coeffs.size()==ucoeffs.size()+1));
+    }
+
+
+    inline int rows() const { return m_coeffs.size(); }
+    inline int cols() const { return m_coeffs.size(); }
+
+    inline const Scalar coeff(int row, int col) const
+    {
+      switch (row-col) {
+          case 0: return m_coeffs.coeff(row); break;
+          case 1: return m_ucoeffs.coeff(col); break;
+          case -1: return m_dcoeffs.coeff(row); break;
+	  default: return static_cast<Scalar>(0); break;
+      };
+    }
+
+  protected:
+    const typename CoeffsVectorType::Nested m_coeffs;
+    const typename UCoeffsVectorType::Nested m_ucoeffs;
+    const typename DCoeffsVectorType::Nested m_dcoeffs;
+};
+
+#endif // EIGEN_TRIDIAGONALMATRIX_H
diff -urN Eigen-svn/src/Core/TridiagonalProduct.h Eigen-new/src/Core/TridiagonalProduct.h
--- Eigen-svn/src/Core/TridiagonalProduct.h	1970-01-01 01:00:00.000000000 +0100
+++ Eigen-new/src/Core/TridiagonalProduct.h	2008-12-13 23:30:59.000000000 +0100
@@ -0,0 +1,117 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@xxxxxxx>
+// Copyright (C) 2006-2008 Benoit Jacob <jacob@xxxxxxxxxxxxxxx>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_TRIDIAGONALPRODUCT_H
+#define EIGEN_TRIDIAGONALPRODUCT_H
+
+template<typename LhsNested, typename RhsNested>
+struct ei_traits<Product<LhsNested, RhsNested, TridiagonalProduct> >
+{
+  // clean the nested types:
+  typedef typename ei_unconst<typename ei_unref<LhsNested>::type>::type _LhsNested;
+  typedef typename ei_unconst<typename ei_unref<RhsNested>::type>::type _RhsNested;
+  typedef typename _LhsNested::Scalar Scalar;
+
+  enum {
+    LhsFlags = _LhsNested::Flags,
+    RhsFlags = _RhsNested::Flags,
+    RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
+    ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
+    MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
+    MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
+
+    LhsIsTridiagonal = (_LhsNested::Flags&Tridiagonal)==Tridiagonal,
+    RhsIsTridiagonal = (_RhsNested::Flags&Tridiagonal)==Tridiagonal,
+
+    CanVectorizeRhs =  (!RhsIsTridiagonal) && (RhsFlags & RowMajorBit) && (RhsFlags & PacketAccessBit)
+                     && (ColsAtCompileTime % ei_packet_traits<Scalar>::size == 0),
+
+    CanVectorizeLhs =  (!LhsIsTridiagonal) && (!(LhsFlags & RowMajorBit)) && (LhsFlags & PacketAccessBit)
+                     && (RowsAtCompileTime % ei_packet_traits<Scalar>::size == 0),
+
+    RemovedBits = ~((RhsFlags & RowMajorBit) && (!CanVectorizeLhs) ? 0 : RowMajorBit),
+
+    Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
+          | (CanVectorizeLhs || CanVectorizeRhs ? PacketAccessBit : 0),
+
+    CoeffReadCost = NumTraits<Scalar>::MulCost + _LhsNested::CoeffReadCost + _RhsNested::CoeffReadCost
+  };
+};
+
+template<typename LhsNested, typename RhsNested> class Product<LhsNested, RhsNested, TridiagonalProduct> : ei_no_assignment_operator,
+  public MatrixBase<Product<LhsNested, RhsNested, TridiagonalProduct> >
+{
+    typedef typename ei_traits<Product>::_LhsNested _LhsNested;
+    typedef typename ei_traits<Product>::_RhsNested _RhsNested;
+
+    enum {
+      RhsIsTridiagonal = (_RhsNested::Flags&Tridiagonal)==Tridiagonal
+    };
+
+  public:
+
+    EIGEN_GENERIC_PUBLIC_INTERFACE(Product)
+
+    template<typename Lhs, typename Rhs>
+    inline Product(const Lhs& lhs, const Rhs& rhs)
+      : m_lhs(lhs), m_rhs(rhs)
+    {
+      ei_assert(lhs.cols() == rhs.rows());
+    }
+
+    inline int rows() const { return m_lhs.rows(); }
+    inline int cols() const { return m_rhs.cols(); }
+
+    const Scalar coeff(int row, int col) const
+    {
+      const int unique = RhsIsTridiagonal ? col : row;
+      return m_lhs.coeff(row, unique) * m_rhs.coeff(unique, col)
+              + m_lhs.coeff(row, unique-1) * m_rhs.coeff(unique-1, col)
+              + m_lhs.coeff(row, unique+1) * m_rhs.coeff(unique+1, col);
+    }
+
+    /* 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));
+      }
+    }
+    */
+
+  protected:
+    const LhsNested m_lhs;
+    const RhsNested m_rhs;
+};
+
+#endif // EIGEN_TRIDIAGONALPRODUCT_H
diff -urN Eigen-svn/src/Core/util/Constants.h Eigen-new/src/Core/util/Constants.h
--- Eigen-svn/src/Core/util/Constants.h	2008-12-14 00:03:04.000000000 +0100
+++ Eigen-new/src/Core/util/Constants.h	2008-12-13 23:38:53.000000000 +0100
@@ -178,6 +178,16 @@
   * means the expression includes sparse matrices and the sparse path has to be taken. */
 const unsigned int SparseBit = 0x1000;
 
+/** \ingroup flags
+  *
+  * means the strictly lower triangular part is 0 */
+const unsigned int UpperHessenbergBit = 0x2000;
+
+/** \ingroup flags
+  *
+  * means the strictly upper triangular part is 0 */
+const unsigned int LowerHessenbergBit = 0x4000;
+
 // list of flags that are inherited by default
 const unsigned int HereditaryBits = RowMajorBit
                                   | EvalBeforeNestingBit
@@ -185,23 +195,26 @@
                                   | SparseBit;
 
 // Possible values for the Mode parameter of part() and of extract()
-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;
 const unsigned int SelfAdjoint = SelfAdjointBit;
+const unsigned int UpperHessenberg = UpperHessenbergBit;
+const unsigned int LowerHessenberg = LowerHessenbergBit;
 
 // additional possible values for the Mode parameter of extract()
 const unsigned int UnitUpper = UpperTriangularBit | UnitDiagBit;
 const unsigned int UnitLower = LowerTriangularBit | UnitDiagBit;
 const unsigned int Diagonal = Upper | Lower;
+const unsigned int Tridiagonal = UpperHessenberg | LowerHessenberg;
 
 enum { Aligned, Unaligned };
 enum { ForceAligned, AsRequested };
 enum { ConditionalJumpCost = 5 };
 enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
 enum DirectionType { Vertical, Horizontal };
-enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct, SparseProduct };
+enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct, TridiagonalProduct, SparseProduct };
 
 enum {
   /** \internal Equivalent to a slice vectorization for fixed-size matrices having good alignment
diff -urN Eigen-svn/src/Core/util/ForwardDeclarations.h Eigen-new/src/Core/util/ForwardDeclarations.h
--- Eigen-svn/src/Core/util/ForwardDeclarations.h	2008-12-14 00:03:04.000000000 +0100
+++ Eigen-new/src/Core/util/ForwardDeclarations.h	2008-12-13 23:47:38.000000000 +0100
@@ -46,6 +46,7 @@
 template<typename Lhs, typename Rhs, int ProductMode> class Product;
 template<typename CoeffsVectorType> class DiagonalMatrix;
 template<typename MatrixType> class DiagonalCoeffs;
+template<typename CoeffsVectorType, typename UCoeffsVectorType, typename DCoeffsVectorType> class TridiagonalMatrix;
 template<typename MatrixType, int PacketAccess = AsRequested> class Map;
 template<typename MatrixType, unsigned int Mode> class Part;
 template<typename MatrixType, unsigned int Mode> class Extract;
diff -urN Eigen-svn/src/QR/SelfAdjointEigenSolver.h Eigen-new/src/QR/SelfAdjointEigenSolver.h
--- Eigen-svn/src/QR/SelfAdjointEigenSolver.h	2008-12-14 00:03:04.000000000 +0100
+++ Eigen-new/src/QR/SelfAdjointEigenSolver.h	2008-12-13 23:57:33.000000000 +0100
@@ -42,14 +42,15 @@
 {
   public:
 
-    enum {Size = _MatrixType::RowsAtCompileTime };
+    enum { Size = _MatrixType::RowsAtCompileTime, Flags = _MatrixType::Flags };
     typedef _MatrixType MatrixType;
+    typedef typename _MatrixType::SquareMatrixType EigenVectorMatrixType;
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<Scalar>::Real RealScalar;
     typedef std::complex<RealScalar> Complex;
     typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
     typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
-    typedef Tridiagonalization<MatrixType> TridiagonalizationType;
+    typedef Tridiagonalization<EigenVectorMatrixType> TridiagonalizationType;
 
     SelfAdjointEigenSolver()
         : m_eivec(Size, Size),
@@ -94,7 +95,7 @@
     void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true);
 
     /** \returns the computed eigen vectors as a matrix of column vectors */
-    MatrixType eigenvectors(void) const
+    EigenVectorMatrixType eigenvectors(void) const
     {
       #ifndef NDEBUG
       ei_assert(m_eigenvectorsOk);
@@ -106,7 +107,7 @@
     RealVectorType eigenvalues(void) const { return m_eivalues; }
 
   protected:
-    MatrixType m_eivec;
+    EigenVectorMatrixType m_eivec;
     RealVectorType m_eivalues;
     #ifndef NDEBUG
     bool m_eigenvectorsOk;
@@ -177,7 +178,16 @@
   // (same for diag and subdiag)
   RealVectorType& diag = m_eivalues;
   typename TridiagonalizationType::SubDiagonalType subdiag(n-1);
-  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);
+  }
 
   int end = n-1;
   int start = 0;