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