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 
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 
00086     PartialPivLU();
00087 
00094     PartialPivLU(Index size);
00095 
00103     PartialPivLU(const MatrixType& matrix);
00104 
00105     PartialPivLU& compute(const MatrixType& matrix);
00106 
00113     inline const MatrixType& matrixLU() const
00114     {
00115       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00116       return m_lu;
00117     }
00118 
00121     inline const PermutationType& permutationP() const
00122     {
00123       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00124       return m_p;
00125     }
00126 
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 
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 
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 
00228 template<typename Scalar, int StorageOrder, typename PivIndex>
00229 struct partial_lu_impl
00230 {
00231   // FIXME add a stride to Map, so that the following mapping becomes easier,
00232   // another option would be to create an expression being able to automatically
00233   // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly
00234   // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix,
00235   // and Block.
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 
00252   static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
00253   {
00254     const Index rows = lu.rows();
00255     const Index cols = lu.cols();
00256     const Index size = (std::min)(rows,cols);
00257     nb_transpositions = 0;
00258     int first_zero_pivot = -1;
00259     for(Index k = 0; k < size; ++k)
00260     {
00261       Index rrows = rows-k-1;
00262       Index rcols = cols-k-1;
00263         
00264       Index row_of_biggest_in_col;
00265       RealScalar biggest_in_corner
00266         = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col);
00267       row_of_biggest_in_col += k;
00268 
00269       row_transpositions[k] = row_of_biggest_in_col;
00270 
00271       if(biggest_in_corner != RealScalar(0))
00272       {
00273         if(k != row_of_biggest_in_col)
00274         {
00275           lu.row(k).swap(lu.row(row_of_biggest_in_col));
00276           ++nb_transpositions;
00277         }
00278 
00279         // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k)
00280         // overflow but not the actual quotient?
00281         lu.col(k).tail(rrows) /= lu.coeff(k,k);
00282       }
00283       else if(first_zero_pivot==-1)
00284       {
00285         // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
00286         // and continue the factorization such we still have A = PLU
00287         first_zero_pivot = k;
00288       }
00289 
00290       if(k<rows-1)
00291         lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
00292     }
00293     return first_zero_pivot;
00294   }
00295 
00311   static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
00312   {
00313     MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
00314     MatrixType lu(lu1,0,0,rows,cols);
00315 
00316     const Index size = (std::min)(rows,cols);
00317 
00318     // if the matrix is too small, no blocking:
00319     if(size<=16)
00320     {
00321       return unblocked_lu(lu, row_transpositions, nb_transpositions);
00322     }
00323 
00324     // automatically adjust the number of subdivisions to the size
00325     // of the matrix so that there is enough sub blocks:
00326     Index blockSize;
00327     {
00328       blockSize = size/8;
00329       blockSize = (blockSize/16)*16;
00330       blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
00331     }
00332 
00333     nb_transpositions = 0;
00334     int first_zero_pivot = -1;
00335     for(Index k = 0; k < size; k+=blockSize)
00336     {
00337       Index bs = (std::min)(size-k,blockSize); // actual size of the block
00338       Index trows = rows - k - bs; // trailing rows
00339       Index tsize = size - k - bs; // trailing size
00340 
00341       // partition the matrix:
00342       //                          A00 | A01 | A02
00343       // lu  = A_0 | A_1 | A_2 =  A10 | A11 | A12
00344       //                          A20 | A21 | A22
00345       BlockType A_0(lu,0,0,rows,k);
00346       BlockType A_2(lu,0,k+bs,rows,tsize);
00347       BlockType A11(lu,k,k,bs,bs);
00348       BlockType A12(lu,k,k+bs,bs,tsize);
00349       BlockType A21(lu,k+bs,k,trows,bs);
00350       BlockType A22(lu,k+bs,k+bs,trows,tsize);
00351 
00352       PivIndex nb_transpositions_in_panel;
00353       // recursively call the blocked LU algorithm on [A11^T A21^T]^T
00354       // with a very small blocking size:
00355       Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
00356                    row_transpositions+k, nb_transpositions_in_panel, 16);
00357       if(ret>=0 && first_zero_pivot==-1)
00358         first_zero_pivot = k+ret;
00359 
00360       nb_transpositions += nb_transpositions_in_panel;
00361       // update permutations and apply them to A_0
00362       for(Index i=k; i<k+bs; ++i)
00363       {
00364         Index piv = (row_transpositions[i] += k);
00365         A_0.row(i).swap(A_0.row(piv));
00366       }
00367 
00368       if(trows)
00369       {
00370         // apply permutations to A_2
00371         for(Index i=k;i<k+bs; ++i)
00372           A_2.row(i).swap(A_2.row(row_transpositions[i]));
00373 
00374         // A12 = A11^-1 A12
00375         A11.template triangularView<UnitLower>().solveInPlace(A12);
00376 
00377         A22.noalias() -= A21 * A12;
00378       }
00379     }
00380     return first_zero_pivot;
00381   }
00382 };
00383 
00386 template<typename MatrixType, typename TranspositionType>
00387 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::Index& nb_transpositions)
00388 {
00389   eigen_assert(lu.cols() == row_transpositions.size());
00390   eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
00391 
00392   partial_lu_impl
00393     <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::Index>
00394     ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
00395 }
00396 
00397 } // end namespace internal
00398 
00399 template<typename MatrixType>
00400 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix)
00401 {
00402   m_lu = matrix;
00403 
00404   eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
00405   const Index size = matrix.rows();
00406 
00407   m_rowsTranspositions.resize(size);
00408 
00409   typename TranspositionType::Index nb_transpositions;
00410   internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
00411   m_det_p = (nb_transpositions%2) ? -1 : 1;
00412 
00413   m_p = m_rowsTranspositions;
00414 
00415   m_isInitialized = true;
00416   return *this;
00417 }
00418 
00419 template<typename MatrixType>
00420 typename internal::traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
00421 {
00422   eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00423   return Scalar(m_det_p) * m_lu.diagonal().prod();
00424 }
00425 
00429 template<typename MatrixType>
00430 MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
00431 {
00432   eigen_assert(m_isInitialized && "LU is not initialized.");
00433   // LU
00434   MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
00435                  * m_lu.template triangularView<Upper>();
00436 
00437   // P^{-1}(LU)
00438   res = m_p.inverse() * res;
00439 
00440   return res;
00441 }
00442 
00443 /***** Implementation of solve() *****************************************************/
00444 
00445 namespace internal {
00446 
00447 template<typename _MatrixType, typename Rhs>
00448 struct solve_retval<PartialPivLU<_MatrixType>, Rhs>
00449   : solve_retval_base<PartialPivLU<_MatrixType>, Rhs>
00450 {
00451   EIGEN_MAKE_SOLVE_HELPERS(PartialPivLU<_MatrixType>,Rhs)
00452 
00453   template<typename Dest> void evalTo(Dest& dst) const
00454   {
00455     /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
00456     * So we proceed as follows:
00457     * Step 1: compute c = Pb.
00458     * Step 2: replace c by the solution x to Lx = c.
00459     * Step 3: replace c by the solution x to Ux = c.
00460     */
00461 
00462     eigen_assert(rhs().rows() == dec().matrixLU().rows());
00463 
00464     // Step 1
00465     dst = dec().permutationP() * rhs();
00466 
00467     // Step 2
00468     dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst);
00469 
00470     // Step 3
00471     dec().matrixLU().template triangularView<Upper>().solveInPlace(dst);
00472   }
00473 };
00474 
00475 } // end namespace internal
00476 
00477 /******** MatrixBase methods *******/
00478 
00485 template<typename Derived>
00486 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
00487 MatrixBase<Derived>::partialPivLu() const
00488 {
00489   return PartialPivLU<PlainObject>(eval());
00490 }
00491 
00492 #if EIGEN2_SUPPORT_STAGE > STAGE20_RESOLVE_API_CONFLICTS
00493 
00501 template<typename Derived>
00502 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
00503 MatrixBase<Derived>::lu() const
00504 {
00505   return PartialPivLU<PlainObject>(eval());
00506 }
00507 #endif
00508 
00509 #endif // EIGEN_PARTIALLU_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:10