00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef EIGEN_PARTIALLU_H
00027 #define EIGEN_PARTIALLU_H
00028
00029 namespace Eigen {
00030
00062 template<typename _MatrixType> class PartialPivLU
00063 {
00064 public:
00065
00066 typedef _MatrixType MatrixType;
00067 enum {
00068 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00069 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00070 Options = MatrixType::Options,
00071 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00072 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00073 };
00074 typedef typename MatrixType::Scalar Scalar;
00075 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00076 typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
00077 typedef typename MatrixType::Index Index;
00078 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
00079 typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
00080
00081
00088 PartialPivLU();
00089
00096 PartialPivLU(Index size);
00097
00105 PartialPivLU(const MatrixType& matrix);
00106
00107 PartialPivLU& compute(const MatrixType& matrix);
00108
00115 inline const MatrixType& matrixLU() const
00116 {
00117 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00118 return m_lu;
00119 }
00120
00123 inline const PermutationType& permutationP() const
00124 {
00125 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00126 return m_p;
00127 }
00128
00146 template<typename Rhs>
00147 inline const internal::solve_retval<PartialPivLU, Rhs>
00148 solve(const MatrixBase<Rhs>& b) const
00149 {
00150 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00151 return internal::solve_retval<PartialPivLU, Rhs>(*this, b.derived());
00152 }
00153
00161 inline const internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> inverse() const
00162 {
00163 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00164 return internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType>
00165 (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
00166 }
00167
00181 typename internal::traits<MatrixType>::Scalar determinant() const;
00182
00183 MatrixType reconstructedMatrix() const;
00184
00185 inline Index rows() const { return m_lu.rows(); }
00186 inline Index cols() const { return m_lu.cols(); }
00187
00188 protected:
00189 MatrixType m_lu;
00190 PermutationType m_p;
00191 TranspositionType m_rowsTranspositions;
00192 Index m_det_p;
00193 bool m_isInitialized;
00194 };
00195
00196 template<typename MatrixType>
00197 PartialPivLU<MatrixType>::PartialPivLU()
00198 : m_lu(),
00199 m_p(),
00200 m_rowsTranspositions(),
00201 m_det_p(0),
00202 m_isInitialized(false)
00203 {
00204 }
00205
00206 template<typename MatrixType>
00207 PartialPivLU<MatrixType>::PartialPivLU(Index size)
00208 : m_lu(size, size),
00209 m_p(size),
00210 m_rowsTranspositions(size),
00211 m_det_p(0),
00212 m_isInitialized(false)
00213 {
00214 }
00215
00216 template<typename MatrixType>
00217 PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
00218 : m_lu(matrix.rows(), matrix.rows()),
00219 m_p(matrix.rows()),
00220 m_rowsTranspositions(matrix.rows()),
00221 m_det_p(0),
00222 m_isInitialized(false)
00223 {
00224 compute(matrix);
00225 }
00226
00227 namespace internal {
00228
00230 template<typename Scalar, int StorageOrder, typename PivIndex>
00231 struct partial_lu_impl
00232 {
00233
00234
00235
00236
00237
00238 typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > MapLU;
00239 typedef Block<MapLU, Dynamic, Dynamic> MatrixType;
00240 typedef Block<MatrixType,Dynamic,Dynamic> BlockType;
00241 typedef typename MatrixType::RealScalar RealScalar;
00242 typedef typename MatrixType::Index Index;
00243
00254 static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
00255 {
00256 const Index rows = lu.rows();
00257 const Index cols = lu.cols();
00258 const Index size = (std::min)(rows,cols);
00259 nb_transpositions = 0;
00260 int first_zero_pivot = -1;
00261 for(Index k = 0; k < size; ++k)
00262 {
00263 Index rrows = rows-k-1;
00264 Index rcols = cols-k-1;
00265
00266 Index row_of_biggest_in_col;
00267 RealScalar biggest_in_corner
00268 = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col);
00269 row_of_biggest_in_col += k;
00270
00271 row_transpositions[k] = row_of_biggest_in_col;
00272
00273 if(biggest_in_corner != RealScalar(0))
00274 {
00275 if(k != row_of_biggest_in_col)
00276 {
00277 lu.row(k).swap(lu.row(row_of_biggest_in_col));
00278 ++nb_transpositions;
00279 }
00280
00281
00282
00283 lu.col(k).tail(rrows) /= lu.coeff(k,k);
00284 }
00285 else if(first_zero_pivot==-1)
00286 {
00287
00288
00289 first_zero_pivot = k;
00290 }
00291
00292 if(k<rows-1)
00293 lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
00294 }
00295 return first_zero_pivot;
00296 }
00297
00313 static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
00314 {
00315 MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
00316 MatrixType lu(lu1,0,0,rows,cols);
00317
00318 const Index size = (std::min)(rows,cols);
00319
00320
00321 if(size<=16)
00322 {
00323 return unblocked_lu(lu, row_transpositions, nb_transpositions);
00324 }
00325
00326
00327
00328 Index blockSize;
00329 {
00330 blockSize = size/8;
00331 blockSize = (blockSize/16)*16;
00332 blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
00333 }
00334
00335 nb_transpositions = 0;
00336 int first_zero_pivot = -1;
00337 for(Index k = 0; k < size; k+=blockSize)
00338 {
00339 Index bs = (std::min)(size-k,blockSize);
00340 Index trows = rows - k - bs;
00341 Index tsize = size - k - bs;
00342
00343
00344
00345
00346
00347 BlockType A_0(lu,0,0,rows,k);
00348 BlockType A_2(lu,0,k+bs,rows,tsize);
00349 BlockType A11(lu,k,k,bs,bs);
00350 BlockType A12(lu,k,k+bs,bs,tsize);
00351 BlockType A21(lu,k+bs,k,trows,bs);
00352 BlockType A22(lu,k+bs,k+bs,trows,tsize);
00353
00354 PivIndex nb_transpositions_in_panel;
00355
00356
00357 Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
00358 row_transpositions+k, nb_transpositions_in_panel, 16);
00359 if(ret>=0 && first_zero_pivot==-1)
00360 first_zero_pivot = k+ret;
00361
00362 nb_transpositions += nb_transpositions_in_panel;
00363
00364 for(Index i=k; i<k+bs; ++i)
00365 {
00366 Index piv = (row_transpositions[i] += k);
00367 A_0.row(i).swap(A_0.row(piv));
00368 }
00369
00370 if(trows)
00371 {
00372
00373 for(Index i=k;i<k+bs; ++i)
00374 A_2.row(i).swap(A_2.row(row_transpositions[i]));
00375
00376
00377 A11.template triangularView<UnitLower>().solveInPlace(A12);
00378
00379 A22.noalias() -= A21 * A12;
00380 }
00381 }
00382 return first_zero_pivot;
00383 }
00384 };
00385
00388 template<typename MatrixType, typename TranspositionType>
00389 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::Index& nb_transpositions)
00390 {
00391 eigen_assert(lu.cols() == row_transpositions.size());
00392 eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
00393
00394 partial_lu_impl
00395 <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::Index>
00396 ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
00397 }
00398
00399 }
00400
00401 template<typename MatrixType>
00402 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix)
00403 {
00404 m_lu = matrix;
00405
00406 eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
00407 const Index size = matrix.rows();
00408
00409 m_rowsTranspositions.resize(size);
00410
00411 typename TranspositionType::Index nb_transpositions;
00412 internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
00413 m_det_p = (nb_transpositions%2) ? -1 : 1;
00414
00415 m_p = m_rowsTranspositions;
00416
00417 m_isInitialized = true;
00418 return *this;
00419 }
00420
00421 template<typename MatrixType>
00422 typename internal::traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
00423 {
00424 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00425 return Scalar(m_det_p) * m_lu.diagonal().prod();
00426 }
00427
00431 template<typename MatrixType>
00432 MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
00433 {
00434 eigen_assert(m_isInitialized && "LU is not initialized.");
00435
00436 MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
00437 * m_lu.template triangularView<Upper>();
00438
00439
00440 res = m_p.inverse() * res;
00441
00442 return res;
00443 }
00444
00445
00446
00447 namespace internal {
00448
00449 template<typename _MatrixType, typename Rhs>
00450 struct solve_retval<PartialPivLU<_MatrixType>, Rhs>
00451 : solve_retval_base<PartialPivLU<_MatrixType>, Rhs>
00452 {
00453 EIGEN_MAKE_SOLVE_HELPERS(PartialPivLU<_MatrixType>,Rhs)
00454
00455 template<typename Dest> void evalTo(Dest& dst) const
00456 {
00457
00458
00459
00460
00461
00462
00463
00464 eigen_assert(rhs().rows() == dec().matrixLU().rows());
00465
00466
00467 dst = dec().permutationP() * rhs();
00468
00469
00470 dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst);
00471
00472
00473 dec().matrixLU().template triangularView<Upper>().solveInPlace(dst);
00474 }
00475 };
00476
00477 }
00478
00479
00480
00487 template<typename Derived>
00488 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
00489 MatrixBase<Derived>::partialPivLu() const
00490 {
00491 return PartialPivLU<PlainObject>(eval());
00492 }
00493
00494 #if EIGEN2_SUPPORT_STAGE > STAGE20_RESOLVE_API_CONFLICTS
00495
00503 template<typename Derived>
00504 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
00505 MatrixBase<Derived>::lu() const
00506 {
00507 return PartialPivLU<PlainObject>(eval());
00508 }
00509 #endif
00510
00511 }
00512
00513 #endif // EIGEN_PARTIALLU_H