Eigen  3.3.3
PartialPivLU.h
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 
00016 namespace internal {
00017 template<typename _MatrixType> struct traits<PartialPivLU<_MatrixType> >
00018  : traits<_MatrixType>
00019 {
00020   typedef MatrixXpr XprKind;
00021   typedef SolverStorage StorageKind;
00022   typedef traits<_MatrixType> BaseTraits;
00023   enum {
00024     Flags = BaseTraits::Flags & RowMajorBit,
00025     CoeffReadCost = Dynamic
00026   };
00027 };
00028 
00029 template<typename T,typename Derived>
00030 struct enable_if_ref;
00031 // {
00032 //   typedef Derived type;
00033 // };
00034 
00035 template<typename T,typename Derived>
00036 struct enable_if_ref<Ref<T>,Derived> {
00037   typedef Derived type;
00038 };
00039 
00040 } // end namespace internal
00041 
00075 template<typename _MatrixType> class PartialPivLU
00076   : public SolverBase<PartialPivLU<_MatrixType> >
00077 {
00078   public:
00079 
00080     typedef _MatrixType MatrixType;
00081     typedef SolverBase<PartialPivLU> Base;
00082     EIGEN_GENERIC_PUBLIC_INTERFACE(PartialPivLU)
00083     // FIXME StorageIndex defined in EIGEN_GENERIC_PUBLIC_INTERFACE should be int
00084     enum {
00085       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00086       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00087     };
00088     typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
00089     typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
00090     typedef typename MatrixType::PlainObject PlainObject;
00091 
00098     PartialPivLU();
00099 
00106     explicit PartialPivLU(Index size);
00107 
00115     template<typename InputType>
00116     explicit PartialPivLU(const EigenBase<InputType>& matrix);
00117 
00125     template<typename InputType>
00126     explicit PartialPivLU(EigenBase<InputType>& matrix);
00127 
00128     template<typename InputType>
00129     PartialPivLU& compute(const EigenBase<InputType>& matrix) {
00130       m_lu = matrix.derived();
00131       compute();
00132       return *this;
00133     }
00134 
00141     inline const MatrixType& matrixLU() const
00142     {
00143       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00144       return m_lu;
00145     }
00146 
00149     inline const PermutationType& permutationP() const
00150     {
00151       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00152       return m_p;
00153     }
00154 
00172     // FIXME this is a copy-paste of the base-class member to add the isInitialized assertion.
00173     template<typename Rhs>
00174     inline const Solve<PartialPivLU, Rhs>
00175     solve(const MatrixBase<Rhs>& b) const
00176     {
00177       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00178       return Solve<PartialPivLU, Rhs>(*this, b.derived());
00179     }
00180 
00184     inline RealScalar rcond() const
00185     {
00186       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00187       return internal::rcond_estimate_helper(m_l1_norm, *this);
00188     }
00189 
00197     inline const Inverse<PartialPivLU> inverse() const
00198     {
00199       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00200       return Inverse<PartialPivLU>(*this);
00201     }
00202 
00216     Scalar determinant() const;
00217 
00218     MatrixType reconstructedMatrix() const;
00219 
00220     inline Index rows() const { return m_lu.rows(); }
00221     inline Index cols() const { return m_lu.cols(); }
00222 
00223     #ifndef EIGEN_PARSED_BY_DOXYGEN
00224     template<typename RhsType, typename DstType>
00225     EIGEN_DEVICE_FUNC
00226     void _solve_impl(const RhsType &rhs, DstType &dst) const {
00227      /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
00228       * So we proceed as follows:
00229       * Step 1: compute c = Pb.
00230       * Step 2: replace c by the solution x to Lx = c.
00231       * Step 3: replace c by the solution x to Ux = c.
00232       */
00233 
00234       eigen_assert(rhs.rows() == m_lu.rows());
00235 
00236       // Step 1
00237       dst = permutationP() * rhs;
00238 
00239       // Step 2
00240       m_lu.template triangularView<UnitLower>().solveInPlace(dst);
00241 
00242       // Step 3
00243       m_lu.template triangularView<Upper>().solveInPlace(dst);
00244     }
00245 
00246     template<bool Conjugate, typename RhsType, typename DstType>
00247     EIGEN_DEVICE_FUNC
00248     void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const {
00249      /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
00250       * So we proceed as follows:
00251       * Step 1: compute c = Pb.
00252       * Step 2: replace c by the solution x to Lx = c.
00253       * Step 3: replace c by the solution x to Ux = c.
00254       */
00255 
00256       eigen_assert(rhs.rows() == m_lu.cols());
00257 
00258       if (Conjugate) {
00259         // Step 1
00260         dst = m_lu.template triangularView<Upper>().adjoint().solve(rhs);
00261         // Step 2
00262         m_lu.template triangularView<UnitLower>().adjoint().solveInPlace(dst);
00263       } else {
00264         // Step 1
00265         dst = m_lu.template triangularView<Upper>().transpose().solve(rhs);
00266         // Step 2
00267         m_lu.template triangularView<UnitLower>().transpose().solveInPlace(dst);
00268       }
00269       // Step 3
00270       dst = permutationP().transpose() * dst;
00271     }
00272     #endif
00273 
00274   protected:
00275 
00276     static void check_template_parameters()
00277     {
00278       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00279     }
00280 
00281     void compute();
00282 
00283     MatrixType m_lu;
00284     PermutationType m_p;
00285     TranspositionType m_rowsTranspositions;
00286     RealScalar m_l1_norm;
00287     signed char m_det_p;
00288     bool m_isInitialized;
00289 };
00290 
00291 template<typename MatrixType>
00292 PartialPivLU<MatrixType>::PartialPivLU()
00293   : m_lu(),
00294     m_p(),
00295     m_rowsTranspositions(),
00296     m_l1_norm(0),
00297     m_det_p(0),
00298     m_isInitialized(false)
00299 {
00300 }
00301 
00302 template<typename MatrixType>
00303 PartialPivLU<MatrixType>::PartialPivLU(Index size)
00304   : m_lu(size, size),
00305     m_p(size),
00306     m_rowsTranspositions(size),
00307     m_l1_norm(0),
00308     m_det_p(0),
00309     m_isInitialized(false)
00310 {
00311 }
00312 
00313 template<typename MatrixType>
00314 template<typename InputType>
00315 PartialPivLU<MatrixType>::PartialPivLU(const EigenBase<InputType>& matrix)
00316   : m_lu(matrix.rows(),matrix.cols()),
00317     m_p(matrix.rows()),
00318     m_rowsTranspositions(matrix.rows()),
00319     m_l1_norm(0),
00320     m_det_p(0),
00321     m_isInitialized(false)
00322 {
00323   compute(matrix.derived());
00324 }
00325 
00326 template<typename MatrixType>
00327 template<typename InputType>
00328 PartialPivLU<MatrixType>::PartialPivLU(EigenBase<InputType>& matrix)
00329   : m_lu(matrix.derived()),
00330     m_p(matrix.rows()),
00331     m_rowsTranspositions(matrix.rows()),
00332     m_l1_norm(0),
00333     m_det_p(0),
00334     m_isInitialized(false)
00335 {
00336   compute();
00337 }
00338 
00339 namespace internal {
00340 
00342 template<typename Scalar, int StorageOrder, typename PivIndex>
00343 struct partial_lu_impl
00344 {
00345   // FIXME add a stride to Map, so that the following mapping becomes easier,
00346   // another option would be to create an expression being able to automatically
00347   // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly
00348   // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix,
00349   // and Block.
00350   typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > MapLU;
00351   typedef Block<MapLU, Dynamic, Dynamic> MatrixType;
00352   typedef Block<MatrixType,Dynamic,Dynamic> BlockType;
00353   typedef typename MatrixType::RealScalar RealScalar;
00354 
00365   static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
00366   {
00367     typedef scalar_score_coeff_op<Scalar> Scoring;
00368     typedef typename Scoring::result_type Score;
00369     const Index rows = lu.rows();
00370     const Index cols = lu.cols();
00371     const Index size = (std::min)(rows,cols);
00372     nb_transpositions = 0;
00373     Index first_zero_pivot = -1;
00374     for(Index k = 0; k < size; ++k)
00375     {
00376       Index rrows = rows-k-1;
00377       Index rcols = cols-k-1;
00378 
00379       Index row_of_biggest_in_col;
00380       Score biggest_in_corner
00381         = lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col);
00382       row_of_biggest_in_col += k;
00383 
00384       row_transpositions[k] = PivIndex(row_of_biggest_in_col);
00385 
00386       if(biggest_in_corner != Score(0))
00387       {
00388         if(k != row_of_biggest_in_col)
00389         {
00390           lu.row(k).swap(lu.row(row_of_biggest_in_col));
00391           ++nb_transpositions;
00392         }
00393 
00394         // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k)
00395         // overflow but not the actual quotient?
00396         lu.col(k).tail(rrows) /= lu.coeff(k,k);
00397       }
00398       else if(first_zero_pivot==-1)
00399       {
00400         // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
00401         // and continue the factorization such we still have A = PLU
00402         first_zero_pivot = k;
00403       }
00404 
00405       if(k<rows-1)
00406         lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
00407     }
00408     return first_zero_pivot;
00409   }
00410 
00426   static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
00427   {
00428     MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
00429     MatrixType lu(lu1,0,0,rows,cols);
00430 
00431     const Index size = (std::min)(rows,cols);
00432 
00433     // if the matrix is too small, no blocking:
00434     if(size<=16)
00435     {
00436       return unblocked_lu(lu, row_transpositions, nb_transpositions);
00437     }
00438 
00439     // automatically adjust the number of subdivisions to the size
00440     // of the matrix so that there is enough sub blocks:
00441     Index blockSize;
00442     {
00443       blockSize = size/8;
00444       blockSize = (blockSize/16)*16;
00445       blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
00446     }
00447 
00448     nb_transpositions = 0;
00449     Index first_zero_pivot = -1;
00450     for(Index k = 0; k < size; k+=blockSize)
00451     {
00452       Index bs = (std::min)(size-k,blockSize); // actual size of the block
00453       Index trows = rows - k - bs; // trailing rows
00454       Index tsize = size - k - bs; // trailing size
00455 
00456       // partition the matrix:
00457       //                          A00 | A01 | A02
00458       // lu  = A_0 | A_1 | A_2 =  A10 | A11 | A12
00459       //                          A20 | A21 | A22
00460       BlockType A_0(lu,0,0,rows,k);
00461       BlockType A_2(lu,0,k+bs,rows,tsize);
00462       BlockType A11(lu,k,k,bs,bs);
00463       BlockType A12(lu,k,k+bs,bs,tsize);
00464       BlockType A21(lu,k+bs,k,trows,bs);
00465       BlockType A22(lu,k+bs,k+bs,trows,tsize);
00466 
00467       PivIndex nb_transpositions_in_panel;
00468       // recursively call the blocked LU algorithm on [A11^T A21^T]^T
00469       // with a very small blocking size:
00470       Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
00471                    row_transpositions+k, nb_transpositions_in_panel, 16);
00472       if(ret>=0 && first_zero_pivot==-1)
00473         first_zero_pivot = k+ret;
00474 
00475       nb_transpositions += nb_transpositions_in_panel;
00476       // update permutations and apply them to A_0
00477       for(Index i=k; i<k+bs; ++i)
00478       {
00479         Index piv = (row_transpositions[i] += internal::convert_index<PivIndex>(k));
00480         A_0.row(i).swap(A_0.row(piv));
00481       }
00482 
00483       if(trows)
00484       {
00485         // apply permutations to A_2
00486         for(Index i=k;i<k+bs; ++i)
00487           A_2.row(i).swap(A_2.row(row_transpositions[i]));
00488 
00489         // A12 = A11^-1 A12
00490         A11.template triangularView<UnitLower>().solveInPlace(A12);
00491 
00492         A22.noalias() -= A21 * A12;
00493       }
00494     }
00495     return first_zero_pivot;
00496   }
00497 };
00498 
00501 template<typename MatrixType, typename TranspositionType>
00502 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::StorageIndex& nb_transpositions)
00503 {
00504   eigen_assert(lu.cols() == row_transpositions.size());
00505   eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
00506 
00507   partial_lu_impl
00508     <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::StorageIndex>
00509     ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
00510 }
00511 
00512 } // end namespace internal
00513 
00514 template<typename MatrixType>
00515 void PartialPivLU<MatrixType>::compute()
00516 {
00517   check_template_parameters();
00518 
00519   // the row permutation is stored as int indices, so just to be sure:
00520   eigen_assert(m_lu.rows()<NumTraits<int>::highest());
00521 
00522   m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
00523 
00524   eigen_assert(m_lu.rows() == m_lu.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
00525   const Index size = m_lu.rows();
00526 
00527   m_rowsTranspositions.resize(size);
00528 
00529   typename TranspositionType::StorageIndex nb_transpositions;
00530   internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
00531   m_det_p = (nb_transpositions%2) ? -1 : 1;
00532 
00533   m_p = m_rowsTranspositions;
00534 
00535   m_isInitialized = true;
00536 }
00537 
00538 template<typename MatrixType>
00539 typename PartialPivLU<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
00540 {
00541   eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00542   return Scalar(m_det_p) * m_lu.diagonal().prod();
00543 }
00544 
00548 template<typename MatrixType>
00549 MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
00550 {
00551   eigen_assert(m_isInitialized && "LU is not initialized.");
00552   // LU
00553   MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
00554                  * m_lu.template triangularView<Upper>();
00555 
00556   // P^{-1}(LU)
00557   res = m_p.inverse() * res;
00558 
00559   return res;
00560 }
00561 
00562 /***** Implementation details *****************************************************/
00563 
00564 namespace internal {
00565 
00566 /***** Implementation of inverse() *****************************************************/
00567 template<typename DstXprType, typename MatrixType>
00568 struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType>::Scalar>, Dense2Dense>
00569 {
00570   typedef PartialPivLU<MatrixType> LuType;
00571   typedef Inverse<LuType> SrcXprType;
00572   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename LuType::Scalar> &)
00573   {
00574     dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
00575   }
00576 };
00577 } // end namespace internal
00578 
00579 /******** MatrixBase methods *******/
00580 
00587 template<typename Derived>
00588 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
00589 MatrixBase<Derived>::partialPivLu() const
00590 {
00591   return PartialPivLU<PlainObject>(eval());
00592 }
00593 
00602 template<typename Derived>
00603 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
00604 MatrixBase<Derived>::lu() const
00605 {
00606   return PartialPivLU<PlainObject>(eval());
00607 }
00608 
00609 } // end namespace Eigen
00610 
00611 #endif // EIGEN_PARTIALLU_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends