Michael Ernst Peter / Eigen

src/LU/PartialPivLU.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 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
ykuroda 0:13a5d365ba16 6 //
ykuroda 0:13a5d365ba16 7 // This Source Code Form is subject to the terms of the Mozilla
ykuroda 0:13a5d365ba16 8 // Public License v. 2.0. If a copy of the MPL was not distributed
ykuroda 0:13a5d365ba16 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
ykuroda 0:13a5d365ba16 10
ykuroda 0:13a5d365ba16 11 #ifndef EIGEN_PARTIALLU_H
ykuroda 0:13a5d365ba16 12 #define EIGEN_PARTIALLU_H
ykuroda 0:13a5d365ba16 13
ykuroda 0:13a5d365ba16 14 namespace Eigen {
ykuroda 0:13a5d365ba16 15
ykuroda 0:13a5d365ba16 16 /** \ingroup LU_Module
ykuroda 0:13a5d365ba16 17 *
ykuroda 0:13a5d365ba16 18 * \class PartialPivLU
ykuroda 0:13a5d365ba16 19 *
ykuroda 0:13a5d365ba16 20 * \brief LU decomposition of a matrix with partial pivoting, and related features
ykuroda 0:13a5d365ba16 21 *
ykuroda 0:13a5d365ba16 22 * \param MatrixType the type of the matrix of which we are computing the LU decomposition
ykuroda 0:13a5d365ba16 23 *
ykuroda 0:13a5d365ba16 24 * This class represents a LU decomposition of a \b square \b invertible matrix, with partial pivoting: the matrix A
ykuroda 0:13a5d365ba16 25 * is decomposed as A = PLU where L is unit-lower-triangular, U is upper-triangular, and P
ykuroda 0:13a5d365ba16 26 * is a permutation matrix.
ykuroda 0:13a5d365ba16 27 *
ykuroda 0:13a5d365ba16 28 * Typically, partial pivoting LU decomposition is only considered numerically stable for square invertible
ykuroda 0:13a5d365ba16 29 * matrices. Thus LAPACK's dgesv and dgesvx require the matrix to be square and invertible. The present class
ykuroda 0:13a5d365ba16 30 * does the same. It will assert that the matrix is square, but it won't (actually it can't) check that the
ykuroda 0:13a5d365ba16 31 * matrix is invertible: it is your task to check that you only use this decomposition on invertible matrices.
ykuroda 0:13a5d365ba16 32 *
ykuroda 0:13a5d365ba16 33 * The guaranteed safe alternative, working for all matrices, is the full pivoting LU decomposition, provided
ykuroda 0:13a5d365ba16 34 * by class FullPivLU.
ykuroda 0:13a5d365ba16 35 *
ykuroda 0:13a5d365ba16 36 * This is \b not a rank-revealing LU decomposition. Many features are intentionally absent from this class,
ykuroda 0:13a5d365ba16 37 * such as rank computation. If you need these features, use class FullPivLU.
ykuroda 0:13a5d365ba16 38 *
ykuroda 0:13a5d365ba16 39 * This LU decomposition is suitable to invert invertible matrices. It is what MatrixBase::inverse() uses
ykuroda 0:13a5d365ba16 40 * in the general case.
ykuroda 0:13a5d365ba16 41 * On the other hand, it is \b not suitable to determine whether a given matrix is invertible.
ykuroda 0:13a5d365ba16 42 *
ykuroda 0:13a5d365ba16 43 * The data of the LU decomposition can be directly accessed through the methods matrixLU(), permutationP().
ykuroda 0:13a5d365ba16 44 *
ykuroda 0:13a5d365ba16 45 * \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU
ykuroda 0:13a5d365ba16 46 */
ykuroda 0:13a5d365ba16 47 template<typename _MatrixType> class PartialPivLU
ykuroda 0:13a5d365ba16 48 {
ykuroda 0:13a5d365ba16 49 public:
ykuroda 0:13a5d365ba16 50
ykuroda 0:13a5d365ba16 51 typedef _MatrixType MatrixType;
ykuroda 0:13a5d365ba16 52 enum {
ykuroda 0:13a5d365ba16 53 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 54 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 55 Options = MatrixType::Options,
ykuroda 0:13a5d365ba16 56 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
ykuroda 0:13a5d365ba16 57 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
ykuroda 0:13a5d365ba16 58 };
ykuroda 0:13a5d365ba16 59 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 60 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
ykuroda 0:13a5d365ba16 61 typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
ykuroda 0:13a5d365ba16 62 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 63 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
ykuroda 0:13a5d365ba16 64 typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
ykuroda 0:13a5d365ba16 65
ykuroda 0:13a5d365ba16 66
ykuroda 0:13a5d365ba16 67 /**
ykuroda 0:13a5d365ba16 68 * \brief Default Constructor.
ykuroda 0:13a5d365ba16 69 *
ykuroda 0:13a5d365ba16 70 * The default constructor is useful in cases in which the user intends to
ykuroda 0:13a5d365ba16 71 * perform decompositions via PartialPivLU::compute(const MatrixType&).
ykuroda 0:13a5d365ba16 72 */
ykuroda 0:13a5d365ba16 73 PartialPivLU();
ykuroda 0:13a5d365ba16 74
ykuroda 0:13a5d365ba16 75 /** \brief Default Constructor with memory preallocation
ykuroda 0:13a5d365ba16 76 *
ykuroda 0:13a5d365ba16 77 * Like the default constructor but with preallocation of the internal data
ykuroda 0:13a5d365ba16 78 * according to the specified problem \a size.
ykuroda 0:13a5d365ba16 79 * \sa PartialPivLU()
ykuroda 0:13a5d365ba16 80 */
ykuroda 0:13a5d365ba16 81 PartialPivLU(Index size);
ykuroda 0:13a5d365ba16 82
ykuroda 0:13a5d365ba16 83 /** Constructor.
ykuroda 0:13a5d365ba16 84 *
ykuroda 0:13a5d365ba16 85 * \param matrix the matrix of which to compute the LU decomposition.
ykuroda 0:13a5d365ba16 86 *
ykuroda 0:13a5d365ba16 87 * \warning The matrix should have full rank (e.g. if it's square, it should be invertible).
ykuroda 0:13a5d365ba16 88 * If you need to deal with non-full rank, use class FullPivLU instead.
ykuroda 0:13a5d365ba16 89 */
ykuroda 0:13a5d365ba16 90 PartialPivLU(const MatrixType& matrix);
ykuroda 0:13a5d365ba16 91
ykuroda 0:13a5d365ba16 92 PartialPivLU& compute(const MatrixType& matrix);
ykuroda 0:13a5d365ba16 93
ykuroda 0:13a5d365ba16 94 /** \returns the LU decomposition matrix: the upper-triangular part is U, the
ykuroda 0:13a5d365ba16 95 * unit-lower-triangular part is L (at least for square matrices; in the non-square
ykuroda 0:13a5d365ba16 96 * case, special care is needed, see the documentation of class FullPivLU).
ykuroda 0:13a5d365ba16 97 *
ykuroda 0:13a5d365ba16 98 * \sa matrixL(), matrixU()
ykuroda 0:13a5d365ba16 99 */
ykuroda 0:13a5d365ba16 100 inline const MatrixType& matrixLU() const
ykuroda 0:13a5d365ba16 101 {
ykuroda 0:13a5d365ba16 102 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
ykuroda 0:13a5d365ba16 103 return m_lu;
ykuroda 0:13a5d365ba16 104 }
ykuroda 0:13a5d365ba16 105
ykuroda 0:13a5d365ba16 106 /** \returns the permutation matrix P.
ykuroda 0:13a5d365ba16 107 */
ykuroda 0:13a5d365ba16 108 inline const PermutationType& permutationP() const
ykuroda 0:13a5d365ba16 109 {
ykuroda 0:13a5d365ba16 110 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
ykuroda 0:13a5d365ba16 111 return m_p;
ykuroda 0:13a5d365ba16 112 }
ykuroda 0:13a5d365ba16 113
ykuroda 0:13a5d365ba16 114 /** This method returns the solution x to the equation Ax=b, where A is the matrix of which
ykuroda 0:13a5d365ba16 115 * *this is the LU decomposition.
ykuroda 0:13a5d365ba16 116 *
ykuroda 0:13a5d365ba16 117 * \param b the right-hand-side of the equation to solve. Can be a vector or a matrix,
ykuroda 0:13a5d365ba16 118 * the only requirement in order for the equation to make sense is that
ykuroda 0:13a5d365ba16 119 * b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition.
ykuroda 0:13a5d365ba16 120 *
ykuroda 0:13a5d365ba16 121 * \returns the solution.
ykuroda 0:13a5d365ba16 122 *
ykuroda 0:13a5d365ba16 123 * Example: \include PartialPivLU_solve.cpp
ykuroda 0:13a5d365ba16 124 * Output: \verbinclude PartialPivLU_solve.out
ykuroda 0:13a5d365ba16 125 *
ykuroda 0:13a5d365ba16 126 * Since this PartialPivLU class assumes anyway that the matrix A is invertible, the solution
ykuroda 0:13a5d365ba16 127 * theoretically exists and is unique regardless of b.
ykuroda 0:13a5d365ba16 128 *
ykuroda 0:13a5d365ba16 129 * \sa TriangularView::solve(), inverse(), computeInverse()
ykuroda 0:13a5d365ba16 130 */
ykuroda 0:13a5d365ba16 131 template<typename Rhs>
ykuroda 0:13a5d365ba16 132 inline const internal::solve_retval<PartialPivLU, Rhs>
ykuroda 0:13a5d365ba16 133 solve(const MatrixBase<Rhs>& b) const
ykuroda 0:13a5d365ba16 134 {
ykuroda 0:13a5d365ba16 135 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
ykuroda 0:13a5d365ba16 136 return internal::solve_retval<PartialPivLU, Rhs>(*this, b.derived());
ykuroda 0:13a5d365ba16 137 }
ykuroda 0:13a5d365ba16 138
ykuroda 0:13a5d365ba16 139 /** \returns the inverse of the matrix of which *this is the LU decomposition.
ykuroda 0:13a5d365ba16 140 *
ykuroda 0:13a5d365ba16 141 * \warning The matrix being decomposed here is assumed to be invertible. If you need to check for
ykuroda 0:13a5d365ba16 142 * invertibility, use class FullPivLU instead.
ykuroda 0:13a5d365ba16 143 *
ykuroda 0:13a5d365ba16 144 * \sa MatrixBase::inverse(), LU::inverse()
ykuroda 0:13a5d365ba16 145 */
ykuroda 0:13a5d365ba16 146 inline const internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> inverse() const
ykuroda 0:13a5d365ba16 147 {
ykuroda 0:13a5d365ba16 148 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
ykuroda 0:13a5d365ba16 149 return internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType>
ykuroda 0:13a5d365ba16 150 (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
ykuroda 0:13a5d365ba16 151 }
ykuroda 0:13a5d365ba16 152
ykuroda 0:13a5d365ba16 153 /** \returns the determinant of the matrix of which
ykuroda 0:13a5d365ba16 154 * *this is the LU decomposition. It has only linear complexity
ykuroda 0:13a5d365ba16 155 * (that is, O(n) where n is the dimension of the square matrix)
ykuroda 0:13a5d365ba16 156 * as the LU decomposition has already been computed.
ykuroda 0:13a5d365ba16 157 *
ykuroda 0:13a5d365ba16 158 * \note For fixed-size matrices of size up to 4, MatrixBase::determinant() offers
ykuroda 0:13a5d365ba16 159 * optimized paths.
ykuroda 0:13a5d365ba16 160 *
ykuroda 0:13a5d365ba16 161 * \warning a determinant can be very big or small, so for matrices
ykuroda 0:13a5d365ba16 162 * of large enough dimension, there is a risk of overflow/underflow.
ykuroda 0:13a5d365ba16 163 *
ykuroda 0:13a5d365ba16 164 * \sa MatrixBase::determinant()
ykuroda 0:13a5d365ba16 165 */
ykuroda 0:13a5d365ba16 166 typename internal::traits<MatrixType>::Scalar determinant() const;
ykuroda 0:13a5d365ba16 167
ykuroda 0:13a5d365ba16 168 MatrixType reconstructedMatrix() const;
ykuroda 0:13a5d365ba16 169
ykuroda 0:13a5d365ba16 170 inline Index rows() const { return m_lu.rows(); }
ykuroda 0:13a5d365ba16 171 inline Index cols() const { return m_lu.cols(); }
ykuroda 0:13a5d365ba16 172
ykuroda 0:13a5d365ba16 173 protected:
ykuroda 0:13a5d365ba16 174
ykuroda 0:13a5d365ba16 175 static void check_template_parameters()
ykuroda 0:13a5d365ba16 176 {
ykuroda 0:13a5d365ba16 177 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
ykuroda 0:13a5d365ba16 178 }
ykuroda 0:13a5d365ba16 179
ykuroda 0:13a5d365ba16 180 MatrixType m_lu;
ykuroda 0:13a5d365ba16 181 PermutationType m_p;
ykuroda 0:13a5d365ba16 182 TranspositionType m_rowsTranspositions;
ykuroda 0:13a5d365ba16 183 Index m_det_p;
ykuroda 0:13a5d365ba16 184 bool m_isInitialized;
ykuroda 0:13a5d365ba16 185 };
ykuroda 0:13a5d365ba16 186
ykuroda 0:13a5d365ba16 187 template<typename MatrixType>
ykuroda 0:13a5d365ba16 188 PartialPivLU<MatrixType>::PartialPivLU()
ykuroda 0:13a5d365ba16 189 : m_lu(),
ykuroda 0:13a5d365ba16 190 m_p(),
ykuroda 0:13a5d365ba16 191 m_rowsTranspositions(),
ykuroda 0:13a5d365ba16 192 m_det_p(0),
ykuroda 0:13a5d365ba16 193 m_isInitialized(false)
ykuroda 0:13a5d365ba16 194 {
ykuroda 0:13a5d365ba16 195 }
ykuroda 0:13a5d365ba16 196
ykuroda 0:13a5d365ba16 197 template<typename MatrixType>
ykuroda 0:13a5d365ba16 198 PartialPivLU<MatrixType>::PartialPivLU(Index size)
ykuroda 0:13a5d365ba16 199 : m_lu(size, size),
ykuroda 0:13a5d365ba16 200 m_p(size),
ykuroda 0:13a5d365ba16 201 m_rowsTranspositions(size),
ykuroda 0:13a5d365ba16 202 m_det_p(0),
ykuroda 0:13a5d365ba16 203 m_isInitialized(false)
ykuroda 0:13a5d365ba16 204 {
ykuroda 0:13a5d365ba16 205 }
ykuroda 0:13a5d365ba16 206
ykuroda 0:13a5d365ba16 207 template<typename MatrixType>
ykuroda 0:13a5d365ba16 208 PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
ykuroda 0:13a5d365ba16 209 : m_lu(matrix.rows(), matrix.rows()),
ykuroda 0:13a5d365ba16 210 m_p(matrix.rows()),
ykuroda 0:13a5d365ba16 211 m_rowsTranspositions(matrix.rows()),
ykuroda 0:13a5d365ba16 212 m_det_p(0),
ykuroda 0:13a5d365ba16 213 m_isInitialized(false)
ykuroda 0:13a5d365ba16 214 {
ykuroda 0:13a5d365ba16 215 compute(matrix);
ykuroda 0:13a5d365ba16 216 }
ykuroda 0:13a5d365ba16 217
ykuroda 0:13a5d365ba16 218 namespace internal {
ykuroda 0:13a5d365ba16 219
ykuroda 0:13a5d365ba16 220 /** \internal This is the blocked version of fullpivlu_unblocked() */
ykuroda 0:13a5d365ba16 221 template<typename Scalar, int StorageOrder, typename PivIndex>
ykuroda 0:13a5d365ba16 222 struct partial_lu_impl
ykuroda 0:13a5d365ba16 223 {
ykuroda 0:13a5d365ba16 224 // FIXME add a stride to Map, so that the following mapping becomes easier,
ykuroda 0:13a5d365ba16 225 // another option would be to create an expression being able to automatically
ykuroda 0:13a5d365ba16 226 // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly
ykuroda 0:13a5d365ba16 227 // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix,
ykuroda 0:13a5d365ba16 228 // and Block.
ykuroda 0:13a5d365ba16 229 typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > MapLU;
ykuroda 0:13a5d365ba16 230 typedef Block<MapLU, Dynamic, Dynamic> MatrixType;
ykuroda 0:13a5d365ba16 231 typedef Block<MatrixType,Dynamic,Dynamic> BlockType;
ykuroda 0:13a5d365ba16 232 typedef typename MatrixType::RealScalar RealScalar;
ykuroda 0:13a5d365ba16 233 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 234
ykuroda 0:13a5d365ba16 235 /** \internal performs the LU decomposition in-place of the matrix \a lu
ykuroda 0:13a5d365ba16 236 * using an unblocked algorithm.
ykuroda 0:13a5d365ba16 237 *
ykuroda 0:13a5d365ba16 238 * In addition, this function returns the row transpositions in the
ykuroda 0:13a5d365ba16 239 * vector \a row_transpositions which must have a size equal to the number
ykuroda 0:13a5d365ba16 240 * of columns of the matrix \a lu, and an integer \a nb_transpositions
ykuroda 0:13a5d365ba16 241 * which returns the actual number of transpositions.
ykuroda 0:13a5d365ba16 242 *
ykuroda 0:13a5d365ba16 243 * \returns The index of the first pivot which is exactly zero if any, or a negative number otherwise.
ykuroda 0:13a5d365ba16 244 */
ykuroda 0:13a5d365ba16 245 static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
ykuroda 0:13a5d365ba16 246 {
ykuroda 0:13a5d365ba16 247 const Index rows = lu.rows();
ykuroda 0:13a5d365ba16 248 const Index cols = lu.cols();
ykuroda 0:13a5d365ba16 249 const Index size = (std::min)(rows,cols);
ykuroda 0:13a5d365ba16 250 nb_transpositions = 0;
ykuroda 0:13a5d365ba16 251 Index first_zero_pivot = -1;
ykuroda 0:13a5d365ba16 252 for(Index k = 0; k < size; ++k)
ykuroda 0:13a5d365ba16 253 {
ykuroda 0:13a5d365ba16 254 Index rrows = rows-k-1;
ykuroda 0:13a5d365ba16 255 Index rcols = cols-k-1;
ykuroda 0:13a5d365ba16 256
ykuroda 0:13a5d365ba16 257 Index row_of_biggest_in_col;
ykuroda 0:13a5d365ba16 258 RealScalar biggest_in_corner
ykuroda 0:13a5d365ba16 259 = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col);
ykuroda 0:13a5d365ba16 260 row_of_biggest_in_col += k;
ykuroda 0:13a5d365ba16 261
ykuroda 0:13a5d365ba16 262 row_transpositions[k] = PivIndex(row_of_biggest_in_col);
ykuroda 0:13a5d365ba16 263
ykuroda 0:13a5d365ba16 264 if(biggest_in_corner != RealScalar(0))
ykuroda 0:13a5d365ba16 265 {
ykuroda 0:13a5d365ba16 266 if(k != row_of_biggest_in_col)
ykuroda 0:13a5d365ba16 267 {
ykuroda 0:13a5d365ba16 268 lu.row(k).swap(lu.row(row_of_biggest_in_col));
ykuroda 0:13a5d365ba16 269 ++nb_transpositions;
ykuroda 0:13a5d365ba16 270 }
ykuroda 0:13a5d365ba16 271
ykuroda 0:13a5d365ba16 272 // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k)
ykuroda 0:13a5d365ba16 273 // overflow but not the actual quotient?
ykuroda 0:13a5d365ba16 274 lu.col(k).tail(rrows) /= lu.coeff(k,k);
ykuroda 0:13a5d365ba16 275 }
ykuroda 0:13a5d365ba16 276 else if(first_zero_pivot==-1)
ykuroda 0:13a5d365ba16 277 {
ykuroda 0:13a5d365ba16 278 // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
ykuroda 0:13a5d365ba16 279 // and continue the factorization such we still have A = PLU
ykuroda 0:13a5d365ba16 280 first_zero_pivot = k;
ykuroda 0:13a5d365ba16 281 }
ykuroda 0:13a5d365ba16 282
ykuroda 0:13a5d365ba16 283 if(k<rows-1)
ykuroda 0:13a5d365ba16 284 lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
ykuroda 0:13a5d365ba16 285 }
ykuroda 0:13a5d365ba16 286 return first_zero_pivot;
ykuroda 0:13a5d365ba16 287 }
ykuroda 0:13a5d365ba16 288
ykuroda 0:13a5d365ba16 289 /** \internal performs the LU decomposition in-place of the matrix represented
ykuroda 0:13a5d365ba16 290 * by the variables \a rows, \a cols, \a lu_data, and \a lu_stride using a
ykuroda 0:13a5d365ba16 291 * recursive, blocked algorithm.
ykuroda 0:13a5d365ba16 292 *
ykuroda 0:13a5d365ba16 293 * In addition, this function returns the row transpositions in the
ykuroda 0:13a5d365ba16 294 * vector \a row_transpositions which must have a size equal to the number
ykuroda 0:13a5d365ba16 295 * of columns of the matrix \a lu, and an integer \a nb_transpositions
ykuroda 0:13a5d365ba16 296 * which returns the actual number of transpositions.
ykuroda 0:13a5d365ba16 297 *
ykuroda 0:13a5d365ba16 298 * \returns The index of the first pivot which is exactly zero if any, or a negative number otherwise.
ykuroda 0:13a5d365ba16 299 *
ykuroda 0:13a5d365ba16 300 * \note This very low level interface using pointers, etc. is to:
ykuroda 0:13a5d365ba16 301 * 1 - reduce the number of instanciations to the strict minimum
ykuroda 0:13a5d365ba16 302 * 2 - avoid infinite recursion of the instanciations with Block<Block<Block<...> > >
ykuroda 0:13a5d365ba16 303 */
ykuroda 0:13a5d365ba16 304 static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
ykuroda 0:13a5d365ba16 305 {
ykuroda 0:13a5d365ba16 306 MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
ykuroda 0:13a5d365ba16 307 MatrixType lu(lu1,0,0,rows,cols);
ykuroda 0:13a5d365ba16 308
ykuroda 0:13a5d365ba16 309 const Index size = (std::min)(rows,cols);
ykuroda 0:13a5d365ba16 310
ykuroda 0:13a5d365ba16 311 // if the matrix is too small, no blocking:
ykuroda 0:13a5d365ba16 312 if(size<=16)
ykuroda 0:13a5d365ba16 313 {
ykuroda 0:13a5d365ba16 314 return unblocked_lu(lu, row_transpositions, nb_transpositions);
ykuroda 0:13a5d365ba16 315 }
ykuroda 0:13a5d365ba16 316
ykuroda 0:13a5d365ba16 317 // automatically adjust the number of subdivisions to the size
ykuroda 0:13a5d365ba16 318 // of the matrix so that there is enough sub blocks:
ykuroda 0:13a5d365ba16 319 Index blockSize;
ykuroda 0:13a5d365ba16 320 {
ykuroda 0:13a5d365ba16 321 blockSize = size/8;
ykuroda 0:13a5d365ba16 322 blockSize = (blockSize/16)*16;
ykuroda 0:13a5d365ba16 323 blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
ykuroda 0:13a5d365ba16 324 }
ykuroda 0:13a5d365ba16 325
ykuroda 0:13a5d365ba16 326 nb_transpositions = 0;
ykuroda 0:13a5d365ba16 327 Index first_zero_pivot = -1;
ykuroda 0:13a5d365ba16 328 for(Index k = 0; k < size; k+=blockSize)
ykuroda 0:13a5d365ba16 329 {
ykuroda 0:13a5d365ba16 330 Index bs = (std::min)(size-k,blockSize); // actual size of the block
ykuroda 0:13a5d365ba16 331 Index trows = rows - k - bs; // trailing rows
ykuroda 0:13a5d365ba16 332 Index tsize = size - k - bs; // trailing size
ykuroda 0:13a5d365ba16 333
ykuroda 0:13a5d365ba16 334 // partition the matrix:
ykuroda 0:13a5d365ba16 335 // A00 | A01 | A02
ykuroda 0:13a5d365ba16 336 // lu = A_0 | A_1 | A_2 = A10 | A11 | A12
ykuroda 0:13a5d365ba16 337 // A20 | A21 | A22
ykuroda 0:13a5d365ba16 338 BlockType A_0(lu,0,0,rows,k);
ykuroda 0:13a5d365ba16 339 BlockType A_2(lu,0,k+bs,rows,tsize);
ykuroda 0:13a5d365ba16 340 BlockType A11(lu,k,k,bs,bs);
ykuroda 0:13a5d365ba16 341 BlockType A12(lu,k,k+bs,bs,tsize);
ykuroda 0:13a5d365ba16 342 BlockType A21(lu,k+bs,k,trows,bs);
ykuroda 0:13a5d365ba16 343 BlockType A22(lu,k+bs,k+bs,trows,tsize);
ykuroda 0:13a5d365ba16 344
ykuroda 0:13a5d365ba16 345 PivIndex nb_transpositions_in_panel;
ykuroda 0:13a5d365ba16 346 // recursively call the blocked LU algorithm on [A11^T A21^T]^T
ykuroda 0:13a5d365ba16 347 // with a very small blocking size:
ykuroda 0:13a5d365ba16 348 Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
ykuroda 0:13a5d365ba16 349 row_transpositions+k, nb_transpositions_in_panel, 16);
ykuroda 0:13a5d365ba16 350 if(ret>=0 && first_zero_pivot==-1)
ykuroda 0:13a5d365ba16 351 first_zero_pivot = k+ret;
ykuroda 0:13a5d365ba16 352
ykuroda 0:13a5d365ba16 353 nb_transpositions += nb_transpositions_in_panel;
ykuroda 0:13a5d365ba16 354 // update permutations and apply them to A_0
ykuroda 0:13a5d365ba16 355 for(Index i=k; i<k+bs; ++i)
ykuroda 0:13a5d365ba16 356 {
ykuroda 0:13a5d365ba16 357 Index piv = (row_transpositions[i] += k);
ykuroda 0:13a5d365ba16 358 A_0.row(i).swap(A_0.row(piv));
ykuroda 0:13a5d365ba16 359 }
ykuroda 0:13a5d365ba16 360
ykuroda 0:13a5d365ba16 361 if(trows)
ykuroda 0:13a5d365ba16 362 {
ykuroda 0:13a5d365ba16 363 // apply permutations to A_2
ykuroda 0:13a5d365ba16 364 for(Index i=k;i<k+bs; ++i)
ykuroda 0:13a5d365ba16 365 A_2.row(i).swap(A_2.row(row_transpositions[i]));
ykuroda 0:13a5d365ba16 366
ykuroda 0:13a5d365ba16 367 // A12 = A11^-1 A12
ykuroda 0:13a5d365ba16 368 A11.template triangularView<UnitLower>().solveInPlace(A12);
ykuroda 0:13a5d365ba16 369
ykuroda 0:13a5d365ba16 370 A22.noalias() -= A21 * A12;
ykuroda 0:13a5d365ba16 371 }
ykuroda 0:13a5d365ba16 372 }
ykuroda 0:13a5d365ba16 373 return first_zero_pivot;
ykuroda 0:13a5d365ba16 374 }
ykuroda 0:13a5d365ba16 375 };
ykuroda 0:13a5d365ba16 376
ykuroda 0:13a5d365ba16 377 /** \internal performs the LU decomposition with partial pivoting in-place.
ykuroda 0:13a5d365ba16 378 */
ykuroda 0:13a5d365ba16 379 template<typename MatrixType, typename TranspositionType>
ykuroda 0:13a5d365ba16 380 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::Index& nb_transpositions)
ykuroda 0:13a5d365ba16 381 {
ykuroda 0:13a5d365ba16 382 eigen_assert(lu.cols() == row_transpositions.size());
ykuroda 0:13a5d365ba16 383 eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
ykuroda 0:13a5d365ba16 384
ykuroda 0:13a5d365ba16 385 partial_lu_impl
ykuroda 0:13a5d365ba16 386 <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::Index>
ykuroda 0:13a5d365ba16 387 ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
ykuroda 0:13a5d365ba16 388 }
ykuroda 0:13a5d365ba16 389
ykuroda 0:13a5d365ba16 390 } // end namespace internal
ykuroda 0:13a5d365ba16 391
ykuroda 0:13a5d365ba16 392 template<typename MatrixType>
ykuroda 0:13a5d365ba16 393 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix)
ykuroda 0:13a5d365ba16 394 {
ykuroda 0:13a5d365ba16 395 check_template_parameters();
ykuroda 0:13a5d365ba16 396
ykuroda 0:13a5d365ba16 397 // the row permutation is stored as int indices, so just to be sure:
ykuroda 0:13a5d365ba16 398 eigen_assert(matrix.rows()<NumTraits<int>::highest());
ykuroda 0:13a5d365ba16 399
ykuroda 0:13a5d365ba16 400 m_lu = matrix;
ykuroda 0:13a5d365ba16 401
ykuroda 0:13a5d365ba16 402 eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
ykuroda 0:13a5d365ba16 403 const Index size = matrix.rows();
ykuroda 0:13a5d365ba16 404
ykuroda 0:13a5d365ba16 405 m_rowsTranspositions.resize(size);
ykuroda 0:13a5d365ba16 406
ykuroda 0:13a5d365ba16 407 typename TranspositionType::Index nb_transpositions;
ykuroda 0:13a5d365ba16 408 internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
ykuroda 0:13a5d365ba16 409 m_det_p = (nb_transpositions%2) ? -1 : 1;
ykuroda 0:13a5d365ba16 410
ykuroda 0:13a5d365ba16 411 m_p = m_rowsTranspositions;
ykuroda 0:13a5d365ba16 412
ykuroda 0:13a5d365ba16 413 m_isInitialized = true;
ykuroda 0:13a5d365ba16 414 return *this;
ykuroda 0:13a5d365ba16 415 }
ykuroda 0:13a5d365ba16 416
ykuroda 0:13a5d365ba16 417 template<typename MatrixType>
ykuroda 0:13a5d365ba16 418 typename internal::traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
ykuroda 0:13a5d365ba16 419 {
ykuroda 0:13a5d365ba16 420 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
ykuroda 0:13a5d365ba16 421 return Scalar(m_det_p) * m_lu.diagonal().prod();
ykuroda 0:13a5d365ba16 422 }
ykuroda 0:13a5d365ba16 423
ykuroda 0:13a5d365ba16 424 /** \returns the matrix represented by the decomposition,
ykuroda 0:13a5d365ba16 425 * i.e., it returns the product: P^{-1} L U.
ykuroda 0:13a5d365ba16 426 * This function is provided for debug purpose. */
ykuroda 0:13a5d365ba16 427 template<typename MatrixType>
ykuroda 0:13a5d365ba16 428 MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
ykuroda 0:13a5d365ba16 429 {
ykuroda 0:13a5d365ba16 430 eigen_assert(m_isInitialized && "LU is not initialized.");
ykuroda 0:13a5d365ba16 431 // LU
ykuroda 0:13a5d365ba16 432 MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
ykuroda 0:13a5d365ba16 433 * m_lu.template triangularView<Upper>();
ykuroda 0:13a5d365ba16 434
ykuroda 0:13a5d365ba16 435 // P^{-1}(LU)
ykuroda 0:13a5d365ba16 436 res = m_p.inverse() * res;
ykuroda 0:13a5d365ba16 437
ykuroda 0:13a5d365ba16 438 return res;
ykuroda 0:13a5d365ba16 439 }
ykuroda 0:13a5d365ba16 440
ykuroda 0:13a5d365ba16 441 /***** Implementation of solve() *****************************************************/
ykuroda 0:13a5d365ba16 442
ykuroda 0:13a5d365ba16 443 namespace internal {
ykuroda 0:13a5d365ba16 444
ykuroda 0:13a5d365ba16 445 template<typename _MatrixType, typename Rhs>
ykuroda 0:13a5d365ba16 446 struct solve_retval<PartialPivLU<_MatrixType>, Rhs>
ykuroda 0:13a5d365ba16 447 : solve_retval_base<PartialPivLU<_MatrixType>, Rhs>
ykuroda 0:13a5d365ba16 448 {
ykuroda 0:13a5d365ba16 449 EIGEN_MAKE_SOLVE_HELPERS(PartialPivLU<_MatrixType>,Rhs)
ykuroda 0:13a5d365ba16 450
ykuroda 0:13a5d365ba16 451 template<typename Dest> void evalTo(Dest& dst) const
ykuroda 0:13a5d365ba16 452 {
ykuroda 0:13a5d365ba16 453 /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
ykuroda 0:13a5d365ba16 454 * So we proceed as follows:
ykuroda 0:13a5d365ba16 455 * Step 1: compute c = Pb.
ykuroda 0:13a5d365ba16 456 * Step 2: replace c by the solution x to Lx = c.
ykuroda 0:13a5d365ba16 457 * Step 3: replace c by the solution x to Ux = c.
ykuroda 0:13a5d365ba16 458 */
ykuroda 0:13a5d365ba16 459
ykuroda 0:13a5d365ba16 460 eigen_assert(rhs().rows() == dec().matrixLU().rows());
ykuroda 0:13a5d365ba16 461
ykuroda 0:13a5d365ba16 462 // Step 1
ykuroda 0:13a5d365ba16 463 dst = dec().permutationP() * rhs();
ykuroda 0:13a5d365ba16 464
ykuroda 0:13a5d365ba16 465 // Step 2
ykuroda 0:13a5d365ba16 466 dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst);
ykuroda 0:13a5d365ba16 467
ykuroda 0:13a5d365ba16 468 // Step 3
ykuroda 0:13a5d365ba16 469 dec().matrixLU().template triangularView<Upper>().solveInPlace(dst);
ykuroda 0:13a5d365ba16 470 }
ykuroda 0:13a5d365ba16 471 };
ykuroda 0:13a5d365ba16 472
ykuroda 0:13a5d365ba16 473 } // end namespace internal
ykuroda 0:13a5d365ba16 474
ykuroda 0:13a5d365ba16 475 /******** MatrixBase methods *******/
ykuroda 0:13a5d365ba16 476
ykuroda 0:13a5d365ba16 477 /** \lu_module
ykuroda 0:13a5d365ba16 478 *
ykuroda 0:13a5d365ba16 479 * \return the partial-pivoting LU decomposition of \c *this.
ykuroda 0:13a5d365ba16 480 *
ykuroda 0:13a5d365ba16 481 * \sa class PartialPivLU
ykuroda 0:13a5d365ba16 482 */
ykuroda 0:13a5d365ba16 483 template<typename Derived>
ykuroda 0:13a5d365ba16 484 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
ykuroda 0:13a5d365ba16 485 MatrixBase<Derived>::partialPivLu() const
ykuroda 0:13a5d365ba16 486 {
ykuroda 0:13a5d365ba16 487 return PartialPivLU<PlainObject>(eval());
ykuroda 0:13a5d365ba16 488 }
ykuroda 0:13a5d365ba16 489
ykuroda 0:13a5d365ba16 490 #if EIGEN2_SUPPORT_STAGE > STAGE20_RESOLVE_API_CONFLICTS
ykuroda 0:13a5d365ba16 491 /** \lu_module
ykuroda 0:13a5d365ba16 492 *
ykuroda 0:13a5d365ba16 493 * Synonym of partialPivLu().
ykuroda 0:13a5d365ba16 494 *
ykuroda 0:13a5d365ba16 495 * \return the partial-pivoting LU decomposition of \c *this.
ykuroda 0:13a5d365ba16 496 *
ykuroda 0:13a5d365ba16 497 * \sa class PartialPivLU
ykuroda 0:13a5d365ba16 498 */
ykuroda 0:13a5d365ba16 499 template<typename Derived>
ykuroda 0:13a5d365ba16 500 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
ykuroda 0:13a5d365ba16 501 MatrixBase<Derived>::lu() const
ykuroda 0:13a5d365ba16 502 {
ykuroda 0:13a5d365ba16 503 return PartialPivLU<PlainObject>(eval());
ykuroda 0:13a5d365ba16 504 }
ykuroda 0:13a5d365ba16 505 #endif
ykuroda 0:13a5d365ba16 506
ykuroda 0:13a5d365ba16 507 } // end namespace Eigen
ykuroda 0:13a5d365ba16 508
ykuroda 0:13a5d365ba16 509 #endif // EIGEN_PARTIALLU_H