diff options
Diffstat (limited to 'Eigen/src/LU')
-rw-r--r-- | Eigen/src/LU/LU.h | 63 |
1 files changed, 26 insertions, 37 deletions
diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index b36dc3026..00a6dcf64 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -45,14 +45,9 @@ * The data of the LU decomposition can be directly accessed through the methods matrixLU(), * permutationP(), permutationQ(). Convenience methods matrixL(), matrixU() are also provided. * - * As an exemple, here is how the original matrix can be retrieved, in the square case: - * \include class_LU_1.cpp - * Output: \verbinclude class_LU_1.out - * - * When the matrix is not square, matrixL() is no longer very useful: if one needs it, one has - * to construct the L matrix by hand, as shown in this example: - * \include class_LU_2.cpp - * Output: \verbinclude class_LU_2.out + * As an exemple, here is how the original matrix can be retrieved: + * \include class_LU.cpp + * Output: \verbinclude class_LU.out * * \sa MatrixBase::lu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse() */ @@ -91,6 +86,14 @@ template<typename MatrixType> class LU MatrixType::MaxColsAtCompileTime // so it has the same number of rows and at most as many columns. > ImageResultType; + typedef Matrix<typename MatrixType::Scalar, + MatrixType::RowsAtCompileTime, + MatrixType::RowsAtCompileTime, + MatrixType::Options, + MatrixType::MaxRowsAtCompileTime, + MatrixType::MaxRowsAtCompileTime + > MatrixLType; + /** Constructor. * * \param matrix the matrix of which to compute the LU decomposition. @@ -108,26 +111,6 @@ template<typename MatrixType> class LU return m_lu; } - /** \returns an expression of the unit-lower-triangular part of the LU matrix. In the square case, - * this is the L matrix. In the non-square, actually obtaining the L matrix takes some - * more care, see the documentation of class LU. - * - * \sa matrixLU(), matrixU() - */ - inline const Part<MatrixType, UnitLowerTriangular> matrixL() const - { - return m_lu; - } - - /** \returns an expression of the U matrix, i.e. the upper-triangular part of the LU matrix. - * - * \sa matrixLU(), matrixL() - */ - inline const Part<MatrixType, UpperTriangular> matrixU() const - { - return m_lu; - } - /** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed, * representing the P permutation i.e. the permutation of the rows. For its precise meaning, * see the examples given in the documentation of class LU. @@ -495,7 +478,7 @@ bool LU<MatrixType>::solve( * Step 4: result = Qd; */ - const int rows = m_lu.rows(); + const int rows = m_lu.rows(), cols = m_lu.cols(); ei_assert(b.rows() == rows); const int smalldim = std::min(rows, m_lu.cols()); @@ -505,14 +488,20 @@ bool LU<MatrixType>::solve( for(int i = 0; i < rows; ++i) c.row(m_p.coeff(i)) = b.row(i); // Step 2 - Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime, - MatrixType::Options, - MatrixType::MaxRowsAtCompileTime, - MatrixType::MaxRowsAtCompileTime> l(rows, rows); - l.setZero(); - l.corner(Eigen::TopLeft,rows,smalldim) - = m_lu.corner(Eigen::TopLeft,rows,smalldim); - l.template marked<UnitLowerTriangular>().solveTriangularInPlace(c); + if(rows <= cols) + m_lu.corner(Eigen::TopLeft,rows,smalldim).template marked<UnitLowerTriangular>().solveTriangularInPlace(c); + else + { + // construct the L matrix. We shouldn't do that everytime, it is a very large overhead in the case of vector solving. + // However the case rows>cols is rather unusual with LU so this is probably not a huge priority. + Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime, + MatrixType::Options, + MatrixType::MaxRowsAtCompileTime, + MatrixType::MaxRowsAtCompileTime> l(rows, rows); + l.setZero(); + l.corner(Eigen::TopLeft,rows,smalldim) = m_lu.corner(Eigen::TopLeft,rows,smalldim); + l.template marked<UnitLowerTriangular>().solveTriangularInPlace(c); + } // Step 3 if(!isSurjective()) |