[136] | 1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
| 2 | // for linear algebra.
|
---|
| 3 | //
|
---|
| 4 | // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
|
---|
| 5 | //
|
---|
| 6 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
| 7 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
| 8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
| 9 |
|
---|
| 10 | #ifndef EIGEN_LU_H
|
---|
| 11 | #define EIGEN_LU_H
|
---|
| 12 |
|
---|
| 13 | namespace Eigen {
|
---|
| 14 |
|
---|
| 15 | /** \ingroup LU_Module
|
---|
| 16 | *
|
---|
| 17 | * \class FullPivLU
|
---|
| 18 | *
|
---|
| 19 | * \brief LU decomposition of a matrix with complete pivoting, and related features
|
---|
| 20 | *
|
---|
| 21 | * \param MatrixType the type of the matrix of which we are computing the LU decomposition
|
---|
| 22 | *
|
---|
| 23 | * This class represents a LU decomposition of any matrix, with complete pivoting: the matrix A is
|
---|
| 24 | * decomposed as \f$ A = P^{-1} L U Q^{-1} \f$ where L is unit-lower-triangular, U is
|
---|
| 25 | * upper-triangular, and P and Q are permutation matrices. This is a rank-revealing LU
|
---|
| 26 | * decomposition. The eigenvalues (diagonal coefficients) of U are sorted in such a way that any
|
---|
| 27 | * zeros are at the end.
|
---|
| 28 | *
|
---|
| 29 | * This decomposition provides the generic approach to solving systems of linear equations, computing
|
---|
| 30 | * the rank, invertibility, inverse, kernel, and determinant.
|
---|
| 31 | *
|
---|
| 32 | * This LU decomposition is very stable and well tested with large matrices. However there are use cases where the SVD
|
---|
| 33 | * decomposition is inherently more stable and/or flexible. For example, when computing the kernel of a matrix,
|
---|
| 34 | * working with the SVD allows to select the smallest singular values of the matrix, something that
|
---|
| 35 | * the LU decomposition doesn't see.
|
---|
| 36 | *
|
---|
| 37 | * The data of the LU decomposition can be directly accessed through the methods matrixLU(),
|
---|
| 38 | * permutationP(), permutationQ().
|
---|
| 39 | *
|
---|
| 40 | * As an exemple, here is how the original matrix can be retrieved:
|
---|
| 41 | * \include class_FullPivLU.cpp
|
---|
| 42 | * Output: \verbinclude class_FullPivLU.out
|
---|
| 43 | *
|
---|
| 44 | * \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse()
|
---|
| 45 | */
|
---|
| 46 | template<typename _MatrixType> class FullPivLU
|
---|
| 47 | {
|
---|
| 48 | public:
|
---|
| 49 | typedef _MatrixType MatrixType;
|
---|
| 50 | enum {
|
---|
| 51 | RowsAtCompileTime = MatrixType::RowsAtCompileTime,
|
---|
| 52 | ColsAtCompileTime = MatrixType::ColsAtCompileTime,
|
---|
| 53 | Options = MatrixType::Options,
|
---|
| 54 | MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
|
---|
| 55 | MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
|
---|
| 56 | };
|
---|
| 57 | typedef typename MatrixType::Scalar Scalar;
|
---|
| 58 | typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
|
---|
| 59 | typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
|
---|
| 60 | typedef typename MatrixType::Index Index;
|
---|
| 61 | typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
|
---|
| 62 | typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
|
---|
| 63 | typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
|
---|
| 64 | typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
|
---|
| 65 |
|
---|
| 66 | /**
|
---|
| 67 | * \brief Default Constructor.
|
---|
| 68 | *
|
---|
| 69 | * The default constructor is useful in cases in which the user intends to
|
---|
| 70 | * perform decompositions via LU::compute(const MatrixType&).
|
---|
| 71 | */
|
---|
| 72 | FullPivLU();
|
---|
| 73 |
|
---|
| 74 | /** \brief Default Constructor with memory preallocation
|
---|
| 75 | *
|
---|
| 76 | * Like the default constructor but with preallocation of the internal data
|
---|
| 77 | * according to the specified problem \a size.
|
---|
| 78 | * \sa FullPivLU()
|
---|
| 79 | */
|
---|
| 80 | FullPivLU(Index rows, Index cols);
|
---|
| 81 |
|
---|
| 82 | /** Constructor.
|
---|
| 83 | *
|
---|
| 84 | * \param matrix the matrix of which to compute the LU decomposition.
|
---|
| 85 | * It is required to be nonzero.
|
---|
| 86 | */
|
---|
| 87 | FullPivLU(const MatrixType& matrix);
|
---|
| 88 |
|
---|
| 89 | /** Computes the LU decomposition of the given matrix.
|
---|
| 90 | *
|
---|
| 91 | * \param matrix the matrix of which to compute the LU decomposition.
|
---|
| 92 | * It is required to be nonzero.
|
---|
| 93 | *
|
---|
| 94 | * \returns a reference to *this
|
---|
| 95 | */
|
---|
| 96 | FullPivLU& compute(const MatrixType& matrix);
|
---|
| 97 |
|
---|
| 98 | /** \returns the LU decomposition matrix: the upper-triangular part is U, the
|
---|
| 99 | * unit-lower-triangular part is L (at least for square matrices; in the non-square
|
---|
| 100 | * case, special care is needed, see the documentation of class FullPivLU).
|
---|
| 101 | *
|
---|
| 102 | * \sa matrixL(), matrixU()
|
---|
| 103 | */
|
---|
| 104 | inline const MatrixType& matrixLU() const
|
---|
| 105 | {
|
---|
| 106 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 107 | return m_lu;
|
---|
| 108 | }
|
---|
| 109 |
|
---|
| 110 | /** \returns the number of nonzero pivots in the LU decomposition.
|
---|
| 111 | * Here nonzero is meant in the exact sense, not in a fuzzy sense.
|
---|
| 112 | * So that notion isn't really intrinsically interesting, but it is
|
---|
| 113 | * still useful when implementing algorithms.
|
---|
| 114 | *
|
---|
| 115 | * \sa rank()
|
---|
| 116 | */
|
---|
| 117 | inline Index nonzeroPivots() const
|
---|
| 118 | {
|
---|
| 119 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 120 | return m_nonzero_pivots;
|
---|
| 121 | }
|
---|
| 122 |
|
---|
| 123 | /** \returns the absolute value of the biggest pivot, i.e. the biggest
|
---|
| 124 | * diagonal coefficient of U.
|
---|
| 125 | */
|
---|
| 126 | RealScalar maxPivot() const { return m_maxpivot; }
|
---|
| 127 |
|
---|
| 128 | /** \returns the permutation matrix P
|
---|
| 129 | *
|
---|
| 130 | * \sa permutationQ()
|
---|
| 131 | */
|
---|
| 132 | inline const PermutationPType& permutationP() const
|
---|
| 133 | {
|
---|
| 134 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 135 | return m_p;
|
---|
| 136 | }
|
---|
| 137 |
|
---|
| 138 | /** \returns the permutation matrix Q
|
---|
| 139 | *
|
---|
| 140 | * \sa permutationP()
|
---|
| 141 | */
|
---|
| 142 | inline const PermutationQType& permutationQ() const
|
---|
| 143 | {
|
---|
| 144 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 145 | return m_q;
|
---|
| 146 | }
|
---|
| 147 |
|
---|
| 148 | /** \returns the kernel of the matrix, also called its null-space. The columns of the returned matrix
|
---|
| 149 | * will form a basis of the kernel.
|
---|
| 150 | *
|
---|
| 151 | * \note If the kernel has dimension zero, then the returned matrix is a column-vector filled with zeros.
|
---|
| 152 | *
|
---|
| 153 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
| 154 | * For that, it uses the threshold value that you can control by calling
|
---|
| 155 | * setThreshold(const RealScalar&).
|
---|
| 156 | *
|
---|
| 157 | * Example: \include FullPivLU_kernel.cpp
|
---|
| 158 | * Output: \verbinclude FullPivLU_kernel.out
|
---|
| 159 | *
|
---|
| 160 | * \sa image()
|
---|
| 161 | */
|
---|
| 162 | inline const internal::kernel_retval<FullPivLU> kernel() const
|
---|
| 163 | {
|
---|
| 164 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 165 | return internal::kernel_retval<FullPivLU>(*this);
|
---|
| 166 | }
|
---|
| 167 |
|
---|
| 168 | /** \returns the image of the matrix, also called its column-space. The columns of the returned matrix
|
---|
| 169 | * will form a basis of the kernel.
|
---|
| 170 | *
|
---|
| 171 | * \param originalMatrix the original matrix, of which *this is the LU decomposition.
|
---|
| 172 | * The reason why it is needed to pass it here, is that this allows
|
---|
| 173 | * a large optimization, as otherwise this method would need to reconstruct it
|
---|
| 174 | * from the LU decomposition.
|
---|
| 175 | *
|
---|
| 176 | * \note If the image has dimension zero, then the returned matrix is a column-vector filled with zeros.
|
---|
| 177 | *
|
---|
| 178 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
| 179 | * For that, it uses the threshold value that you can control by calling
|
---|
| 180 | * setThreshold(const RealScalar&).
|
---|
| 181 | *
|
---|
| 182 | * Example: \include FullPivLU_image.cpp
|
---|
| 183 | * Output: \verbinclude FullPivLU_image.out
|
---|
| 184 | *
|
---|
| 185 | * \sa kernel()
|
---|
| 186 | */
|
---|
| 187 | inline const internal::image_retval<FullPivLU>
|
---|
| 188 | image(const MatrixType& originalMatrix) const
|
---|
| 189 | {
|
---|
| 190 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 191 | return internal::image_retval<FullPivLU>(*this, originalMatrix);
|
---|
| 192 | }
|
---|
| 193 |
|
---|
| 194 | /** \return a solution x to the equation Ax=b, where A is the matrix of which
|
---|
| 195 | * *this is the LU decomposition.
|
---|
| 196 | *
|
---|
| 197 | * \param b the right-hand-side of the equation to solve. Can be a vector or a matrix,
|
---|
| 198 | * the only requirement in order for the equation to make sense is that
|
---|
| 199 | * b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition.
|
---|
| 200 | *
|
---|
| 201 | * \returns a solution.
|
---|
| 202 | *
|
---|
| 203 | * \note_about_checking_solutions
|
---|
| 204 | *
|
---|
| 205 | * \note_about_arbitrary_choice_of_solution
|
---|
| 206 | * \note_about_using_kernel_to_study_multiple_solutions
|
---|
| 207 | *
|
---|
| 208 | * Example: \include FullPivLU_solve.cpp
|
---|
| 209 | * Output: \verbinclude FullPivLU_solve.out
|
---|
| 210 | *
|
---|
| 211 | * \sa TriangularView::solve(), kernel(), inverse()
|
---|
| 212 | */
|
---|
| 213 | template<typename Rhs>
|
---|
| 214 | inline const internal::solve_retval<FullPivLU, Rhs>
|
---|
| 215 | solve(const MatrixBase<Rhs>& b) const
|
---|
| 216 | {
|
---|
| 217 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 218 | return internal::solve_retval<FullPivLU, Rhs>(*this, b.derived());
|
---|
| 219 | }
|
---|
| 220 |
|
---|
| 221 | /** \returns the determinant of the matrix of which
|
---|
| 222 | * *this is the LU decomposition. It has only linear complexity
|
---|
| 223 | * (that is, O(n) where n is the dimension of the square matrix)
|
---|
| 224 | * as the LU decomposition has already been computed.
|
---|
| 225 | *
|
---|
| 226 | * \note This is only for square matrices.
|
---|
| 227 | *
|
---|
| 228 | * \note For fixed-size matrices of size up to 4, MatrixBase::determinant() offers
|
---|
| 229 | * optimized paths.
|
---|
| 230 | *
|
---|
| 231 | * \warning a determinant can be very big or small, so for matrices
|
---|
| 232 | * of large enough dimension, there is a risk of overflow/underflow.
|
---|
| 233 | *
|
---|
| 234 | * \sa MatrixBase::determinant()
|
---|
| 235 | */
|
---|
| 236 | typename internal::traits<MatrixType>::Scalar determinant() const;
|
---|
| 237 |
|
---|
| 238 | /** Allows to prescribe a threshold to be used by certain methods, such as rank(),
|
---|
| 239 | * who need to determine when pivots are to be considered nonzero. This is not used for the
|
---|
| 240 | * LU decomposition itself.
|
---|
| 241 | *
|
---|
| 242 | * When it needs to get the threshold value, Eigen calls threshold(). By default, this
|
---|
| 243 | * uses a formula to automatically determine a reasonable threshold.
|
---|
| 244 | * Once you have called the present method setThreshold(const RealScalar&),
|
---|
| 245 | * your value is used instead.
|
---|
| 246 | *
|
---|
| 247 | * \param threshold The new value to use as the threshold.
|
---|
| 248 | *
|
---|
| 249 | * A pivot will be considered nonzero if its absolute value is strictly greater than
|
---|
| 250 | * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
|
---|
| 251 | * where maxpivot is the biggest pivot.
|
---|
| 252 | *
|
---|
| 253 | * If you want to come back to the default behavior, call setThreshold(Default_t)
|
---|
| 254 | */
|
---|
| 255 | FullPivLU& setThreshold(const RealScalar& threshold)
|
---|
| 256 | {
|
---|
| 257 | m_usePrescribedThreshold = true;
|
---|
| 258 | m_prescribedThreshold = threshold;
|
---|
| 259 | return *this;
|
---|
| 260 | }
|
---|
| 261 |
|
---|
| 262 | /** Allows to come back to the default behavior, letting Eigen use its default formula for
|
---|
| 263 | * determining the threshold.
|
---|
| 264 | *
|
---|
| 265 | * You should pass the special object Eigen::Default as parameter here.
|
---|
| 266 | * \code lu.setThreshold(Eigen::Default); \endcode
|
---|
| 267 | *
|
---|
| 268 | * See the documentation of setThreshold(const RealScalar&).
|
---|
| 269 | */
|
---|
| 270 | FullPivLU& setThreshold(Default_t)
|
---|
| 271 | {
|
---|
| 272 | m_usePrescribedThreshold = false;
|
---|
| 273 | return *this;
|
---|
| 274 | }
|
---|
| 275 |
|
---|
| 276 | /** Returns the threshold that will be used by certain methods such as rank().
|
---|
| 277 | *
|
---|
| 278 | * See the documentation of setThreshold(const RealScalar&).
|
---|
| 279 | */
|
---|
| 280 | RealScalar threshold() const
|
---|
| 281 | {
|
---|
| 282 | eigen_assert(m_isInitialized || m_usePrescribedThreshold);
|
---|
| 283 | return m_usePrescribedThreshold ? m_prescribedThreshold
|
---|
| 284 | // this formula comes from experimenting (see "LU precision tuning" thread on the list)
|
---|
| 285 | // and turns out to be identical to Higham's formula used already in LDLt.
|
---|
| 286 | : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize();
|
---|
| 287 | }
|
---|
| 288 |
|
---|
| 289 | /** \returns the rank of the matrix of which *this is the LU decomposition.
|
---|
| 290 | *
|
---|
| 291 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
| 292 | * For that, it uses the threshold value that you can control by calling
|
---|
| 293 | * setThreshold(const RealScalar&).
|
---|
| 294 | */
|
---|
| 295 | inline Index rank() const
|
---|
| 296 | {
|
---|
| 297 | using std::abs;
|
---|
| 298 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 299 | RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
|
---|
| 300 | Index result = 0;
|
---|
| 301 | for(Index i = 0; i < m_nonzero_pivots; ++i)
|
---|
| 302 | result += (abs(m_lu.coeff(i,i)) > premultiplied_threshold);
|
---|
| 303 | return result;
|
---|
| 304 | }
|
---|
| 305 |
|
---|
| 306 | /** \returns the dimension of the kernel of the matrix of which *this is the LU decomposition.
|
---|
| 307 | *
|
---|
| 308 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
| 309 | * For that, it uses the threshold value that you can control by calling
|
---|
| 310 | * setThreshold(const RealScalar&).
|
---|
| 311 | */
|
---|
| 312 | inline Index dimensionOfKernel() const
|
---|
| 313 | {
|
---|
| 314 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 315 | return cols() - rank();
|
---|
| 316 | }
|
---|
| 317 |
|
---|
| 318 | /** \returns true if the matrix of which *this is the LU decomposition represents an injective
|
---|
| 319 | * linear map, i.e. has trivial kernel; false otherwise.
|
---|
| 320 | *
|
---|
| 321 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
| 322 | * For that, it uses the threshold value that you can control by calling
|
---|
| 323 | * setThreshold(const RealScalar&).
|
---|
| 324 | */
|
---|
| 325 | inline bool isInjective() const
|
---|
| 326 | {
|
---|
| 327 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 328 | return rank() == cols();
|
---|
| 329 | }
|
---|
| 330 |
|
---|
| 331 | /** \returns true if the matrix of which *this is the LU decomposition represents a surjective
|
---|
| 332 | * linear map; false otherwise.
|
---|
| 333 | *
|
---|
| 334 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
| 335 | * For that, it uses the threshold value that you can control by calling
|
---|
| 336 | * setThreshold(const RealScalar&).
|
---|
| 337 | */
|
---|
| 338 | inline bool isSurjective() const
|
---|
| 339 | {
|
---|
| 340 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 341 | return rank() == rows();
|
---|
| 342 | }
|
---|
| 343 |
|
---|
| 344 | /** \returns true if the matrix of which *this is the LU decomposition is invertible.
|
---|
| 345 | *
|
---|
| 346 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
| 347 | * For that, it uses the threshold value that you can control by calling
|
---|
| 348 | * setThreshold(const RealScalar&).
|
---|
| 349 | */
|
---|
| 350 | inline bool isInvertible() const
|
---|
| 351 | {
|
---|
| 352 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 353 | return isInjective() && (m_lu.rows() == m_lu.cols());
|
---|
| 354 | }
|
---|
| 355 |
|
---|
| 356 | /** \returns the inverse of the matrix of which *this is the LU decomposition.
|
---|
| 357 | *
|
---|
| 358 | * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
|
---|
| 359 | * Use isInvertible() to first determine whether this matrix is invertible.
|
---|
| 360 | *
|
---|
| 361 | * \sa MatrixBase::inverse()
|
---|
| 362 | */
|
---|
| 363 | inline const internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType> inverse() const
|
---|
| 364 | {
|
---|
| 365 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 366 | eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
|
---|
| 367 | return internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType>
|
---|
| 368 | (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
|
---|
| 369 | }
|
---|
| 370 |
|
---|
| 371 | MatrixType reconstructedMatrix() const;
|
---|
| 372 |
|
---|
| 373 | inline Index rows() const { return m_lu.rows(); }
|
---|
| 374 | inline Index cols() const { return m_lu.cols(); }
|
---|
| 375 |
|
---|
| 376 | protected:
|
---|
| 377 |
|
---|
| 378 | static void check_template_parameters()
|
---|
| 379 | {
|
---|
| 380 | EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
|
---|
| 381 | }
|
---|
| 382 |
|
---|
| 383 | MatrixType m_lu;
|
---|
| 384 | PermutationPType m_p;
|
---|
| 385 | PermutationQType m_q;
|
---|
| 386 | IntColVectorType m_rowsTranspositions;
|
---|
| 387 | IntRowVectorType m_colsTranspositions;
|
---|
| 388 | Index m_det_pq, m_nonzero_pivots;
|
---|
| 389 | RealScalar m_maxpivot, m_prescribedThreshold;
|
---|
| 390 | bool m_isInitialized, m_usePrescribedThreshold;
|
---|
| 391 | };
|
---|
| 392 |
|
---|
| 393 | template<typename MatrixType>
|
---|
| 394 | FullPivLU<MatrixType>::FullPivLU()
|
---|
| 395 | : m_isInitialized(false), m_usePrescribedThreshold(false)
|
---|
| 396 | {
|
---|
| 397 | }
|
---|
| 398 |
|
---|
| 399 | template<typename MatrixType>
|
---|
| 400 | FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
|
---|
| 401 | : m_lu(rows, cols),
|
---|
| 402 | m_p(rows),
|
---|
| 403 | m_q(cols),
|
---|
| 404 | m_rowsTranspositions(rows),
|
---|
| 405 | m_colsTranspositions(cols),
|
---|
| 406 | m_isInitialized(false),
|
---|
| 407 | m_usePrescribedThreshold(false)
|
---|
| 408 | {
|
---|
| 409 | }
|
---|
| 410 |
|
---|
| 411 | template<typename MatrixType>
|
---|
| 412 | FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
|
---|
| 413 | : m_lu(matrix.rows(), matrix.cols()),
|
---|
| 414 | m_p(matrix.rows()),
|
---|
| 415 | m_q(matrix.cols()),
|
---|
| 416 | m_rowsTranspositions(matrix.rows()),
|
---|
| 417 | m_colsTranspositions(matrix.cols()),
|
---|
| 418 | m_isInitialized(false),
|
---|
| 419 | m_usePrescribedThreshold(false)
|
---|
| 420 | {
|
---|
| 421 | compute(matrix);
|
---|
| 422 | }
|
---|
| 423 |
|
---|
| 424 | template<typename MatrixType>
|
---|
| 425 | FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
|
---|
| 426 | {
|
---|
| 427 | check_template_parameters();
|
---|
| 428 |
|
---|
| 429 | // the permutations are stored as int indices, so just to be sure:
|
---|
| 430 | eigen_assert(matrix.rows()<=NumTraits<int>::highest() && matrix.cols()<=NumTraits<int>::highest());
|
---|
| 431 |
|
---|
| 432 | m_isInitialized = true;
|
---|
| 433 | m_lu = matrix;
|
---|
| 434 |
|
---|
| 435 | const Index size = matrix.diagonalSize();
|
---|
| 436 | const Index rows = matrix.rows();
|
---|
| 437 | const Index cols = matrix.cols();
|
---|
| 438 |
|
---|
| 439 | // will store the transpositions, before we accumulate them at the end.
|
---|
| 440 | // can't accumulate on-the-fly because that will be done in reverse order for the rows.
|
---|
| 441 | m_rowsTranspositions.resize(matrix.rows());
|
---|
| 442 | m_colsTranspositions.resize(matrix.cols());
|
---|
| 443 | Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
|
---|
| 444 |
|
---|
| 445 | m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
|
---|
| 446 | m_maxpivot = RealScalar(0);
|
---|
| 447 |
|
---|
| 448 | for(Index k = 0; k < size; ++k)
|
---|
| 449 | {
|
---|
| 450 | // First, we need to find the pivot.
|
---|
| 451 |
|
---|
| 452 | // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
|
---|
| 453 | Index row_of_biggest_in_corner, col_of_biggest_in_corner;
|
---|
| 454 | RealScalar biggest_in_corner;
|
---|
| 455 | biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
|
---|
| 456 | .cwiseAbs()
|
---|
| 457 | .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
|
---|
| 458 | row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
|
---|
| 459 | col_of_biggest_in_corner += k; // need to add k to them.
|
---|
| 460 |
|
---|
| 461 | if(biggest_in_corner==RealScalar(0))
|
---|
| 462 | {
|
---|
| 463 | // before exiting, make sure to initialize the still uninitialized transpositions
|
---|
| 464 | // in a sane state without destroying what we already have.
|
---|
| 465 | m_nonzero_pivots = k;
|
---|
| 466 | for(Index i = k; i < size; ++i)
|
---|
| 467 | {
|
---|
| 468 | m_rowsTranspositions.coeffRef(i) = i;
|
---|
| 469 | m_colsTranspositions.coeffRef(i) = i;
|
---|
| 470 | }
|
---|
| 471 | break;
|
---|
| 472 | }
|
---|
| 473 |
|
---|
| 474 | if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
|
---|
| 475 |
|
---|
| 476 | // Now that we've found the pivot, we need to apply the row/col swaps to
|
---|
| 477 | // bring it to the location (k,k).
|
---|
| 478 |
|
---|
| 479 | m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
|
---|
| 480 | m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
|
---|
| 481 | if(k != row_of_biggest_in_corner) {
|
---|
| 482 | m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
|
---|
| 483 | ++number_of_transpositions;
|
---|
| 484 | }
|
---|
| 485 | if(k != col_of_biggest_in_corner) {
|
---|
| 486 | m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
|
---|
| 487 | ++number_of_transpositions;
|
---|
| 488 | }
|
---|
| 489 |
|
---|
| 490 | // Now that the pivot is at the right location, we update the remaining
|
---|
| 491 | // bottom-right corner by Gaussian elimination.
|
---|
| 492 |
|
---|
| 493 | if(k<rows-1)
|
---|
| 494 | m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
|
---|
| 495 | if(k<size-1)
|
---|
| 496 | m_lu.block(k+1,k+1,rows-k-1,cols-k-1).noalias() -= m_lu.col(k).tail(rows-k-1) * m_lu.row(k).tail(cols-k-1);
|
---|
| 497 | }
|
---|
| 498 |
|
---|
| 499 | // the main loop is over, we still have to accumulate the transpositions to find the
|
---|
| 500 | // permutations P and Q
|
---|
| 501 |
|
---|
| 502 | m_p.setIdentity(rows);
|
---|
| 503 | for(Index k = size-1; k >= 0; --k)
|
---|
| 504 | m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
|
---|
| 505 |
|
---|
| 506 | m_q.setIdentity(cols);
|
---|
| 507 | for(Index k = 0; k < size; ++k)
|
---|
| 508 | m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
|
---|
| 509 |
|
---|
| 510 | m_det_pq = (number_of_transpositions%2) ? -1 : 1;
|
---|
| 511 | return *this;
|
---|
| 512 | }
|
---|
| 513 |
|
---|
| 514 | template<typename MatrixType>
|
---|
| 515 | typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
|
---|
| 516 | {
|
---|
| 517 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 518 | eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
|
---|
| 519 | return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
|
---|
| 520 | }
|
---|
| 521 |
|
---|
| 522 | /** \returns the matrix represented by the decomposition,
|
---|
| 523 | * i.e., it returns the product: \f$ P^{-1} L U Q^{-1} \f$.
|
---|
| 524 | * This function is provided for debug purposes. */
|
---|
| 525 | template<typename MatrixType>
|
---|
| 526 | MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
|
---|
| 527 | {
|
---|
| 528 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
| 529 | const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
|
---|
| 530 | // LU
|
---|
| 531 | MatrixType res(m_lu.rows(),m_lu.cols());
|
---|
| 532 | // FIXME the .toDenseMatrix() should not be needed...
|
---|
| 533 | res = m_lu.leftCols(smalldim)
|
---|
| 534 | .template triangularView<UnitLower>().toDenseMatrix()
|
---|
| 535 | * m_lu.topRows(smalldim)
|
---|
| 536 | .template triangularView<Upper>().toDenseMatrix();
|
---|
| 537 |
|
---|
| 538 | // P^{-1}(LU)
|
---|
| 539 | res = m_p.inverse() * res;
|
---|
| 540 |
|
---|
| 541 | // (P^{-1}LU)Q^{-1}
|
---|
| 542 | res = res * m_q.inverse();
|
---|
| 543 |
|
---|
| 544 | return res;
|
---|
| 545 | }
|
---|
| 546 |
|
---|
| 547 | /********* Implementation of kernel() **************************************************/
|
---|
| 548 |
|
---|
| 549 | namespace internal {
|
---|
| 550 | template<typename _MatrixType>
|
---|
| 551 | struct kernel_retval<FullPivLU<_MatrixType> >
|
---|
| 552 | : kernel_retval_base<FullPivLU<_MatrixType> >
|
---|
| 553 | {
|
---|
| 554 | EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>)
|
---|
| 555 |
|
---|
| 556 | enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
|
---|
| 557 | MatrixType::MaxColsAtCompileTime,
|
---|
| 558 | MatrixType::MaxRowsAtCompileTime)
|
---|
| 559 | };
|
---|
| 560 |
|
---|
| 561 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
| 562 | {
|
---|
| 563 | using std::abs;
|
---|
| 564 | const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
|
---|
| 565 | if(dimker == 0)
|
---|
| 566 | {
|
---|
| 567 | // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
|
---|
| 568 | // avoid crashing/asserting as that depends on floating point calculations. Let's
|
---|
| 569 | // just return a single column vector filled with zeros.
|
---|
| 570 | dst.setZero();
|
---|
| 571 | return;
|
---|
| 572 | }
|
---|
| 573 |
|
---|
| 574 | /* Let us use the following lemma:
|
---|
| 575 | *
|
---|
| 576 | * Lemma: If the matrix A has the LU decomposition PAQ = LU,
|
---|
| 577 | * then Ker A = Q(Ker U).
|
---|
| 578 | *
|
---|
| 579 | * Proof: trivial: just keep in mind that P, Q, L are invertible.
|
---|
| 580 | */
|
---|
| 581 |
|
---|
| 582 | /* Thus, all we need to do is to compute Ker U, and then apply Q.
|
---|
| 583 | *
|
---|
| 584 | * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
|
---|
| 585 | * Thus, the diagonal of U ends with exactly
|
---|
| 586 | * dimKer zero's. Let us use that to construct dimKer linearly
|
---|
| 587 | * independent vectors in Ker U.
|
---|
| 588 | */
|
---|
| 589 |
|
---|
| 590 | Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
|
---|
| 591 | RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
|
---|
| 592 | Index p = 0;
|
---|
| 593 | for(Index i = 0; i < dec().nonzeroPivots(); ++i)
|
---|
| 594 | if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
|
---|
| 595 | pivots.coeffRef(p++) = i;
|
---|
| 596 | eigen_internal_assert(p == rank());
|
---|
| 597 |
|
---|
| 598 | // we construct a temporaty trapezoid matrix m, by taking the U matrix and
|
---|
| 599 | // permuting the rows and cols to bring the nonnegligible pivots to the top of
|
---|
| 600 | // the main diagonal. We need that to be able to apply our triangular solvers.
|
---|
| 601 | // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
|
---|
| 602 | Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
|
---|
| 603 | MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
|
---|
| 604 | m(dec().matrixLU().block(0, 0, rank(), cols));
|
---|
| 605 | for(Index i = 0; i < rank(); ++i)
|
---|
| 606 | {
|
---|
| 607 | if(i) m.row(i).head(i).setZero();
|
---|
| 608 | m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
|
---|
| 609 | }
|
---|
| 610 | m.block(0, 0, rank(), rank());
|
---|
| 611 | m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
|
---|
| 612 | for(Index i = 0; i < rank(); ++i)
|
---|
| 613 | m.col(i).swap(m.col(pivots.coeff(i)));
|
---|
| 614 |
|
---|
| 615 | // ok, we have our trapezoid matrix, we can apply the triangular solver.
|
---|
| 616 | // notice that the math behind this suggests that we should apply this to the
|
---|
| 617 | // negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
|
---|
| 618 | m.topLeftCorner(rank(), rank())
|
---|
| 619 | .template triangularView<Upper>().solveInPlace(
|
---|
| 620 | m.topRightCorner(rank(), dimker)
|
---|
| 621 | );
|
---|
| 622 |
|
---|
| 623 | // now we must undo the column permutation that we had applied!
|
---|
| 624 | for(Index i = rank()-1; i >= 0; --i)
|
---|
| 625 | m.col(i).swap(m.col(pivots.coeff(i)));
|
---|
| 626 |
|
---|
| 627 | // see the negative sign in the next line, that's what we were talking about above.
|
---|
| 628 | for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
|
---|
| 629 | for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
|
---|
| 630 | for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
|
---|
| 631 | }
|
---|
| 632 | };
|
---|
| 633 |
|
---|
| 634 | /***** Implementation of image() *****************************************************/
|
---|
| 635 |
|
---|
| 636 | template<typename _MatrixType>
|
---|
| 637 | struct image_retval<FullPivLU<_MatrixType> >
|
---|
| 638 | : image_retval_base<FullPivLU<_MatrixType> >
|
---|
| 639 | {
|
---|
| 640 | EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
|
---|
| 641 |
|
---|
| 642 | enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
|
---|
| 643 | MatrixType::MaxColsAtCompileTime,
|
---|
| 644 | MatrixType::MaxRowsAtCompileTime)
|
---|
| 645 | };
|
---|
| 646 |
|
---|
| 647 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
| 648 | {
|
---|
| 649 | using std::abs;
|
---|
| 650 | if(rank() == 0)
|
---|
| 651 | {
|
---|
| 652 | // The Image is just {0}, so it doesn't have a basis properly speaking, but let's
|
---|
| 653 | // avoid crashing/asserting as that depends on floating point calculations. Let's
|
---|
| 654 | // just return a single column vector filled with zeros.
|
---|
| 655 | dst.setZero();
|
---|
| 656 | return;
|
---|
| 657 | }
|
---|
| 658 |
|
---|
| 659 | Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
|
---|
| 660 | RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
|
---|
| 661 | Index p = 0;
|
---|
| 662 | for(Index i = 0; i < dec().nonzeroPivots(); ++i)
|
---|
| 663 | if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
|
---|
| 664 | pivots.coeffRef(p++) = i;
|
---|
| 665 | eigen_internal_assert(p == rank());
|
---|
| 666 |
|
---|
| 667 | for(Index i = 0; i < rank(); ++i)
|
---|
| 668 | dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
|
---|
| 669 | }
|
---|
| 670 | };
|
---|
| 671 |
|
---|
| 672 | /***** Implementation of solve() *****************************************************/
|
---|
| 673 |
|
---|
| 674 | template<typename _MatrixType, typename Rhs>
|
---|
| 675 | struct solve_retval<FullPivLU<_MatrixType>, Rhs>
|
---|
| 676 | : solve_retval_base<FullPivLU<_MatrixType>, Rhs>
|
---|
| 677 | {
|
---|
| 678 | EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs)
|
---|
| 679 |
|
---|
| 680 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
| 681 | {
|
---|
| 682 | /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
|
---|
| 683 | * So we proceed as follows:
|
---|
| 684 | * Step 1: compute c = P * rhs.
|
---|
| 685 | * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
|
---|
| 686 | * Step 3: replace c by the solution x to Ux = c. May or may not exist.
|
---|
| 687 | * Step 4: result = Q * c;
|
---|
| 688 | */
|
---|
| 689 |
|
---|
| 690 | const Index rows = dec().rows(), cols = dec().cols(),
|
---|
| 691 | nonzero_pivots = dec().rank();
|
---|
| 692 | eigen_assert(rhs().rows() == rows);
|
---|
| 693 | const Index smalldim = (std::min)(rows, cols);
|
---|
| 694 |
|
---|
| 695 | if(nonzero_pivots == 0)
|
---|
| 696 | {
|
---|
| 697 | dst.setZero();
|
---|
| 698 | return;
|
---|
| 699 | }
|
---|
| 700 |
|
---|
| 701 | typename Rhs::PlainObject c(rhs().rows(), rhs().cols());
|
---|
| 702 |
|
---|
| 703 | // Step 1
|
---|
| 704 | c = dec().permutationP() * rhs();
|
---|
| 705 |
|
---|
| 706 | // Step 2
|
---|
| 707 | dec().matrixLU()
|
---|
| 708 | .topLeftCorner(smalldim,smalldim)
|
---|
| 709 | .template triangularView<UnitLower>()
|
---|
| 710 | .solveInPlace(c.topRows(smalldim));
|
---|
| 711 | if(rows>cols)
|
---|
| 712 | {
|
---|
| 713 | c.bottomRows(rows-cols)
|
---|
| 714 | -= dec().matrixLU().bottomRows(rows-cols)
|
---|
| 715 | * c.topRows(cols);
|
---|
| 716 | }
|
---|
| 717 |
|
---|
| 718 | // Step 3
|
---|
| 719 | dec().matrixLU()
|
---|
| 720 | .topLeftCorner(nonzero_pivots, nonzero_pivots)
|
---|
| 721 | .template triangularView<Upper>()
|
---|
| 722 | .solveInPlace(c.topRows(nonzero_pivots));
|
---|
| 723 |
|
---|
| 724 | // Step 4
|
---|
| 725 | for(Index i = 0; i < nonzero_pivots; ++i)
|
---|
| 726 | dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i);
|
---|
| 727 | for(Index i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
|
---|
| 728 | dst.row(dec().permutationQ().indices().coeff(i)).setZero();
|
---|
| 729 | }
|
---|
| 730 | };
|
---|
| 731 |
|
---|
| 732 | } // end namespace internal
|
---|
| 733 |
|
---|
| 734 | /******* MatrixBase methods *****************************************************************/
|
---|
| 735 |
|
---|
| 736 | /** \lu_module
|
---|
| 737 | *
|
---|
| 738 | * \return the full-pivoting LU decomposition of \c *this.
|
---|
| 739 | *
|
---|
| 740 | * \sa class FullPivLU
|
---|
| 741 | */
|
---|
| 742 | template<typename Derived>
|
---|
| 743 | inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
|
---|
| 744 | MatrixBase<Derived>::fullPivLu() const
|
---|
| 745 | {
|
---|
| 746 | return FullPivLU<PlainObject>(eval());
|
---|
| 747 | }
|
---|
| 748 |
|
---|
| 749 | } // end namespace Eigen
|
---|
| 750 |
|
---|
| 751 | #endif // EIGEN_LU_H
|
---|