Michael Ernst Peter / Eigen

Eigen libary for mbed

src/LU/FullPivLU.h@1:3b8049da21b8, 2019-09-24 (annotated)

Committer:
jsoh91
Date:
Tue Sep 24 00:18:23 2019 +0000
Revision:
1:3b8049da21b8
Parent:
0:13a5d365ba16 
ignore and revise some of error parts

Who changed what in which revision?

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