From ee982709d37e219e9d67beaf777a5bf2131e835e Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Sat, 15 Aug 2009 23:12:39 -0400
Subject: [PATCH] in all decs, make the compute() methods return *this
 (implements feature request #18)

---
 Eigen/src/Cholesky/LDLT.h             |  7 ++++---
 Eigen/src/Cholesky/LLT.h              |  8 ++++++--
 Eigen/src/LU/LU.h                     |  7 +++++--
 Eigen/src/LU/PartialLU.h              |  5 +++--
 Eigen/src/QR/EigenSolver.h            |  5 +++--
 Eigen/src/QR/QR.h                     |  5 +++--
 Eigen/src/QR/SelfAdjointEigenSolver.h | 12 +++++++-----
 Eigen/src/SVD/JacobiSquareSVD.h       |  5 +++--
 Eigen/src/SVD/SVD.h                   |  9 ++++++---
 9 files changed, 40 insertions(+), 23 deletions(-)

diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 1cb1294ac..adca9fe2e 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -123,7 +123,7 @@ template<typename MatrixType> class LDLT
     template<typename Derived>
     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
 
-    void compute(const MatrixType& matrix);
+    LDLT& compute(const MatrixType& matrix);
 
   protected:
     /** \internal
@@ -142,7 +142,7 @@ template<typename MatrixType> class LDLT
 /** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix
   */
 template<typename MatrixType>
-void LDLT<MatrixType>::compute(const MatrixType& a)
+LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
 {
   ei_assert(a.rows()==a.cols());
   const int size = a.rows();
@@ -154,7 +154,7 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
     m_transpositions.setZero();
     m_sign = ei_real(a.coeff(0,0))>0 ? 1:-1;
     m_isInitialized = true;
-    return;
+    return *this;
   }
 
   RealScalar cutoff = 0, biggest_in_corner;
@@ -235,6 +235,7 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
   }
 
   m_isInitialized = true;
+  return *this;
 }
 
 /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index ef04b8fe4..3ce85b2d9 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -105,7 +105,7 @@ template<typename MatrixType, int _UpLo> class LLT
     template<typename Derived>
     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
 
-    void compute(const MatrixType& matrix);
+    LLT& compute(const MatrixType& matrix);
 
   protected:
     /** \internal
@@ -213,9 +213,12 @@ template<typename MatrixType> struct LLT_Traits<MatrixType,UpperTriangular>
 };
 
 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
+  *
+  *
+  * \returns a reference to *this
   */
 template<typename MatrixType, int _UpLo>
-void LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
+LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
 {
   assert(a.rows()==a.cols());
   const int size = a.rows();
@@ -223,6 +226,7 @@ void LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
   m_matrix = a;
 
   m_isInitialized = Traits::inplace_decomposition(m_matrix);
+  return *this;
 }
 
 /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h
index 733fa0cbc..f072ff906 100644
--- a/Eigen/src/LU/LU.h
+++ b/Eigen/src/LU/LU.h
@@ -111,8 +111,10 @@ template<typename MatrixType> class LU
       *
       * \param matrix the matrix of which to compute the LU decomposition.
       *               It is required to be nonzero.
+      *
+      * \returns a reference to *this
       */
-    void compute(const MatrixType& matrix);
+    LU& compute(const MatrixType& matrix);
 
     /** \returns the LU decomposition matrix: the upper-triangular part is U, the
       * unit-lower-triangular part is L (at least for square matrices; in the non-square
@@ -377,7 +379,7 @@ LU<MatrixType>::LU(const MatrixType& matrix)
 }
 
 template<typename MatrixType>
-void LU<MatrixType>::compute(const MatrixType& matrix)
+LU<MatrixType>& LU<MatrixType>::compute(const MatrixType& matrix)
 {
   m_originalMatrix = &matrix;
   m_lu = matrix;
@@ -447,6 +449,7 @@ void LU<MatrixType>::compute(const MatrixType& matrix)
     std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k)));
 
   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
+  return *this;
 }
 
 template<typename MatrixType>
diff --git a/Eigen/src/LU/PartialLU.h b/Eigen/src/LU/PartialLU.h
index f5c9b8ae9..eb5c3d34b 100644
--- a/Eigen/src/LU/PartialLU.h
+++ b/Eigen/src/LU/PartialLU.h
@@ -87,7 +87,7 @@ template<typename MatrixType> class PartialLU
       */
     PartialLU(const MatrixType& matrix);
 
