diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 5d5b147d9..2c25576a9 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -189,7 +189,7 @@ template class MatrixBase /** \returns the size of the inner dimension according to the storage order, * i.e., the number of rows for a columns major matrix, and the number of cols otherwise */ int innerSize() const { return (int(Flags)&RowMajorBit) ? this->cols() : this->rows(); } - + /** Only plain matrices, not expressions may be resized; therefore the only useful resize method is * Matrix::resize(). The present method only asserts that the new size equals the old size, and does * nothing else. @@ -202,7 +202,7 @@ template class MatrixBase */ void resize(int rows, int cols) { ei_assert(rows == this->rows() && cols == this->cols() && "MatrixBase::resize() does not actually allow to resize."); } - + #ifndef EIGEN_PARSED_BY_DOXYGEN /** \internal the plain matrix type corresponding to this expression. Note that is not necessarily * exactly the return type of eval(): in the case of plain matrices, the return type of eval() is a const diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h index 248b48044..4bf7ffdc9 100644 --- a/Eigen/src/LU/Inverse.h +++ b/Eigen/src/LU/Inverse.h @@ -200,7 +200,7 @@ struct ei_compute_inverse { static inline void run(const MatrixType& matrix, MatrixType* result) { - matrix.partialLu().computeInverse(result); + result = matrix.partialLu().inverse(); } }; @@ -281,9 +281,7 @@ inline void MatrixBase::computeInverse(PlainMatrixType *result) const template inline const typename MatrixBase::PlainMatrixType MatrixBase::inverse() const { - PlainMatrixType result(rows(), cols()); - computeInverse(&result); - return result; + return inverse(*this); } @@ -299,7 +297,7 @@ struct ei_compute_inverse_with_check typedef typename MatrixType::Scalar Scalar; LU lu( matrix ); if( !lu.isInvertible() ) return false; - lu.computeInverse(result); + *result = lu.inverse(); return true; } }; diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index 300c7152f..fb81713a3 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -58,7 +58,7 @@ template struct ei_lu_image_impl; * \include class_LU.cpp * Output: \verbinclude class_LU.out * - * \sa MatrixBase::lu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse() + * \sa MatrixBase::lu(), MatrixBase::determinant(), MatrixBase::inverse() */ template class LU { @@ -193,7 +193,7 @@ template class LU * Example: \include LU_solve.cpp * Output: \verbinclude LU_solve.out * - * \sa TriangularView::solve(), kernel(), inverse(), computeInverse() + * \sa TriangularView::solve(), kernel(), inverse() */ template inline const ei_lu_solve_impl @@ -277,34 +277,19 @@ template class LU return isInjective() && isSurjective(); } - /** Computes the inverse of the matrix of which *this is the LU decomposition. - * - * \param result a pointer to the matrix into which to store the inverse. Resized if needed. - * - * \note If this matrix is not invertible, *result is left with undefined coefficients. - * Use isInvertible() to first determine whether this matrix is invertible. - * - * \sa MatrixBase::computeInverse(), inverse() - */ - inline void computeInverse(MatrixType *result) const - { - ei_assert(m_originalMatrix != 0 && "LU is not initialized."); - ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!"); - *result = solve(MatrixType::Identity(m_lu.rows(), m_lu.cols())); - } - /** \returns the inverse of the matrix of which *this is the LU decomposition. * * \note If this matrix is not invertible, the returned matrix has undefined coefficients. * Use isInvertible() to first determine whether this matrix is invertible. * - * \sa computeInverse(), MatrixBase::inverse() + * \sa MatrixBase::inverse() */ - inline MatrixType inverse() const + inline const ei_lu_solve_impl > inverse() const { - MatrixType result; - computeInverse(&result); - return result; + ei_assert(m_originalMatrix != 0 && "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 > + (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue()); } protected: diff --git a/test/lu.cpp b/test/lu.cpp index a7217f5b9..442b87f82 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -40,6 +40,7 @@ template void lu_non_invertible() typename ei_lu_kernel_impl::ReturnMatrixType m1kernel = lu.kernel(); typename ei_lu_image_impl ::ReturnMatrixType m1image = lu.image(); + // std::cerr << rank << " " << lu.rank() << std::endl; VERIFY(rank == lu.rank()); VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); VERIFY(!lu.isInjective()); @@ -54,7 +55,7 @@ template void lu_non_invertible() m3 = m1*m2; m2 = MatrixType::Random(cols,cols2); // test that the code, which does resize(), may be applied to an xpr - m2.block(0,0,cols,cols2) = lu.solve(m3); + m2.block(0,0,m2.rows(),m2.cols()) = lu.solve(m3); VERIFY_IS_APPROX(m3, m1*m2); typedef Matrix SquareMatrixType; @@ -111,7 +112,6 @@ template void lu_verify_assert() VERIFY_RAISES_ASSERT(lu.isInjective()) VERIFY_RAISES_ASSERT(lu.isSurjective()) VERIFY_RAISES_ASSERT(lu.isInvertible()) - VERIFY_RAISES_ASSERT(lu.computeInverse(&tmp)) VERIFY_RAISES_ASSERT(lu.inverse()) PartialLU plu; @@ -119,7 +119,6 @@ template void lu_verify_assert() VERIFY_RAISES_ASSERT(plu.permutationP()) VERIFY_RAISES_ASSERT(plu.solve(tmp,&tmp)) VERIFY_RAISES_ASSERT(plu.determinant()) - VERIFY_RAISES_ASSERT(plu.computeInverse(&tmp)) VERIFY_RAISES_ASSERT(plu.inverse()) }