![]() |
Eigen
3.3.3
|
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