-    void compute(const MatrixType& matrix);
+    PartialLU& compute(const MatrixType& matrix);
 
     /** \returns the LU decomposition matrix: the upper-triangular part is U, the
       * unit-lower-triangular part is L (at least for square matrices; in the non-square
@@ -350,7 +350,7 @@ void ei_partial_lu_inplace(MatrixType& lu, IntVector& row_transpositions, int& n
 }
 
 template<typename MatrixType>
-void PartialLU<MatrixType>::compute(const MatrixType& matrix)
+PartialLU<MatrixType>& PartialLU<MatrixType>::compute(const MatrixType& matrix)
 {
   m_lu = matrix;
   m_p.resize(matrix.rows());
@@ -369,6 +369,7 @@ void PartialLU<MatrixType>::compute(const MatrixType& matrix)
     std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k)));
 
   m_isInitialized = true;
+  return *this;
 }
 
 template<typename MatrixType>
diff --git a/Eigen/src/QR/EigenSolver.h b/Eigen/src/QR/EigenSolver.h
index dbe1aca89..bd7a76c49 100644
--- a/Eigen/src/QR/EigenSolver.h
+++ b/Eigen/src/QR/EigenSolver.h
@@ -118,7 +118,7 @@ template<typename _MatrixType> class EigenSolver
       return m_eivalues; 
     }
 
-    void compute(const MatrixType& matrix);
+    EigenSolver& compute(const MatrixType& matrix);
 
   private:
 
@@ -189,7 +189,7 @@ typename EigenSolver<MatrixType>::EigenvectorType EigenSolver<MatrixType>::eigen
 }
 
 template<typename MatrixType>
-void EigenSolver<MatrixType>::compute(const MatrixType& matrix)
+EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix)
 {
   assert(matrix.cols() == matrix.rows());
   int n = matrix.cols();
@@ -205,6 +205,7 @@ void EigenSolver<MatrixType>::compute(const MatrixType& matrix)
   hqr2(matH);
 
   m_isInitialized = true;
+  return *this;
 }
 
 // Nonsymmetric reduction to Hessenberg form.
diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h
index 98ef4da25..d37a019ff 100644
--- a/Eigen/src/QR/QR.h
+++ b/Eigen/src/QR/QR.h
@@ -101,7 +101,7 @@ template<typename MatrixType> class HouseholderQR
       */
     const MatrixType& matrixQR() const { return m_qr; }
 
-    void compute(const MatrixType& matrix);
+    HouseholderQR& compute(const MatrixType& matrix);
 
   protected:
     MatrixType m_qr;
@@ -112,7 +112,7 @@ template<typename MatrixType> class HouseholderQR
 #ifndef EIGEN_HIDE_HEAVY_CODE
 
 template<typename MatrixType>
-void HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
+HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
 {
   m_qr = matrix;
   m_hCoeffs.resize(matrix.cols());
@@ -175,6 +175,7 @@ void HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
     }
   }
   m_isInitialized = true;
+  return *this;
 }
 
 template<typename MatrixType>
diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h
index d07d28424..809a717b9 100644
--- a/Eigen/src/QR/SelfAdjointEigenSolver.h
+++ b/Eigen/src/QR/SelfAdjointEigenSolver.h
@@ -90,9 +90,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       compute(matA, matB, computeEigenvectors);
     }
 
-    void compute(const MatrixType& matrix, bool computeEigenvectors = true);
+    SelfAdjointEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
 
-    void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true);
+    SelfAdjointEigenSolver& 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
@@ -182,7 +182,7 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int st
   * \sa SelfAdjointEigenSolver(MatrixType,bool), compute(MatrixType,MatrixType,bool)
   */
 template<typename MatrixType>
