in all decs, make the compute() methods return *this

(implements feature request #18)
This commit is contained in:
Benoit Jacob 2009-08-15 23:12:39 -04:00
parent 65fe5f76fd
commit ee982709d3
9 changed files with 40 additions and 23 deletions

View File

@ -123,7 +123,7 @@ template<typename MatrixType> class LDLT
template<typename Derived> template<typename Derived>
bool solveInPlace(MatrixBase<Derived> &bAndX) const; bool solveInPlace(MatrixBase<Derived> &bAndX) const;
void compute(const MatrixType& matrix); LDLT& compute(const MatrixType& matrix);
protected: protected:
/** \internal /** \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 /** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix
*/ */
template<typename MatrixType> template<typename MatrixType>
void LDLT<MatrixType>::compute(const MatrixType& a) LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
{ {
ei_assert(a.rows()==a.cols()); ei_assert(a.rows()==a.cols());
const int size = a.rows(); const int size = a.rows();
@ -154,7 +154,7 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
m_transpositions.setZero(); m_transpositions.setZero();
m_sign = ei_real(a.coeff(0,0))>0 ? 1:-1; m_sign = ei_real(a.coeff(0,0))>0 ? 1:-1;
m_isInitialized = true; m_isInitialized = true;
return; return *this;
} }
RealScalar cutoff = 0, biggest_in_corner; RealScalar cutoff = 0, biggest_in_corner;
@ -235,6 +235,7 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
} }
m_isInitialized = true; m_isInitialized = true;
return *this;
} }
/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A. /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.

View File

@ -105,7 +105,7 @@ template<typename MatrixType, int _UpLo> class LLT
template<typename Derived> template<typename Derived>
bool solveInPlace(MatrixBase<Derived> &bAndX) const; bool solveInPlace(MatrixBase<Derived> &bAndX) const;
void compute(const MatrixType& matrix); LLT& compute(const MatrixType& matrix);
protected: protected:
/** \internal /** \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 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
*
*
* \returns a reference to *this
*/ */
template<typename MatrixType, int _UpLo> 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()); assert(a.rows()==a.cols());
const int size = a.rows(); const int size = a.rows();
@ -223,6 +226,7 @@ void LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
m_matrix = a; m_matrix = a;
m_isInitialized = Traits::inplace_decomposition(m_matrix); 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. /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.

View File

@ -111,8 +111,10 @@ template<typename MatrixType> class LU
* *
* \param matrix the matrix of which to compute the LU decomposition. * \param matrix the matrix of which to compute the LU decomposition.
* It is required to be nonzero. * 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 /** \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 * 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> template<typename MatrixType>
void LU<MatrixType>::compute(const MatrixType& matrix) LU<MatrixType>& LU<MatrixType>::compute(const MatrixType& matrix)
{ {
m_originalMatrix = &matrix; m_originalMatrix = &matrix;
m_lu = 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))); std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k)));
m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_det_pq = (number_of_transpositions%2) ? -1 : 1;
return *this;
} }
template<typename MatrixType> template<typename MatrixType>

View File

@ -87,7 +87,7 @@ template<typename MatrixType> class PartialLU
*/ */
PartialLU(const MatrixType& matrix); 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 /** \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 * 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> template<typename MatrixType>
void PartialLU<MatrixType>::compute(const MatrixType& matrix) PartialLU<MatrixType>& PartialLU<MatrixType>::compute(const MatrixType& matrix)
{ {
m_lu = matrix; m_lu = matrix;
m_p.resize(matrix.rows()); 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))); std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k)));
m_isInitialized = true; m_isInitialized = true;
return *this;
} }
template<typename MatrixType> template<typename MatrixType>

View File

@ -118,7 +118,7 @@ template<typename _MatrixType> class EigenSolver
return m_eivalues; return m_eivalues;
} }
void compute(const MatrixType& matrix); EigenSolver& compute(const MatrixType& matrix);
private: private:
@ -189,7 +189,7 @@ typename EigenSolver<MatrixType>::EigenvectorType EigenSolver<MatrixType>::eigen
} }
template<typename MatrixType> template<typename MatrixType>
void EigenSolver<MatrixType>::compute(const MatrixType& matrix) EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix)
{ {
assert(matrix.cols() == matrix.rows()); assert(matrix.cols() == matrix.rows());
int n = matrix.cols(); int n = matrix.cols();
@ -205,6 +205,7 @@ void EigenSolver<MatrixType>::compute(const MatrixType& matrix)
hqr2(matH); hqr2(matH);
m_isInitialized = true; m_isInitialized = true;
return *this;
} }
// Nonsymmetric reduction to Hessenberg form. // Nonsymmetric reduction to Hessenberg form.

View File

@ -101,7 +101,7 @@ template<typename MatrixType> class HouseholderQR
*/ */
const MatrixType& matrixQR() const { return m_qr; } const MatrixType& matrixQR() const { return m_qr; }
void compute(const MatrixType& matrix); HouseholderQR& compute(const MatrixType& matrix);
protected: protected:
MatrixType m_qr; MatrixType m_qr;
@ -112,7 +112,7 @@ template<typename MatrixType> class HouseholderQR
#ifndef EIGEN_HIDE_HEAVY_CODE #ifndef EIGEN_HIDE_HEAVY_CODE
template<typename MatrixType> template<typename MatrixType>
void HouseholderQR<MatrixType>::compute(const MatrixType& matrix) HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
{ {
m_qr = matrix; m_qr = matrix;
m_hCoeffs.resize(matrix.cols()); m_hCoeffs.resize(matrix.cols());
@ -175,6 +175,7 @@ void HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
} }
} }
m_isInitialized = true; m_isInitialized = true;
return *this;
} }
template<typename MatrixType> template<typename MatrixType>

