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
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060 template<typename _MatrixType> class PartialPivLU
00061 {
00062 public:
00063
00064 typedef _MatrixType MatrixType;
00065 enum {
00066 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00067 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00068 Options = MatrixType::Options,
00069 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00070 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00071 };
00072 typedef typename MatrixType::Scalar Scalar;
00073 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00074 typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
00075 typedef typename MatrixType::Index Index;
00076 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
00077 typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
00078
00079
00080
00081
00082
00083
00084
00085
00086 PartialPivLU();
00087
00088
00089
00090
00091
00092
00093
00094 PartialPivLU(Index size);
00095
00096
00097
00098
00099
00100
00101
00102
00103 PartialPivLU(const MatrixType& matrix);
00104
00105 PartialPivLU& compute(const MatrixType& matrix);
00106
00107
00108
00109
00110
00111
00112
00113 inline const MatrixType& matrixLU() const
00114 {
00115 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00116 return m_lu;
00117 }
00118
00119
00120
00121 inline const PermutationType& permutationP() const
00122 {
00123 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00124 return m_p;
00125 }
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144 template<typename Rhs>
00145 inline const internal::solve_retval<PartialPivLU, Rhs>
00146 solve(const MatrixBase<Rhs>& b) const
00147 {
00148 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00149 return internal::solve_retval<PartialPivLU, Rhs>(*this, b.derived());
00150 }
00151
00152
00153
00154
00155
00156
00157
00158
00159 inline const internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> inverse() const
00160 {
00161 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00162 return internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType>
00163 (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
00164 }
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179 typename internal::traits<MatrixType>::Scalar determinant() const;
00180
00181 MatrixType reconstructedMatrix() const;
00182
00183 inline Index rows() const { return m_lu.rows(); }
00184 inline Index cols() const { return m_lu.cols(); }
00185
00186 protected:
00187 MatrixType m_lu;
00188 PermutationType m_p;
00189 TranspositionType m_rowsTranspositions;
00190 Index m_det_p;
00191 bool m_isInitialized;
00192 };
00193
00194 template<typename MatrixType>
00195 PartialPivLU<MatrixType>::PartialPivLU()
00196 : m_lu(),
00197 m_p(),
00198 m_rowsTranspositions(),
00199 m_det_p(0),
00200 m_isInitialized(false)
00201 {
00202 }
00203
00204 template<typename MatrixType>
00205 PartialPivLU<MatrixType>::PartialPivLU(Index size)
00206 : m_lu(size, size),
00207 m_p(size),
00208 m_rowsTranspositions(size),
00209 m_det_p(0),
00210 m_isInitialized(false)
00211 {
00212 }
00213
00214 template<typename MatrixType>
00215 PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
00216 : m_lu(matrix.rows(), matrix.rows()),
00217 m_p(matrix.rows()),
00218 m_rowsTranspositions(matrix.rows()),
00219 m_det_p(0),
00220 m_isInitialized(false)
00221 {
00222 compute(matrix);
00223 }
00224
00225 namespace internal {
00226
00227
00228 template<typename Scalar, int StorageOrder>
00229 struct partial_lu_impl
00230 {
00231
00232
00233
00234
00235
00236 typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > MapLU;
00237 typedef Block<MapLU, Dynamic, Dynamic> MatrixType;
00238 typedef Block<MatrixType,Dynamic,Dynamic> BlockType;
00239 typedef typename MatrixType::RealScalar RealScalar;
00240 typedef typename MatrixType::Index Index;
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254 static bool unblocked_lu(MatrixType& lu, Index* row_transpositions, Index& nb_transpositions)
00255 {
00256 const Index rows = lu.rows();
00257 const Index size = std::min(lu.rows(),lu.cols());
00258 nb_transpositions = 0;
00259 for(Index k = 0; k < size; ++k)
00260 {
00261 Index row_of_biggest_in_col;
00262 RealScalar biggest_in_corner
00263 = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col);
00264 row_of_biggest_in_col += k;
00265
00266 if(biggest_in_corner == 0)
00267 {
00268
00269
00270
00271
00272
00273 for(Index i = k; i < size; i++)
00274 row_transpositions[i] = i;
00275 return false;
00276 }
00277
00278 row_transpositions[k] = row_of_biggest_in_col;
00279
00280 if(k != row_of_biggest_in_col)
00281 {
00282 lu.row(k).swap(lu.row(row_of_biggest_in_col));
00283 ++nb_transpositions;
00284 }
00285
00286 if(k<rows-1)
00287 {
00288 Index rrows = rows-k-1;
00289 Index rsize = size-k-1;
00290 lu.col(k).tail(rrows) /= lu.coeff(k,k);
00291 lu.bottomRightCorner(rrows,rsize).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rsize);
00292 }
00293 }
00294 return true;
00295 }
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314 static bool blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, Index* row_transpositions, Index& nb_transpositions, Index maxBlockSize=256)
00315 {
00316 MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
00317 MatrixType lu(lu1,0,0,rows,cols);
00318
00319 const Index size = std::min(rows,cols);
00320
00321
00322 if(size<=16)
00323 {
00324 return unblocked_lu(lu, row_transpositions, nb_transpositions);
00325 }
00326
00327
00328
00329 Index blockSize;
00330 {
00331 blockSize = size/8;
00332 blockSize = (blockSize/16)*16;
00333 blockSize = std::min(std::max(blockSize,Index(8)), maxBlockSize);
00334 }
00335
00336 nb_transpositions = 0;
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 Index nb_transpositions_in_panel;
00355
00356
00357 if(!blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
00358 row_transpositions+k, nb_transpositions_in_panel, 16))
00359 {
00360
00361
00362
00363 for(Index i=k; i<size; ++i)
00364 row_transpositions[i] = i;
00365 return false;
00366 }
00367 nb_transpositions += nb_transpositions_in_panel;
00368
00369
00370 for(Index i=k; i<k+bs; ++i)
00371 {
00372 Index piv = (row_transpositions[i] += k);
00373 A_0.row(i).swap(A_0.row(piv));
00374 }
00375
00376 if(trows)
00377 {
00378
00379 for(Index i=k;i<k+bs; ++i)
00380 A_2.row(i).swap(A_2.row(row_transpositions[i]));
00381
00382
00383 A11.template triangularView<UnitLower>().solveInPlace(A12);
00384
00385 A22.noalias() -= A21 * A12;
00386 }
00387 }
00388 return true;
00389 }
00390 };
00391
00392
00393
00394 template<typename MatrixType, typename TranspositionType>
00395 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename MatrixType::Index& nb_transpositions)
00396 {
00397 eigen_assert(lu.cols() == row_transpositions.size());
00398 eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
00399
00400 partial_lu_impl
00401 <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor>
00402 ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
00403 }
00404
00405 }
00406
00407 template<typename MatrixType>
00408 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix)
00409 {
00410 m_lu = matrix;
00411
00412 eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
00413 const Index size = matrix.rows();
00414
00415 m_rowsTranspositions.resize(size);
00416
00417 Index nb_transpositions;
00418 internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
00419 m_det_p = (nb_transpositions%2) ? -1 : 1;
00420
00421 m_p = m_rowsTranspositions;
00422
00423 m_isInitialized = true;
00424 return *this;
00425 }
00426
00427 template<typename MatrixType>
00428 typename internal::traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
00429 {
00430 eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00431 return Scalar(m_det_p) * m_lu.diagonal().prod();
00432 }
00433
00434
00435
00436
00437 template<typename MatrixType>
00438 MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
00439 {
00440 eigen_assert(m_isInitialized && "LU is not initialized.");
00441
00442 MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
00443 * m_lu.template triangularView<Upper>();
00444
00445
00446 res = m_p.inverse() * res;
00447
00448 return res;
00449 }
00450
00451
00452
00453 namespace internal {
00454
00455 template<typename _MatrixType, typename Rhs>
00456 struct solve_retval<PartialPivLU<_MatrixType>, Rhs>
00457 : solve_retval_base<PartialPivLU<_MatrixType>, Rhs>
00458 {
00459 EIGEN_MAKE_SOLVE_HELPERS(PartialPivLU<_MatrixType>,Rhs)
00460
00461 template<typename Dest> void evalTo(Dest& dst) const
00462 {
00463
00464
00465
00466
00467
00468
00469
00470 eigen_assert(rhs().rows() == dec().matrixLU().rows());
00471
00472
00473 dst = dec().permutationP() * rhs();
00474
00475
00476 dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst);
00477
00478
00479 dec().matrixLU().template triangularView<Upper>().solveInPlace(dst);
00480 }
00481 };
00482
00483 }
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493 template<typename Derived>
00494 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
00495 MatrixBase<Derived>::partialPivLu() const
00496 {
00497 return PartialPivLU<PlainObject>(eval());
00498 }
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508 template<typename Derived>
00509 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
00510 MatrixBase<Derived>::lu() const
00511 {
00512 return PartialPivLU<PlainObject>(eval());
00513 }
00514
00515 #endif // EIGEN_PARTIALLU_H