-void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
+SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
 {
   #ifndef NDEBUG
   m_eigenvectorsOk = computeEigenvectors;
@@ -195,7 +195,7 @@ void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool
   {
     m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0));
     m_eivec.setOnes();
-    return;
+    return *this;
   }
 
   m_eivec = matrix;
@@ -240,6 +240,7 @@ void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool
       m_eivec.col(i).swap(m_eivec.col(k+i));
     }
   }
+  return *this;
 }
 
 /** Computes the eigenvalues of the generalized eigen problem
@@ -250,7 +251,7 @@ void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool
   * \sa SelfAdjointEigenSolver(MatrixType,MatrixType,bool), compute(MatrixType,bool)
   */
 template<typename MatrixType>
-void SelfAdjointEigenSolver<MatrixType>::
+SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::
 compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors)
 {
   ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
@@ -282,6 +283,7 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors
     for (int i=0; i<m_eivec.cols(); ++i)
       m_eivec.col(i) = m_eivec.col(i).normalized();
   }
+  return *this;
 }
 
 #endif // EIGEN_HIDE_HEAVY_CODE
diff --git a/Eigen/src/SVD/JacobiSquareSVD.h b/Eigen/src/SVD/JacobiSquareSVD.h
index 32adb1301..f21c9edca 100644
--- a/Eigen/src/SVD/JacobiSquareSVD.h
+++ b/Eigen/src/SVD/JacobiSquareSVD.h
@@ -67,7 +67,7 @@ template<typename MatrixType, bool ComputeU, bool ComputeV> class JacobiSquareSV
       compute(matrix);
     }
     
-    void compute(const MatrixType& matrix);
+    JacobiSquareSVD& compute(const MatrixType& matrix);
     
     const MatrixUType& matrixU() const
     {
@@ -95,7 +95,7 @@ template<typename MatrixType, bool ComputeU, bool ComputeV> class JacobiSquareSV
 };
 
 template<typename MatrixType, bool ComputeU, bool ComputeV>
-void JacobiSquareSVD<MatrixType, ComputeU, ComputeV>::compute(const MatrixType& matrix)
+JacobiSquareSVD<MatrixType, ComputeU, ComputeV>& JacobiSquareSVD<MatrixType, ComputeU, ComputeV>::compute(const MatrixType& matrix)
 {
   MatrixType work_matrix(matrix);
   int size = matrix.rows();
@@ -164,5 +164,6 @@ sweep_again:
   }
 
   m_isInitialized = true;
+  return *this;
 }
 #endif // EIGEN_JACOBISQUARESVD_H
diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h
index 2948087a2..1a7e6c49a 100644
--- a/Eigen/src/SVD/SVD.h
+++ b/Eigen/src/SVD/SVD.h
@@ -97,7 +97,7 @@ template<typename MatrixType> class SVD
       return m_matV;
     }
 
-    void compute(const MatrixType& matrix);
+    SVD& compute(const MatrixType& matrix);
 
     template<typename UnitaryType, typename PositiveType>
     void computeUnitaryPositive(UnitaryType *unitary, PositiveType *positive) const;
@@ -138,9 +138,11 @@ template<typename MatrixType> class SVD
 /** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix
   *
   * \note this code has been adapted from Numerical Recipes, third edition.
+  *
+  * \returns a reference to *this
   */
 template<typename MatrixType>
-void SVD<MatrixType>::compute(const MatrixType& matrix)
+SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
 {
   const int m = matrix.rows();
   const int n = matrix.cols();
@@ -157,7 +159,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
   SingularValuesType& W = m_sigma;
 
   bool flag;
-  int i,its,j,jj,k,l,nm;
+  int i,its,j,k,l,nm;
   Scalar anorm, c, f, g, h, s, scale, x, y, z;
   bool convergence = true;
   Scalar eps = precision<Scalar>();
@@ -392,6 +394,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
     m_matU = A.block(0,0,m,m);
 
   m_isInitialized = true;
+  return *this;
 }
 
 /** \returns the solution of \f$ A x = b \f$ using the current SVD decomposition of A.