View File

@ -90,9 +90,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
compute(matA, matB, computeEigenvectors); 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 */ /** \returns the computed eigen vectors as a matrix of column vectors */
MatrixType eigenvectors(void) const 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) * \sa SelfAdjointEigenSolver(MatrixType,bool), compute(MatrixType,MatrixType,bool)
*/ */
template<typename MatrixType> template<typename MatrixType>
void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
{ {
#ifndef NDEBUG #ifndef NDEBUG
m_eigenvectorsOk = computeEigenvectors; 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_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0));
m_eivec.setOnes(); m_eivec.setOnes();
return; return *this;
} }
m_eivec = matrix; 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)); m_eivec.col(i).swap(m_eivec.col(k+i));
} }
} }
return *this;
} }
/** Computes the eigenvalues of the generalized eigen problem /** 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) * \sa SelfAdjointEigenSolver(MatrixType,MatrixType,bool), compute(MatrixType,bool)
*/ */
template<typename MatrixType> template<typename MatrixType>
void SelfAdjointEigenSolver<MatrixType>:: SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::
compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors) compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors)
{ {
ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); 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) for (int i=0; i<m_eivec.cols(); ++i)
m_eivec.col(i) = m_eivec.col(i).normalized(); m_eivec.col(i) = m_eivec.col(i).normalized();
} }
return *this;
} }
#endif // EIGEN_HIDE_HEAVY_CODE #endif // EIGEN_HIDE_HEAVY_CODE

View File

@ -67,7 +67,7 @@ template<typename MatrixType, bool ComputeU, bool ComputeV> class JacobiSquareSV
compute(matrix); compute(matrix);
} }
void compute(const MatrixType& matrix); JacobiSquareSVD& compute(const MatrixType& matrix);
const MatrixUType& matrixU() const const MatrixUType& matrixU() const
{ {
@ -95,7 +95,7 @@ template<typename MatrixType, bool ComputeU, bool ComputeV> class JacobiSquareSV
}; };
template<typename MatrixType, bool ComputeU, bool ComputeV> 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); MatrixType work_matrix(matrix);
int size = matrix.rows(); int size = matrix.rows();
@ -164,5 +164,6 @@ sweep_again:
} }
m_isInitialized = true; m_isInitialized = true;
return *this;
} }
#endif // EIGEN_JACOBISQUARESVD_H #endif // EIGEN_JACOBISQUARESVD_H

View File

@ -97,7 +97,7 @@ template<typename MatrixType> class SVD
return m_matV; return m_matV;
} }
void compute(const MatrixType& matrix); SVD& compute(const MatrixType& matrix);
template<typename UnitaryType, typename PositiveType> template<typename UnitaryType, typename PositiveType>
void computeUnitaryPositive(UnitaryType *unitary, PositiveType *positive) const; 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 /** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix
* *
* \note this code has been adapted from Numerical Recipes, third edition. * \note this code has been adapted from Numerical Recipes, third edition.
*
* \returns a reference to *this
*/ */
template<typename MatrixType> 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 m = matrix.rows();
const int n = matrix.cols(); const int n = matrix.cols();
@ -157,7 +159,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
SingularValuesType& W = m_sigma; SingularValuesType& W = m_sigma;
bool flag; 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; Scalar anorm, c, f, g, h, s, scale, x, y, z;
bool convergence = true; bool convergence = true;
Scalar eps = precision<Scalar>(); Scalar eps = precision<Scalar>();
@ -392,6 +394,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
m_matU = A.block(0,0,m,m); m_matU = A.block(0,0,m,m);
m_isInitialized = true; m_isInitialized = true;
return *this;
} }
/** \returns the solution of \f$ A x = b \f$ using the current SVD decomposition of A. /** \returns the solution of \f$ A x = b \f$ using the current SVD decomposition of A.