PartialPivLU.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
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   // FIXME add a stride to Map, so that the following mapping becomes easier,
00234   // another option would be to create an expression being able to automatically
00235   // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly
00236   // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix,
00237   // and Block.
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         // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k)
00282         // overflow but not the actual quotient?
00283         lu.col(k).tail(rrows) /= lu.coeff(k,k);
00284       }
00285       else if(first_zero_pivot==-1)
00286       {
00287         // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
00288         // and continue the factorization such we still have A = PLU
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     // if the matrix is too small, no blocking:
00321     if(size<=16)
00322     {
00323       return unblocked_lu(lu, row_transpositions, nb_transpositions);
00324     }
00325 
00326     // automatically adjust the number of subdivisions to the size
00327     // of the matrix so that there is enough sub blocks:
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); // actual size of the block
00340       Index trows = rows - k - bs; // trailing rows
00341       Index tsize = size - k - bs; // trailing size
00342 
00343       // partition the matrix:
00344       //                          A00 | A01 | A02
00345       // lu  = A_0 | A_1 | A_2 =  A10 | A11 | A12
00346       //                          A20 | A21 | A22
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       // recursively call the blocked LU algorithm on [A11^T A21^T]^T
00356       // with a very small blocking size:
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       // update permutations and apply them to A_0
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         // apply permutations to A_2
00373         for(Index i=k;i<k+bs; ++i)
00374           A_2.row(i).swap(A_2.row(row_transpositions[i]));
00375 
00376         // A12 = A11^-1 A12
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 } // end namespace internal
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   // LU
00436   MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
00437                  * m_lu.template triangularView<Upper>();
00438 
00439   // P^{-1}(LU)
00440   res = m_p.inverse() * res;
00441 
00442   return res;
00443 }
00444 
00445 /***** Implementation of solve() *****************************************************/
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     /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
00458     * So we proceed as follows:
00459     * Step 1: compute c = Pb.
00460     * Step 2: replace c by the solution x to Lx = c.
00461     * Step 3: replace c by the solution x to Ux = c.
00462     */
00463 
00464     eigen_assert(rhs().rows() == dec().matrixLU().rows());
00465 
00466     // Step 1
00467     dst = dec().permutationP() * rhs();
00468 
00469     // Step 2
00470     dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst);
00471 
00472     // Step 3
00473     dec().matrixLU().template triangularView<Upper>().solveInPlace(dst);
00474   }
00475 };
00476 
00477 } // end namespace internal
00478 
00479 /******** MatrixBase methods *******/
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 } // end namespace Eigen
00512 
00513 #endif // EIGEN_PARTIALLU_H