diff options
| -rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 88 | ||||
| -rw-r--r-- | Eigen/src/Cholesky/LLT.h | 86 | ||||
| -rw-r--r-- | Eigen/src/LU/FullPivLU.h | 40 | ||||
| -rw-r--r-- | Eigen/src/LU/PartialPivLU.h | 26 | ||||
| -rw-r--r-- | doc/Doxyfile.in | 2 | ||||
| -rw-r--r-- | doc/snippets/LLT_solve.cpp | 4 | ||||
| -rw-r--r-- | test/cholesky.cpp | 16 | ||||
| -rw-r--r-- | test/lu.cpp | 4 |
8 files changed, 169 insertions, 97 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index f95e10935..04c02d4fc 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -27,6 +27,8 @@ #ifndef EIGEN_LDLT_H #define EIGEN_LDLT_H +template<typename MatrixType, typename Rhs> struct ei_ldlt_solve_impl; + /** \ingroup cholesky_Module * * \class LDLT @@ -43,8 +45,8 @@ * zeros in the bottom right rank(A) - n submatrix. Avoiding the square root * on D also stabilizes the computation. * - * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky decomposition to determine - * whether a system of equations has a solution. + * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky + * decomposition to determine whether a system of equations has a solution. * * \sa MatrixBase::ldlt(), class LLT */ @@ -117,14 +119,37 @@ template<typename MatrixType> class LDLT return m_sign == -1; } - template<typename RhsDerived, typename ResultType> - bool solve(const MatrixBase<RhsDerived> &b, ResultType *result) const; - + /** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A. + * + * \note_about_checking_solutions + * + * \sa solveInPlace(), MatrixBase::ldlt() + */ + template<typename Rhs> + inline const ei_ldlt_solve_impl<MatrixType, Rhs> + solve(const MatrixBase<Rhs>& b) const + { + ei_assert(m_isInitialized && "LDLT is not initialized."); + ei_assert(m_matrix.rows()==b.rows() + && "LDLT::solve(): invalid number of rows of the right hand side matrix b"); + return ei_ldlt_solve_impl<MatrixType, Rhs>(*this, b.derived()); + } + template<typename Derived> bool solveInPlace(MatrixBase<Derived> &bAndX) const; LDLT& compute(const MatrixType& matrix); + /** \returns the LDLT decomposition matrix + * + * TODO: document the storage layout + */ + inline const MatrixType& matrixLDLT() const + { + ei_assert(m_isInitialized && "LDLT is not initialized."); + return m_matrix; + } + protected: /** \internal * Used to compute and store the Cholesky decomposition A = L D L^* = U^* D U. @@ -134,7 +159,7 @@ template<typename MatrixType> class LDLT */ MatrixType m_matrix; IntColVectorType m_p; - IntColVectorType m_transpositions; + IntColVectorType m_transpositions; // FIXME do we really need to store permanently the transpositions? int m_sign; bool m_isInitialized; }; @@ -238,27 +263,38 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a) return *this; } -/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A. - * The result is stored in \a result - * - * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD. - * - * In other words, it computes \f$ b = A^{-1} b \f$ with - * \f$ P^T{L^{*}}^{-1} D^{-1} L^{-1} P b \f$ from right to left. - * - * \sa LDLT::solveInPlace(), MatrixBase::ldlt() - */ -template<typename MatrixType> -template<typename RhsDerived, typename ResultType> -bool LDLT<MatrixType> -::solve(const MatrixBase<RhsDerived> &b, ResultType *result) const +template<typename MatrixType,typename Rhs> +struct ei_traits<ei_ldlt_solve_impl<MatrixType,Rhs> > { - ei_assert(m_isInitialized && "LDLT is not initialized."); - const int size = m_matrix.rows(); - ei_assert(size==b.rows() && "LDLT::solve(): invalid number of rows of the right hand side matrix b"); - *result = b; - return solveInPlace(*result); -} + typedef Matrix<typename Rhs::Scalar, + MatrixType::ColsAtCompileTime, + Rhs::ColsAtCompileTime, + Rhs::PlainMatrixType::Options, + MatrixType::MaxColsAtCompileTime, + Rhs::MaxColsAtCompileTime> ReturnMatrixType; +}; + +template<typename MatrixType, typename Rhs> +struct ei_ldlt_solve_impl : public ReturnByValue<ei_ldlt_solve_impl<MatrixType, Rhs> > +{ + typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested; + typedef LDLT<MatrixType> LDLTType; + const LDLTType& m_ldlt; + const typename Rhs::Nested m_rhs; + + ei_ldlt_solve_impl(const LDLTType& ldlt, const Rhs& rhs) + : m_ldlt(ldlt), m_rhs(rhs) + {} + + inline int rows() const { return m_ldlt.matrixLDLT().cols(); } + inline int cols() const { return m_rhs.cols(); } + + template<typename Dest> void evalTo(Dest& dst) const + { + dst = m_rhs; + m_ldlt.solveInPlace(dst); + } +}; /** This is the \em in-place version of solve(). * diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index ec7b8123c..d6fd514cc 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -26,6 +26,7 @@ #define EIGEN_LLT_H template<typename MatrixType, int UpLo> struct LLT_Traits; +template<typename MatrixType, int UpLo, typename Rhs> struct ei_llt_solve_impl; /** \ingroup cholesky_Module * @@ -99,14 +100,41 @@ template<typename MatrixType, int _UpLo> class LLT return Traits::getL(m_matrix); } - template<typename RhsDerived, typename ResultType> - bool solve(const MatrixBase<RhsDerived> &b, ResultType *result) const; - + /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. + * + * Since this LLT class assumes anyway that the matrix A is invertible, the solution + * theoretically exists and is unique regardless of b. + * + * Example: \include LLT_solve.cpp + * Output: \verbinclude LLT_solve.out + * + * \sa solveInPlace(), MatrixBase::llt() + */ + template<typename Rhs> + inline const ei_llt_solve_impl<MatrixType, UpLo, Rhs> + solve(const MatrixBase<Rhs>& b) const + { + ei_assert(m_isInitialized && "LLT is not initialized."); + ei_assert(m_matrix.rows()==b.rows() + && "LLT::solve(): invalid number of rows of the right hand side matrix b"); + return ei_llt_solve_impl<MatrixType, UpLo, Rhs>(*this, b.derived()); + } + template<typename Derived> bool solveInPlace(MatrixBase<Derived> &bAndX) const; LLT& compute(const MatrixType& matrix); + /** \returns the LLT decomposition matrix + * + * TODO: document the storage layout + */ + inline const MatrixType& matrixLLT() const + { + ei_assert(m_isInitialized && "LLT is not initialized."); + return m_matrix; + } + protected: /** \internal * Used to compute and store L @@ -229,28 +257,38 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a) return *this; } -/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A. - * The result is stored in \a result - * - * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD. - * - * In other words, it computes \f$ b = A^{-1} b \f$ with - * \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left. - * - * Example: \include LLT_solve.cpp - * Output: \verbinclude LLT_solve.out - * - * \sa LLT::solveInPlace(), MatrixBase::llt() - */ -template<typename MatrixType, int _UpLo> -template<typename RhsDerived, typename ResultType> -bool LLT<MatrixType,_UpLo>::solve(const MatrixBase<RhsDerived> &b, ResultType *result) const +template<typename MatrixType, int UpLo, typename Rhs> +struct ei_traits<ei_llt_solve_impl<MatrixType,UpLo,Rhs> > { - ei_assert(m_isInitialized && "LLT is not initialized."); - const int size = m_matrix.rows(); - ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b"); - return solveInPlace((*result) = b); -} + typedef Matrix<typename Rhs::Scalar, + MatrixType::ColsAtCompileTime, + Rhs::ColsAtCompileTime, + Rhs::PlainMatrixType::Options, + MatrixType::MaxColsAtCompileTime, + Rhs::MaxColsAtCompileTime> ReturnMatrixType; +}; + +template<typename MatrixType, int UpLo, typename Rhs> +struct ei_llt_solve_impl : public ReturnByValue<ei_llt_solve_impl<MatrixType,UpLo,Rhs> > +{ + typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested; + typedef LLT<MatrixType,UpLo> LLTType; + const LLTType& m_llt; + const typename Rhs::Nested m_rhs; + + ei_llt_solve_impl(const LLTType& llt, const Rhs& rhs) + : m_llt(llt), m_rhs(rhs) + {} + + inline int rows() const { return m_llt.matrixLLT().cols(); } + inline int cols() const { return m_rhs.cols(); } + + template<typename Dest> void evalTo(Dest& dst) const + { + dst = m_rhs; + m_llt.solveInPlace(dst); + } +}; /** This is the \em in-place version of solve(). * diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 8743dac92..a28a536b6 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -25,9 +25,9 @@ #ifndef EIGEN_LU_H #define EIGEN_LU_H -template<typename MatrixType, typename Rhs> struct ei_lu_solve_impl; -template<typename MatrixType> struct ei_lu_kernel_impl; -template<typename MatrixType> struct ei_lu_image_impl; +template<typename MatrixType, typename Rhs> struct ei_fullpivlu_solve_impl; +template<typename MatrixType> struct ei_fullpivlu_kernel_impl; +template<typename MatrixType> struct ei_fullpivlu_image_impl; /** \ingroup LU_Module * @@ -167,10 +167,10 @@ template<typename MatrixType> class FullPivLU * * \sa image() */ - inline const ei_lu_kernel_impl<MatrixType> kernel() const + inline const ei_fullpivlu_kernel_impl<MatrixType> kernel() const { ei_assert(m_isInitialized && "LU is not initialized."); - return ei_lu_kernel_impl<MatrixType>(*this); + return ei_fullpivlu_kernel_impl<MatrixType>(*this); } /** \returns the image of the matrix, also called its column-space. The columns of the returned matrix @@ -193,11 +193,11 @@ template<typename MatrixType> class FullPivLU * \sa kernel() */ template<typename OriginalMatrixType> - inline const ei_lu_image_impl<MatrixType> + inline const ei_fullpivlu_image_impl<MatrixType> image(const MatrixBase<OriginalMatrixType>& originalMatrix) const { ei_assert(m_isInitialized && "LU is not initialized."); - return ei_lu_image_impl<MatrixType>(*this, originalMatrix.derived()); + return ei_fullpivlu_image_impl<MatrixType>(*this, originalMatrix.derived()); } /** This method returns a solution x to the equation Ax=b, where A is the matrix of which @@ -220,11 +220,11 @@ template<typename MatrixType> class FullPivLU * \sa TriangularView::solve(), kernel(), inverse() */ template<typename Rhs> - inline const ei_lu_solve_impl<MatrixType, Rhs> + inline const ei_fullpivlu_solve_impl<MatrixType, Rhs> solve(const MatrixBase<Rhs>& b) const { ei_assert(m_isInitialized && "LU is not initialized."); - return ei_lu_solve_impl<MatrixType, Rhs>(*this, b.derived()); + return ei_fullpivlu_solve_impl<MatrixType, Rhs>(*this, b.derived()); } /** \returns the determinant of the matrix of which @@ -365,11 +365,11 @@ template<typename MatrixType> class FullPivLU * * \sa MatrixBase::inverse() */ - inline const ei_lu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const + inline const ei_fullpivlu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const { ei_assert(m_isInitialized && "LU is not initialized."); ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!"); - return ei_lu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > + return ei_fullpivlu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue()); } @@ -493,7 +493,7 @@ typename ei_traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() cons /********* Implementation of kernel() **************************************************/ template<typename MatrixType> -struct ei_traits<ei_lu_kernel_impl<MatrixType> > +struct ei_traits<ei_fullpivlu_kernel_impl<MatrixType> > { typedef Matrix< typename MatrixType::Scalar, @@ -509,7 +509,7 @@ struct ei_traits<ei_lu_kernel_impl<MatrixType> > }; template<typename MatrixType> -struct ei_lu_kernel_impl : public ReturnByValue<ei_lu_kernel_impl<MatrixType> > +struct ei_fullpivlu_kernel_impl : public ReturnByValue<ei_fullpivlu_kernel_impl<MatrixType> > { typedef FullPivLU<MatrixType> LUType; typedef typename MatrixType::Scalar Scalar; @@ -517,7 +517,7 @@ struct ei_lu_kernel_impl : public ReturnByValue<ei_lu_kernel_impl<MatrixType> > const LUType& m_lu; int m_rank, m_cols; - ei_lu_kernel_impl(const LUType& lu) + ei_fullpivlu_kernel_impl(const LUType& lu) : m_lu(lu), m_rank(lu.rank()), m_cols(m_rank==lu.matrixLU().cols() ? 1 : lu.matrixLU().cols() - m_rank){} @@ -599,7 +599,7 @@ struct ei_lu_kernel_impl : public ReturnByValue<ei_lu_kernel_impl<MatrixType> > /***** Implementation of image() *****************************************************/ template<typename MatrixType> -struct ei_traits<ei_lu_image_impl<MatrixType> > +struct ei_traits<ei_fullpivlu_image_impl<MatrixType> > { typedef Matrix< typename MatrixType::Scalar, @@ -613,7 +613,7 @@ struct ei_traits<ei_lu_image_impl<MatrixType> > }; template<typename MatrixType> -struct ei_lu_image_impl : public ReturnByValue<ei_lu_image_impl<MatrixType> > +struct ei_fullpivlu_image_impl : public ReturnByValue<ei_fullpivlu_image_impl<MatrixType> > { typedef FullPivLU<MatrixType> LUType; typedef typename MatrixType::RealScalar RealScalar; @@ -621,7 +621,7 @@ struct ei_lu_image_impl : public ReturnByValue<ei_lu_image_impl<MatrixType> > int m_rank, m_cols; const MatrixType& m_originalMatrix; - ei_lu_image_impl(const LUType& lu, const MatrixType& originalMatrix) + ei_fullpivlu_image_impl(const LUType& lu, const MatrixType& originalMatrix) : m_lu(lu), m_rank(lu.rank()), m_cols(m_rank == 0 ? 1 : m_rank), m_originalMatrix(originalMatrix) {} @@ -656,7 +656,7 @@ struct ei_lu_image_impl : public ReturnByValue<ei_lu_image_impl<MatrixType> > /***** Implementation of solve() *****************************************************/ template<typename MatrixType,typename Rhs> -struct ei_traits<ei_lu_solve_impl<MatrixType,Rhs> > +struct ei_traits<ei_fullpivlu_solve_impl<MatrixType,Rhs> > { typedef Matrix<typename Rhs::Scalar, MatrixType::ColsAtCompileTime, @@ -667,14 +667,14 @@ struct ei_traits<ei_lu_solve_impl<MatrixType,Rhs> > }; template<typename MatrixType, typename Rhs> -struct ei_lu_solve_impl : public ReturnByValue<ei_lu_solve_impl<MatrixType, Rhs> > +struct ei_fullpivlu_solve_impl : public ReturnByValue<ei_fullpivlu_solve_impl<MatrixType, Rhs> > { typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested; typedef FullPivLU<MatrixType> LUType; const LUType& m_lu; const typename Rhs::Nested m_rhs; - ei_lu_solve_impl(const LUType& lu, const Rhs& rhs) + ei_fullpivlu_solve_impl(const LUType& lu, const Rhs& rhs) : m_lu(lu), m_rhs(rhs) {} diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 647ada38f..2a04ec4a9 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -26,7 +26,7 @@ #ifndef EIGEN_PARTIALLU_H #define EIGEN_PARTIALLU_H -template<typename MatrixType, typename Rhs> struct ei_partiallu_solve_impl; +template<typename MatrixType, typename Rhs> struct ei_partialpivlu_solve_impl; /** \ingroup LU_Module * @@ -116,7 +116,7 @@ template<typename MatrixType> class PartialPivLU return m_p; } - /** This method returns a solution x to the equation Ax=b, where A is the matrix of which + /** This method returns the solution x to the equation Ax=b, where A is the matrix of which * *this is the LU decomposition. * * \param b the right-hand-side of the equation to solve. Can be a vector or a matrix, @@ -131,16 +131,14 @@ template<typename MatrixType> class PartialPivLU * Since this PartialPivLU class assumes anyway that the matrix A is invertible, the solution * theoretically exists and is unique regardless of b. * - * \note_about_checking_solutions - * * \sa TriangularView::solve(), inverse(), computeInverse() */ template<typename Rhs> - inline const ei_partiallu_solve_impl<MatrixType, Rhs> + inline const ei_partialpivlu_solve_impl<MatrixType, Rhs> solve(const MatrixBase<Rhs>& b) const { - ei_assert(m_isInitialized && "LU is not initialized."); - return ei_partiallu_solve_impl<MatrixType, Rhs>(*this, b.derived()); + ei_assert(m_isInitialized && "PartialPivLU is not initialized."); + return ei_partialpivlu_solve_impl<MatrixType, Rhs>(*this, b.derived()); } /** \returns the inverse of the matrix of which *this is the LU decomposition. @@ -150,10 +148,10 @@ template<typename MatrixType> class PartialPivLU * * \sa MatrixBase::inverse(), LU::inverse() */ - inline const ei_partiallu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const + inline const ei_partialpivlu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const { - ei_assert(m_isInitialized && "LU is not initialized."); - return ei_partiallu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > + ei_assert(m_isInitialized && "PartialPivLU is not initialized."); + return ei_partialpivlu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue()); } @@ -200,7 +198,7 @@ PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix) -/** This is the blocked version of ei_lu_unblocked() */ +/** This is the blocked version of ei_fullpivlu_unblocked() */ template<typename Scalar, int StorageOrder> struct ei_partial_lu_impl { @@ -410,7 +408,7 @@ typename ei_traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() c /***** Implementation of solve() *****************************************************/ template<typename MatrixType,typename Rhs> -struct ei_traits<ei_partiallu_solve_impl<MatrixType,Rhs> > +struct ei_traits<ei_partialpivlu_solve_impl<MatrixType,Rhs> > { typedef Matrix<typename Rhs::Scalar, MatrixType::ColsAtCompileTime, @@ -421,14 +419,14 @@ struct ei_traits<ei_partiallu_solve_impl<MatrixType,Rhs> > }; template<typename MatrixType, typename Rhs> -struct ei_partiallu_solve_impl : public ReturnByValue<ei_partiallu_solve_impl<MatrixType, Rhs> > +struct ei_partialpivlu_solve_impl : public ReturnByValue<ei_partialpivlu_solve_impl<MatrixType, Rhs> > { typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested; typedef PartialPivLU<MatrixType> LUType; const LUType& m_lu; const typename Rhs::Nested m_rhs; - ei_partiallu_solve_impl(const LUType& lu, const Rhs& rhs) + ei_partialpivlu_solve_impl(const LUType& lu, const Rhs& rhs) : m_lu(lu), m_rhs(rhs) {} diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in index fb00a3cc5..412831e14 100644 --- a/doc/Doxyfile.in +++ b/doc/Doxyfile.in @@ -218,7 +218,7 @@ ALIASES = "only_for_vectors=This is only for vectors (either row- "nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\"" \ "note_about_arbitrary_choice_of_solution=If there exists more than one solution, this method will arbitrarily choose one." \ "note_about_using_kernel_to_study_multiple_solutions=If you need a complete analysis of the space of solutions, take the one solution obtained by this method and add to it elements of the kernel, as determined by kernel()." \ - "note_about_checking_solutions=This method just tries to find as good a solution as possible. If you want to check whether a solution "exists" or if it is accurate, just call this function to get a solution and then compute the error margin, or use MatrixBase::isApprox() directly, for instance like this: \code bool a_solution_exists = (A*result).isApprox(b, precision); \endcode The non-existence of a solution doesn't by itself mean that you'll get \c inf or \c nan values." + "note_about_checking_solutions=This method just tries to find as good a solution as possible. If you want to check whether a solution exists or if it is accurate, just call this function to get a result and then compute the error of this result, or use MatrixBase::isApprox() directly, for instance like this: \code bool a_solution_exists = (A*result).isApprox(b, precision); \endcode This method avoids dividing by zero, so that the non-existence of a solution doesn't by itself mean that you'll get \c inf or \c nan values." # Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C # sources only. Doxygen will then generate output that is more tailored for C. diff --git a/doc/snippets/LLT_solve.cpp b/doc/snippets/LLT_solve.cpp index 76ab09ec5..7095d2cc3 100644 --- a/doc/snippets/LLT_solve.cpp +++ b/doc/snippets/LLT_solve.cpp @@ -3,6 +3,6 @@ typedef Matrix<float,Dynamic,2> DataMatrix; DataMatrix samples = DataMatrix::Random(12,2); VectorXf elevations = 2*samples.col(0) + 3*samples.col(1) + VectorXf::Random(12)*0.1; // and let's solve samples * [x y]^T = elevations in least square sense: -Matrix<float,2,1> xy; -(samples.adjoint() * samples).llt().solve((samples.adjoint()*elevations), &xy); +Matrix<float,2,1> xy + = (samples.adjoint() * samples).llt().solve((samples.adjoint()*elevations)); cout << xy << endl; diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 77d025e59..83f7e0492 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -93,17 +93,17 @@ template<typename MatrixType> void cholesky(const MatrixType& m) { LLT<SquareMatrixType,LowerTriangular> chollo(symmLo); VERIFY_IS_APPROX(symm, chollo.matrixL().toDense() * chollo.matrixL().adjoint().toDense()); - chollo.solve(vecB, &vecX); + vecX = chollo.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); - chollo.solve(matB, &matX); + matX = chollo.solve(matB); VERIFY_IS_APPROX(symm * matX, matB); // test the upper mode LLT<SquareMatrixType,UpperTriangular> cholup(symmUp); VERIFY_IS_APPROX(symm, cholup.matrixL().toDense() * cholup.matrixL().adjoint().toDense()); - cholup.solve(vecB, &vecX); + vecX = cholup.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); - cholup.solve(matB, &matX); + matX = cholup.solve(matB); VERIFY_IS_APPROX(symm * matX, matB); } @@ -118,9 +118,9 @@ template<typename MatrixType> void cholesky(const MatrixType& m) LDLT<SquareMatrixType> ldlt(symm); // TODO(keir): This doesn't make sense now that LDLT pivots. //VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * ldlt.matrixL().adjoint()); - ldlt.solve(vecB, &vecX); + vecX = ldlt.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); - ldlt.solve(matB, &matX); + matX = ldlt.solve(matB); VERIFY_IS_APPROX(symm * matX, matB); } @@ -132,7 +132,7 @@ template<typename MatrixType> void cholesky_verify_assert() LLT<MatrixType> llt; VERIFY_RAISES_ASSERT(llt.matrixL()) - VERIFY_RAISES_ASSERT(llt.solve(tmp,&tmp)) + VERIFY_RAISES_ASSERT(llt.solve(tmp)) VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp)) LDLT<MatrixType> ldlt; @@ -141,7 +141,7 @@ template<typename MatrixType> void cholesky_verify_assert() VERIFY_RAISES_ASSERT(ldlt.vectorD()) VERIFY_RAISES_ASSERT(ldlt.isPositive()) VERIFY_RAISES_ASSERT(ldlt.isNegative()) - VERIFY_RAISES_ASSERT(ldlt.solve(tmp,&tmp)) + VERIFY_RAISES_ASSERT(ldlt.solve(tmp)) VERIFY_RAISES_ASSERT(ldlt.solveInPlace(&tmp)) } diff --git a/test/lu.cpp b/test/lu.cpp index b7214ae12..c46ca9130 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -49,8 +49,8 @@ template<typename MatrixType> void lu_non_invertible() cols2 = cols = MatrixType::ColsAtCompileTime; } - typedef typename ei_lu_kernel_impl<MatrixType>::ReturnMatrixType KernelMatrixType; - typedef typename ei_lu_image_impl <MatrixType>::ReturnMatrixType ImageMatrixType; + typedef typename ei_fullpivlu_kernel_impl<MatrixType>::ReturnMatrixType KernelMatrixType; + typedef typename ei_fullpivlu_image_impl <MatrixType>::ReturnMatrixType ImageMatrixType; typedef Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> DynamicMatrixType; typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> CMatrixType; |