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 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_PARTIALLU_H
00012 #define EIGEN_PARTIALLU_H
00013 
00014 namespace Eigen { 
00015 
00047 template<typename _MatrixType> class PartialPivLU
00048 {
00049   public:
00050 
00051     typedef _MatrixType MatrixType;
00052     enum {
00053       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00054       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00055       Options = MatrixType::Options,
00056       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00057       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00058     };
00059     typedef typename MatrixType::Scalar Scalar;
00060     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00061     typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
00062     typedef typename MatrixType::Index Index;
00063     typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
00064     typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
00065 
00066 
00073     PartialPivLU();
00074 
00081     PartialPivLU(Index size);
00082 
00090     PartialPivLU(const MatrixType& matrix);
00091 
00092     PartialPivLU& compute(const MatrixType& matrix);
00093 
00100     inline const MatrixType& matrixLU() const
00101     {
00102       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00103       return m_lu;
00104     }
00105 
00108     inline const PermutationType& permutationP() const
00109     {
00110       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00111       return m_p;
00112     }
00113 
00131     template<typename Rhs>
00132     inline const internal::solve_retval<PartialPivLU, Rhs>
00133     solve(const MatrixBase<Rhs>& b) const
00134     {
00135       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00136       return internal::solve_retval<PartialPivLU, Rhs>(*this, b.derived());
00137     }
00138 
00146     inline const internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> inverse() const
00147     {
00148       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00149       return internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType>
00150                (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
00151     }
00152 
00166     typename internal::traits<MatrixType>::Scalar determinant() const;
00167 
00168     MatrixType reconstructedMatrix() const;
00169 
00170     inline Index rows() const { return m_lu.rows(); }
00171     inline Index cols() const { return m_lu.cols(); }
00172 
00173   protected:
00174     MatrixType m_lu;
00175     PermutationType m_p;
00176     TranspositionType m_rowsTranspositions;
00177     Index m_det_p;
00178     bool m_isInitialized;
00179 };
00180 
00181 template<typename MatrixType>
00182 PartialPivLU<MatrixType>::PartialPivLU()
00183   : m_lu(),
00184     m_p(),
00185     m_rowsTranspositions(),
00186     m_det_p(0),
00187     m_isInitialized(false)
00188 {
00189 }
00190 
00191 template<typename MatrixType>
00192 PartialPivLU<MatrixType>::PartialPivLU(Index size)
00193   : m_lu(size, size),
00194     m_p(size),
00195     m_rowsTranspositions(size),
00196     m_det_p(0),
00197     m_isInitialized(false)
00198 {
00199 }
00200 
00201 template<typename MatrixType>
00202 PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
00203   : m_lu(matrix.rows(), matrix.rows()),
00204     m_p(matrix.rows()),
00205     m_rowsTranspositions(matrix.rows()),
00206     m_det_p(0),
00207     m_isInitialized(false)
00208 {
00209   compute(matrix);
00210 }
00211 
00212 namespace internal {
00213 
00215 template<typename Scalar, int StorageOrder, typename PivIndex>
00216 struct partial_lu_impl
00217 {
00218   // FIXME add a stride to Map, so that the following mapping becomes easier,
00219   // another option would be to create an expression being able to automatically
00220   // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly
00221   // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix,
00222   // and Block.
00223   typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > MapLU;
00224   typedef Block<MapLU, Dynamic, Dynamic> MatrixType;
00225   typedef Block<MatrixType,Dynamic,Dynamic> BlockType;
00226   typedef typename MatrixType::RealScalar RealScalar;
00227   typedef typename MatrixType::Index Index;
00228 
00239   static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
00240   {
00241     const Index rows = lu.rows();
00242     const Index cols = lu.cols();
00243     const Index size = (std::min)(rows,cols);
00244     nb_transpositions = 0;
00245     int first_zero_pivot = -1;
00246     for(Index k = 0; k < size; ++k)
00247     {
00248       Index rrows = rows-k-1;
00249       Index rcols = cols-k-1;
00250         
00251       Index row_of_biggest_in_col;
00252       RealScalar biggest_in_corner
00253         = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col);
00254       row_of_biggest_in_col += k;
00255 
00256       row_transpositions[k] = row_of_biggest_in_col;
00257 
00258       if(biggest_in_corner != RealScalar(0))
00259       {
00260         if(k != row_of_biggest_in_col)
00261         {
00262           lu.row(k).swap(lu.row(row_of_biggest_in_col));
00263           ++nb_transpositions;
00264         }
00265 
00266         // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k)
00267         // overflow but not the actual quotient?
00268         lu.col(k).tail(rrows) /= lu.coeff(k,k);
00269       }
00270       else if(first_zero_pivot==-1)
00271       {
00272         // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
00273         // and continue the factorization such we still have A = PLU
00274         first_zero_pivot = k;
00275       }
00276 
00277       if(k<rows-1)
00278         lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
00279     }
00280     return first_zero_pivot;
00281   }
00282 
00298   static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
00299   {
00300     MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
00301     MatrixType lu(lu1,0,0,rows,cols);
00302 
00303     const Index size = (std::min)(rows,cols);
00304 
00305     // if the matrix is too small, no blocking:
00306     if(size<=16)
00307     {
00308       return unblocked_lu(lu, row_transpositions, nb_transpositions);
00309     }
00310 
00311     // automatically adjust the number of subdivisions to the size
00312     // of the matrix so that there is enough sub blocks:
00313     Index blockSize;
00314     {
00315       blockSize = size/8;
00316       blockSize = (blockSize/16)*16;
00317       blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
00318     }
00319 
00320     nb_transpositions = 0;
00321     int first_zero_pivot = -1;
00322     for(Index k = 0; k < size; k+=blockSize)
00323     {
00324       Index bs = (std::min)(size-k,blockSize); // actual size of the block
00325       Index trows = rows - k - bs; // trailing rows
00326       Index tsize = size - k - bs; // trailing size
00327 
00328       // partition the matrix:
00329       //                          A00 | A01 | A02
00330       // lu  = A_0 | A_1 | A_2 =  A10 | A11 | A12
00331       //                          A20 | A21 | A22
00332       BlockType A_0(lu,0,0,rows,k);
00333       BlockType A_2(lu,0,k+bs,rows,tsize);
00334       BlockType A11(lu,k,k,bs,bs);
00335       BlockType A12(lu,k,k+bs,bs,tsize);
00336       BlockType A21(lu,k+bs,k,trows,bs);
00337       BlockType A22(lu,k+bs,k+bs,trows,tsize);
00338 
00339       PivIndex nb_transpositions_in_panel;
00340       // recursively call the blocked LU algorithm on [A11^T A21^T]^T
00341       // with a very small blocking size:
00342       Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
00343                    row_transpositions+k, nb_transpositions_in_panel, 16);
00344       if(ret>=0 && first_zero_pivot==-1)
00345         first_zero_pivot = k+ret;
00346 
00347       nb_transpositions += nb_transpositions_in_panel;
00348       // update permutations and apply them to A_0
00349       for(Index i=k; i<k+bs; ++i)
00350       {
00351         Index piv = (row_transpositions[i] += k);
00352         A_0.row(i).swap(A_0.row(piv));
00353       }
00354 
00355       if(trows)
00356       {
00357         // apply permutations to A_2
00358         for(Index i=k;i<k+bs; ++i)
00359           A_2.row(i).swap(A_2.row(row_transpositions[i]));
00360 
00361         // A12 = A11^-1 A12
00362         A11.template triangularView<UnitLower>().solveInPlace(A12);
00363 
00364         A22.noalias() -= A21 * A12;
00365       }
00366     }
00367     return first_zero_pivot;
00368   }
00369 };
00370 
00373 template<typename MatrixType, typename TranspositionType>
00374 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::Index& nb_transpositions)
00375 {
00376   eigen_assert(lu.cols() == row_transpositions.size());
00377   eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
00378 
00379   partial_lu_impl
00380     <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::Index>
00381     ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
00382 }
00383 
00384 } // end namespace internal
00385 
00386 template<typename MatrixType>
00387 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix)
00388 {
00389   m_lu = matrix;
00390 
00391   eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
00392   const Index size = matrix.rows();
00393 
00394   m_rowsTranspositions.resize(size);
00395 
00396   typename TranspositionType::Index nb_transpositions;
00397   internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
00398   m_det_p = (nb_transpositions%2) ? -1 : 1;
00399 
00400   m_p = m_rowsTranspositions;
00401 
00402   m_isInitialized = true;
00403   return *this;
00404 }
00405 
00406 template<typename MatrixType>
00407 typename internal::traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
00408 {
00409   eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00410   return Scalar(m_det_p) * m_lu.diagonal().prod();
00411 }
00412 
00416 template<typename MatrixType>
00417 MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
00418 {
00419   eigen_assert(m_isInitialized && "LU is not initialized.");
00420   // LU
00421   MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
00422                  * m_lu.template triangularView<Upper>();
00423 
00424   // P^{-1}(LU)
00425   res = m_p.inverse() * res;
00426 
00427   return res;
00428 }
00429 
00430 /***** Implementation of solve() *****************************************************/
00431 
00432 namespace internal {
00433 
00434 template<typename _MatrixType, typename Rhs>
00435 struct solve_retval<PartialPivLU<_MatrixType>, Rhs>
00436   : solve_retval_base<PartialPivLU<_MatrixType>, Rhs>
00437 {
00438   EIGEN_MAKE_SOLVE_HELPERS(PartialPivLU<_MatrixType>,Rhs)
00439 
00440   template<typename Dest> void evalTo(Dest& dst) const
00441   {
00442     /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
00443     * So we proceed as follows:
00444     * Step 1: compute c = Pb.
00445     * Step 2: replace c by the solution x to Lx = c.
00446     * Step 3: replace c by the solution x to Ux = c.
00447     */
00448 
00449     eigen_assert(rhs().rows() == dec().matrixLU().rows());
00450 
00451     // Step 1
00452     dst = dec().permutationP() * rhs();
00453 
00454     // Step 2
00455     dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst);
00456 
00457     // Step 3
00458     dec().matrixLU().template triangularView<Upper>().solveInPlace(dst);
00459   }
00460 };
00461 
00462 } // end namespace internal
00463 
00464 /******** MatrixBase methods *******/
00465 
00472 template<typename Derived>
00473 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
00474 MatrixBase<Derived>::partialPivLu() const
00475 {
00476   return PartialPivLU<PlainObject>(eval());
00477 }
00478 
00479 #if EIGEN2_SUPPORT_STAGE > STAGE20_RESOLVE_API_CONFLICTS
00480 
00488 template<typename Derived>
00489 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
00490 MatrixBase<Derived>::lu() const
00491 {
00492   return PartialPivLU<PlainObject>(eval());
00493 }
00494 #endif
00495 
00496 } // end namespace Eigen
00497 
00498 #endif // EIGEN_PARTIALLU_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:25:39