Eigen  3.3.3
IncompleteLUT.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
00005 // Copyright (C) 2014 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_INCOMPLETE_LUT_H
00012 #define EIGEN_INCOMPLETE_LUT_H
00013 
00014 
00015 namespace Eigen { 
00016 
00017 namespace internal {
00018     
00028 template <typename VectorV, typename VectorI>
00029 Index QuickSplit(VectorV &row, VectorI &ind, Index ncut)
00030 {
00031   typedef typename VectorV::RealScalar RealScalar;
00032   using std::swap;
00033   using std::abs;
00034   Index mid;
00035   Index n = row.size(); /* length of the vector */
00036   Index first, last ;
00037   
00038   ncut--; /* to fit the zero-based indices */
00039   first = 0; 
00040   last = n-1; 
00041   if (ncut < first || ncut > last ) return 0;
00042   
00043   do {
00044     mid = first; 
00045     RealScalar abskey = abs(row(mid)); 
00046     for (Index j = first + 1; j <= last; j++) {
00047       if ( abs(row(j)) > abskey) {
00048         ++mid;
00049         swap(row(mid), row(j));
00050         swap(ind(mid), ind(j));
00051       }
00052     }
00053     /* Interchange for the pivot element */
00054     swap(row(mid), row(first));
00055     swap(ind(mid), ind(first));
00056     
00057     if (mid > ncut) last = mid - 1;
00058     else if (mid < ncut ) first = mid + 1; 
00059   } while (mid != ncut );
00060   
00061   return 0; /* mid is equal to ncut */ 
00062 }
00063 
00064 }// end namespace internal
00065 
00098 template <typename _Scalar, typename _StorageIndex = int>
00099 class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar, _StorageIndex> >
00100 {
00101   protected:
00102     typedef SparseSolverBase<IncompleteLUT> Base;
00103     using Base::m_isInitialized;
00104   public:
00105     typedef _Scalar Scalar;
00106     typedef _StorageIndex StorageIndex;
00107     typedef typename NumTraits<Scalar>::Real RealScalar;
00108     typedef Matrix<Scalar,Dynamic,1> Vector;
00109     typedef Matrix<StorageIndex,Dynamic,1> VectorI;
00110     typedef SparseMatrix<Scalar,RowMajor,StorageIndex> FactorType;
00111 
00112     enum {
00113       ColsAtCompileTime = Dynamic,
00114       MaxColsAtCompileTime = Dynamic
00115     };
00116 
00117   public:
00118     
00119     IncompleteLUT()
00120       : m_droptol(NumTraits<Scalar>::dummy_precision()), m_fillfactor(10),
00121         m_analysisIsOk(false), m_factorizationIsOk(false)
00122     {}
00123     
00124     template<typename MatrixType>
00125     explicit IncompleteLUT(const MatrixType& mat, const RealScalar& droptol=NumTraits<Scalar>::dummy_precision(), int fillfactor = 10)
00126       : m_droptol(droptol),m_fillfactor(fillfactor),
00127         m_analysisIsOk(false),m_factorizationIsOk(false)
00128     {
00129       eigen_assert(fillfactor != 0);
00130       compute(mat); 
00131     }
00132     
00133     Index rows() const { return m_lu.rows(); }
00134     
00135     Index cols() const { return m_lu.cols(); }
00136 
00142     ComputationInfo info() const
00143     {
00144       eigen_assert(m_isInitialized && "IncompleteLUT is not initialized.");
00145       return m_info;
00146     }
00147     
00148     template<typename MatrixType>
00149     void analyzePattern(const MatrixType& amat);
00150     
00151     template<typename MatrixType>
00152     void factorize(const MatrixType& amat);
00153     
00159     template<typename MatrixType>
00160     IncompleteLUT& compute(const MatrixType& amat)
00161     {
00162       analyzePattern(amat); 
00163       factorize(amat);
00164       return *this;
00165     }
00166 
00167     void setDroptol(const RealScalar& droptol); 
00168     void setFillfactor(int fillfactor); 
00169     
00170     template<typename Rhs, typename Dest>
00171     void _solve_impl(const Rhs& b, Dest& x) const
00172     {
00173       x = m_Pinv * b;
00174       x = m_lu.template triangularView<UnitLower>().solve(x);
00175       x = m_lu.template triangularView<Upper>().solve(x);
00176       x = m_P * x; 
00177     }
00178 
00179 protected:
00180 
00182     struct keep_diag {
00183       inline bool operator() (const Index& row, const Index& col, const Scalar&) const
00184       {
00185         return row!=col;
00186       }
00187     };
00188 
00189 protected:
00190 
00191     FactorType m_lu;
00192     RealScalar m_droptol;
00193     int m_fillfactor;
00194     bool m_analysisIsOk;
00195     bool m_factorizationIsOk;
00196     ComputationInfo m_info;
00197     PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_P;     // Fill-reducing permutation
00198     PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv;  // Inverse permutation
00199 };
00200 
00205 template<typename Scalar, typename StorageIndex>
00206 void IncompleteLUT<Scalar,StorageIndex>::setDroptol(const RealScalar& droptol)
00207 {
00208   this->m_droptol = droptol;   
00209 }
00210 
00215 template<typename Scalar, typename StorageIndex>
00216 void IncompleteLUT<Scalar,StorageIndex>::setFillfactor(int fillfactor)
00217 {
00218   this->m_fillfactor = fillfactor;   
00219 }
00220 
00221 template <typename Scalar, typename StorageIndex>
00222 template<typename _MatrixType>
00223 void IncompleteLUT<Scalar,StorageIndex>::analyzePattern(const _MatrixType& amat)
00224 {
00225   // Compute the Fill-reducing permutation
00226   // Since ILUT does not perform any numerical pivoting,
00227   // it is highly preferable to keep the diagonal through symmetric permutations.
00228 #ifndef EIGEN_MPL2_ONLY
00229   // To this end, let's symmetrize the pattern and perform AMD on it.
00230   SparseMatrix<Scalar,ColMajor, StorageIndex> mat1 = amat;
00231   SparseMatrix<Scalar,ColMajor, StorageIndex> mat2 = amat.transpose();
00232   // FIXME for a matrix with nearly symmetric pattern, mat2+mat1 is the appropriate choice.
00233   //       on the other hand for a really non-symmetric pattern, mat2*mat1 should be prefered...
00234   SparseMatrix<Scalar,ColMajor, StorageIndex> AtA = mat2 + mat1;
00235   AMDOrdering<StorageIndex> ordering;
00236   ordering(AtA,m_P);
00237   m_Pinv  = m_P.inverse(); // cache the inverse permutation
00238 #else
00239   // If AMD is not available, (MPL2-only), then let's use the slower COLAMD routine.
00240   SparseMatrix<Scalar,ColMajor, StorageIndex> mat1 = amat;
00241   COLAMDOrdering<StorageIndex> ordering;
00242   ordering(mat1,m_Pinv);
00243   m_P = m_Pinv.inverse();
00244 #endif
00245 
00246   m_analysisIsOk = true;
00247   m_factorizationIsOk = false;
00248   m_isInitialized = true;
00249 }
00250 
00251 template <typename Scalar, typename StorageIndex>
00252 template<typename _MatrixType>
00253 void IncompleteLUT<Scalar,StorageIndex>::factorize(const _MatrixType& amat)
00254 {
00255   using std::sqrt;
00256   using std::swap;
00257   using std::abs;
00258   using internal::convert_index;
00259 
00260   eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix");
00261   Index n = amat.cols();  // Size of the matrix
00262   m_lu.resize(n,n);
00263   // Declare Working vectors and variables
00264   Vector u(n) ;     // real values of the row -- maximum size is n --
00265   VectorI ju(n);   // column position of the values in u -- maximum size  is n
00266   VectorI jr(n);   // Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1
00267 
00268   // Apply the fill-reducing permutation
00269   eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00270   SparseMatrix<Scalar,RowMajor, StorageIndex> mat;
00271   mat = amat.twistedBy(m_Pinv);
00272 
00273   // Initialization
00274   jr.fill(-1);
00275   ju.fill(0);
00276   u.fill(0);
00277 
00278   // number of largest elements to keep in each row:
00279   Index fill_in = (amat.nonZeros()*m_fillfactor)/n + 1;
00280   if (fill_in > n) fill_in = n;
00281 
00282   // number of largest nonzero elements to keep in the L and the U part of the current row:
00283   Index nnzL = fill_in/2;
00284   Index nnzU = nnzL;
00285   m_lu.reserve(n * (nnzL + nnzU + 1));
00286 
00287   // global loop over the rows of the sparse matrix
00288   for (Index ii = 0; ii < n; ii++)
00289   {
00290     // 1 - copy the lower and the upper part of the row i of mat in the working vector u
00291 
00292     Index sizeu = 1; // number of nonzero elements in the upper part of the current row
00293     Index sizel = 0; // number of nonzero elements in the lower part of the current row
00294     ju(ii)    = convert_index<StorageIndex>(ii);
00295     u(ii)     = 0;
00296     jr(ii)    = convert_index<StorageIndex>(ii);
00297     RealScalar rownorm = 0;
00298 
00299     typename FactorType::InnerIterator j_it(mat, ii); // Iterate through the current row ii
00300     for (; j_it; ++j_it)
00301     {
00302       Index k = j_it.index();
00303       if (k < ii)
00304       {
00305         // copy the lower part
00306         ju(sizel) = convert_index<StorageIndex>(k);
00307         u(sizel) = j_it.value();
00308         jr(k) = convert_index<StorageIndex>(sizel);
00309         ++sizel;
00310       }
00311       else if (k == ii)
00312       {
00313         u(ii) = j_it.value();
00314       }
00315       else
00316       {
00317         // copy the upper part
00318         Index jpos = ii + sizeu;
00319         ju(jpos) = convert_index<StorageIndex>(k);
00320         u(jpos) = j_it.value();
00321         jr(k) = convert_index<StorageIndex>(jpos);
00322         ++sizeu;
00323       }
00324       rownorm += numext::abs2(j_it.value());
00325     }
00326 
00327     // 2 - detect possible zero row
00328     if(rownorm==0)
00329     {
00330       m_info = NumericalIssue;
00331       return;
00332     }
00333     // Take the 2-norm of the current row as a relative tolerance
00334     rownorm = sqrt(rownorm);
00335 
00336     // 3 - eliminate the previous nonzero rows
00337     Index jj = 0;
00338     Index len = 0;
00339     while (jj < sizel)
00340     {
00341       // In order to eliminate in the correct order,
00342       // we must select first the smallest column index among  ju(jj:sizel)
00343       Index k;
00344       Index minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment
00345       k += jj;
00346       if (minrow != ju(jj))
00347       {
00348         // swap the two locations
00349         Index j = ju(jj);
00350         swap(ju(jj), ju(k));
00351         jr(minrow) = convert_index<StorageIndex>(jj);
00352         jr(j) = convert_index<StorageIndex>(k);
00353         swap(u(jj), u(k));
00354       }
00355       // Reset this location
00356       jr(minrow) = -1;
00357 
00358       // Start elimination
00359       typename FactorType::InnerIterator ki_it(m_lu, minrow);
00360       while (ki_it && ki_it.index() < minrow) ++ki_it;
00361       eigen_internal_assert(ki_it && ki_it.col()==minrow);
00362       Scalar fact = u(jj) / ki_it.value();
00363 
00364       // drop too small elements
00365       if(abs(fact) <= m_droptol)
00366       {
00367         jj++;
00368         continue;
00369       }
00370 
00371       // linear combination of the current row ii and the row minrow
00372       ++ki_it;
00373       for (; ki_it; ++ki_it)
00374       {
00375         Scalar prod = fact * ki_it.value();
00376         Index j     = ki_it.index();
00377         Index jpos  = jr(j);
00378         if (jpos == -1) // fill-in element
00379         {
00380           Index newpos;
00381           if (j >= ii) // dealing with the upper part
00382           {
00383             newpos = ii + sizeu;
00384             sizeu++;
00385             eigen_internal_assert(sizeu<=n);
00386           }
00387           else // dealing with the lower part
00388           {
00389             newpos = sizel;
00390             sizel++;
00391             eigen_internal_assert(sizel<=ii);
00392           }
00393           ju(newpos) = convert_index<StorageIndex>(j);
00394           u(newpos) = -prod;
00395           jr(j) = convert_index<StorageIndex>(newpos);
00396         }
00397         else
00398           u(jpos) -= prod;
00399       }
00400       // store the pivot element
00401       u(len)  = fact;
00402       ju(len) = convert_index<StorageIndex>(minrow);
00403       ++len;
00404 
00405       jj++;
00406     } // end of the elimination on the row ii
00407 
00408     // reset the upper part of the pointer jr to zero
00409     for(Index k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1;
00410 
00411     // 4 - partially sort and insert the elements in the m_lu matrix
00412 
00413     // sort the L-part of the row
00414     sizel = len;
00415     len = (std::min)(sizel, nnzL);
00416     typename Vector::SegmentReturnType ul(u.segment(0, sizel));
00417     typename VectorI::SegmentReturnType jul(ju.segment(0, sizel));
00418     internal::QuickSplit(ul, jul, len);
00419 
00420     // store the largest m_fill elements of the L part
00421     m_lu.startVec(ii);
00422     for(Index k = 0; k < len; k++)
00423       m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
00424 
00425     // store the diagonal element
00426     // apply a shifting rule to avoid zero pivots (we are doing an incomplete factorization)
00427     if (u(ii) == Scalar(0))
00428       u(ii) = sqrt(m_droptol) * rownorm;
00429     m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii);
00430 
00431     // sort the U-part of the row
00432     // apply the dropping rule first
00433     len = 0;
00434     for(Index k = 1; k < sizeu; k++)
00435     {
00436       if(abs(u(ii+k)) > m_droptol * rownorm )
00437       {
00438         ++len;
00439         u(ii + len)  = u(ii + k);
00440         ju(ii + len) = ju(ii + k);
00441       }
00442     }
00443     sizeu = len + 1; // +1 to take into account the diagonal element
00444     len = (std::min)(sizeu, nnzU);
00445     typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
00446     typename VectorI::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
00447     internal::QuickSplit(uu, juu, len);
00448 
00449     // store the largest elements of the U part
00450     for(Index k = ii + 1; k < ii + len; k++)
00451       m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
00452   }
00453   m_lu.finalize();
00454   m_lu.makeCompressed();
00455 
00456   m_factorizationIsOk = true;
00457   m_info = Success;
00458 }
00459 
00460 } // end namespace Eigen
00461 
00462 #endif // EIGEN_INCOMPLETE_LUT_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends