1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_INCOMPLETE_LUT_H 11 #define EIGEN_INCOMPLETE_LUT_H 12 13 14 namespace Eigen { 15 16 namespace internal { 17 18 /** \internal 19 * Compute a quick-sort split of a vector 20 * On output, the vector row is permuted such that its elements satisfy 21 * abs(row(i)) >= abs(row(ncut)) if i<ncut 22 * abs(row(i)) <= abs(row(ncut)) if i>ncut 23 * \param row The vector of values 24 * \param ind The array of index for the elements in @p row 25 * \param ncut The number of largest elements to keep 26 **/ 27 template <typename VectorV, typename VectorI, typename Index> 28 Index QuickSplit(VectorV &row, VectorI &ind, Index ncut) 29 { 30 typedef typename VectorV::RealScalar RealScalar; 31 using std::swap; 32 using std::abs; 33 Index mid; 34 Index n = row.size(); /* length of the vector */ 35 Index first, last ; 36 37 ncut--; /* to fit the zero-based indices */ 38 first = 0; 39 last = n-1; 40 if (ncut < first || ncut > last ) return 0; 41 42 do { 43 mid = first; 44 RealScalar abskey = abs(row(mid)); 45 for (Index j = first + 1; j <= last; j++) { 46 if ( abs(row(j)) > abskey) { 47 ++mid; 48 swap(row(mid), row(j)); 49 swap(ind(mid), ind(j)); 50 } 51 } 52 /* Interchange for the pivot element */ 53 swap(row(mid), row(first)); 54 swap(ind(mid), ind(first)); 55 56 if (mid > ncut) last = mid - 1; 57 else if (mid < ncut ) first = mid + 1; 58 } while (mid != ncut ); 59 60 return 0; /* mid is equal to ncut */ 61 } 62 63 }// end namespace internal 64 65 /** \ingroup IterativeLinearSolvers_Module 66 * \class IncompleteLUT 67 * \brief Incomplete LU factorization with dual-threshold strategy 68 * 69 * During the numerical factorization, two dropping rules are used : 70 * 1) any element whose magnitude is less than some tolerance is dropped. 71 * This tolerance is obtained by multiplying the input tolerance @p droptol 72 * by the average magnitude of all the original elements in the current row. 73 * 2) After the elimination of the row, only the @p fill largest elements in 74 * the L part and the @p fill largest elements in the U part are kept 75 * (in addition to the diagonal element ). Note that @p fill is computed from 76 * the input parameter @p fillfactor which is used the ratio to control the fill_in 77 * relatively to the initial number of nonzero elements. 78 * 79 * The two extreme cases are when @p droptol=0 (to keep all the @p fill*2 largest elements) 80 * and when @p fill=n/2 with @p droptol being different to zero. 81 * 82 * References : Yousef Saad, ILUT: A dual threshold incomplete LU factorization, 83 * Numerical Linear Algebra with Applications, 1(4), pp 387-402, 1994. 84 * 85 * NOTE : The following implementation is derived from the ILUT implementation 86 * in the SPARSKIT package, Copyright (C) 2005, the Regents of the University of Minnesota 87 * released under the terms of the GNU LGPL: 88 * http://www-users.cs.umn.edu/~saad/software/SPARSKIT/README 89 * However, Yousef Saad gave us permission to relicense his ILUT code to MPL2. 90 * See the Eigen mailing list archive, thread: ILUT, date: July 8, 2012: 91 * http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2012/07/msg00064.html 92 * alternatively, on GMANE: 93 * http://comments.gmane.org/gmane.comp.lib.eigen/3302 94 */ 95 template <typename _Scalar> 96 class IncompleteLUT : internal::noncopyable 97 { 98 typedef _Scalar Scalar; 99 typedef typename NumTraits<Scalar>::Real RealScalar; 100 typedef Matrix<Scalar,Dynamic,1> Vector; 101 typedef SparseMatrix<Scalar,RowMajor> FactorType; 102 typedef SparseMatrix<Scalar,ColMajor> PermutType; 103 typedef typename FactorType::Index Index; 104 105 public: 106 typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; 107 108 IncompleteLUT() 109 : m_droptol(NumTraits<Scalar>::dummy_precision()), m_fillfactor(10), 110 m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false) 111 {} 112 113 template<typename MatrixType> 114 IncompleteLUT(const MatrixType& mat, const RealScalar& droptol=NumTraits<Scalar>::dummy_precision(), int fillfactor = 10) 115 : m_droptol(droptol),m_fillfactor(fillfactor), 116 m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false) 117 { 118 eigen_assert(fillfactor != 0); 119 compute(mat); 120 } 121 122 Index rows() const { return m_lu.rows(); } 123 124 Index cols() const { return m_lu.cols(); } 125 126 /** \brief Reports whether previous computation was successful. 127 * 128 * \returns \c Success if computation was succesful, 129 * \c NumericalIssue if the matrix.appears to be negative. 130 */ 131 ComputationInfo info() const 132 { 133 eigen_assert(m_isInitialized && "IncompleteLUT is not initialized."); 134 return m_info; 135 } 136 137 template<typename MatrixType> 138 void analyzePattern(const MatrixType& amat); 139 140 template<typename MatrixType> 141 void factorize(const MatrixType& amat); 142 143 /** 144 * Compute an incomplete LU factorization with dual threshold on the matrix mat 145 * No pivoting is done in this version 146 * 147 **/ 148 template<typename MatrixType> 149 IncompleteLUT<Scalar>& compute(const MatrixType& amat) 150 { 151 analyzePattern(amat); 152 factorize(amat); 153 return *this; 154 } 155 156 void setDroptol(const RealScalar& droptol); 157 void setFillfactor(int fillfactor); 158 159 template<typename Rhs, typename Dest> 160 void _solve(const Rhs& b, Dest& x) const 161 { 162 x = m_Pinv * b; 163 x = m_lu.template triangularView<UnitLower>().solve(x); 164 x = m_lu.template triangularView<Upper>().solve(x); 165 x = m_P * x; 166 } 167 168 template<typename Rhs> inline const internal::solve_retval<IncompleteLUT, Rhs> 169 solve(const MatrixBase<Rhs>& b) const 170 { 171 eigen_assert(m_isInitialized && "IncompleteLUT is not initialized."); 172 eigen_assert(cols()==b.rows() 173 && "IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b"); 174 return internal::solve_retval<IncompleteLUT, Rhs>(*this, b.derived()); 175 } 176 177 protected: 178 179 /** keeps off-diagonal entries; drops diagonal entries */ 180 struct keep_diag { 181 inline bool operator() (const Index& row, const Index& col, const Scalar&) const 182 { 183 return row!=col; 184 } 185 }; 186 187 protected: 188 189 FactorType m_lu; 190 RealScalar m_droptol; 191 int m_fillfactor; 192 bool m_analysisIsOk; 193 bool m_factorizationIsOk; 194 bool m_isInitialized; 195 ComputationInfo m_info; 196 PermutationMatrix<Dynamic,Dynamic,Index> m_P; // Fill-reducing permutation 197 PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // Inverse permutation 198 }; 199 200 /** 201 * Set control parameter droptol 202 * \param droptol Drop any element whose magnitude is less than this tolerance 203 **/ 204 template<typename Scalar> 205 void IncompleteLUT<Scalar>::setDroptol(const RealScalar& droptol) 206 { 207 this->m_droptol = droptol; 208 } 209 210 /** 211 * Set control parameter fillfactor 212 * \param fillfactor This is used to compute the number @p fill_in of largest elements to keep on each row. 213 **/ 214 template<typename Scalar> 215 void IncompleteLUT<Scalar>::setFillfactor(int fillfactor) 216 { 217 this->m_fillfactor = fillfactor; 218 } 219 220 template <typename Scalar> 221 template<typename _MatrixType> 222 void IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat) 223 { 224 // Compute the Fill-reducing permutation 225 SparseMatrix<Scalar,ColMajor, Index> mat1 = amat; 226 SparseMatrix<Scalar,ColMajor, Index> mat2 = amat.transpose(); 227 // Symmetrize the pattern 228 // FIXME for a matrix with nearly symmetric pattern, mat2+mat1 is the appropriate choice. 229 // on the other hand for a really non-symmetric pattern, mat2*mat1 should be prefered... 230 SparseMatrix<Scalar,ColMajor, Index> AtA = mat2 + mat1; 231 AtA.prune(keep_diag()); 232 internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P); // Then compute the AMD ordering... 233 234 m_Pinv = m_P.inverse(); // ... and the inverse permutation 235 236 m_analysisIsOk = true; 237 m_factorizationIsOk = false; 238 m_isInitialized = false; 239 } 240 241 template <typename Scalar> 242 template<typename _MatrixType> 243 void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) 244 { 245 using std::sqrt; 246 using std::swap; 247 using std::abs; 248 249 eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix"); 250 Index n = amat.cols(); // Size of the matrix 251 m_lu.resize(n,n); 252 // Declare Working vectors and variables 253 Vector u(n) ; // real values of the row -- maximum size is n -- 254 VectorXi ju(n); // column position of the values in u -- maximum size is n 255 VectorXi jr(n); // Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1 256 257 // Apply the fill-reducing permutation 258 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 259 SparseMatrix<Scalar,RowMajor, Index> mat; 260 mat = amat.twistedBy(m_Pinv); 261 262 // Initialization 263 jr.fill(-1); 264 ju.fill(0); 265 u.fill(0); 266 267 // number of largest elements to keep in each row: 268 Index fill_in = static_cast<Index> (amat.nonZeros()*m_fillfactor)/n+1; 269 if (fill_in > n) fill_in = n; 270 271 // number of largest nonzero elements to keep in the L and the U part of the current row: 272 Index nnzL = fill_in/2; 273 Index nnzU = nnzL; 274 m_lu.reserve(n * (nnzL + nnzU + 1)); 275 276 // global loop over the rows of the sparse matrix 277 for (Index ii = 0; ii < n; ii++) 278 { 279 // 1 - copy the lower and the upper part of the row i of mat in the working vector u 280 281 Index sizeu = 1; // number of nonzero elements in the upper part of the current row 282 Index sizel = 0; // number of nonzero elements in the lower part of the current row 283 ju(ii) = ii; 284 u(ii) = 0; 285 jr(ii) = ii; 286 RealScalar rownorm = 0; 287 288 typename FactorType::InnerIterator j_it(mat, ii); // Iterate through the current row ii 289 for (; j_it; ++j_it) 290 { 291 Index k = j_it.index(); 292 if (k < ii) 293 { 294 // copy the lower part 295 ju(sizel) = k; 296 u(sizel) = j_it.value(); 297 jr(k) = sizel; 298 ++sizel; 299 } 300 else if (k == ii) 301 { 302 u(ii) = j_it.value(); 303 } 304 else 305 { 306 // copy the upper part 307 Index jpos = ii + sizeu; 308 ju(jpos) = k; 309 u(jpos) = j_it.value(); 310 jr(k) = jpos; 311 ++sizeu; 312 } 313 rownorm += numext::abs2(j_it.value()); 314 } 315 316 // 2 - detect possible zero row 317 if(rownorm==0) 318 { 319 m_info = NumericalIssue; 320 return; 321 } 322 // Take the 2-norm of the current row as a relative tolerance 323 rownorm = sqrt(rownorm); 324 325 // 3 - eliminate the previous nonzero rows 326 Index jj = 0; 327 Index len = 0; 328 while (jj < sizel) 329 { 330 // In order to eliminate in the correct order, 331 // we must select first the smallest column index among ju(jj:sizel) 332 Index k; 333 Index minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment 334 k += jj; 335 if (minrow != ju(jj)) 336 { 337 // swap the two locations 338 Index j = ju(jj); 339 swap(ju(jj), ju(k)); 340 jr(minrow) = jj; jr(j) = k; 341 swap(u(jj), u(k)); 342 } 343 // Reset this location 344 jr(minrow) = -1; 345 346 // Start elimination 347 typename FactorType::InnerIterator ki_it(m_lu, minrow); 348 while (ki_it && ki_it.index() < minrow) ++ki_it; 349 eigen_internal_assert(ki_it && ki_it.col()==minrow); 350 Scalar fact = u(jj) / ki_it.value(); 351 352 // drop too small elements 353 if(abs(fact) <= m_droptol) 354 { 355 jj++; 356 continue; 357 } 358 359 // linear combination of the current row ii and the row minrow 360 ++ki_it; 361 for (; ki_it; ++ki_it) 362 { 363 Scalar prod = fact * ki_it.value(); 364 Index j = ki_it.index(); 365 Index jpos = jr(j); 366 if (jpos == -1) // fill-in element 367 { 368 Index newpos; 369 if (j >= ii) // dealing with the upper part 370 { 371 newpos = ii + sizeu; 372 sizeu++; 373 eigen_internal_assert(sizeu<=n); 374 } 375 else // dealing with the lower part 376 { 377 newpos = sizel; 378 sizel++; 379 eigen_internal_assert(sizel<=ii); 380 } 381 ju(newpos) = j; 382 u(newpos) = -prod; 383 jr(j) = newpos; 384 } 385 else 386 u(jpos) -= prod; 387 } 388 // store the pivot element 389 u(len) = fact; 390 ju(len) = minrow; 391 ++len; 392 393 jj++; 394 } // end of the elimination on the row ii 395 396 // reset the upper part of the pointer jr to zero 397 for(Index k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1; 398 399 // 4 - partially sort and insert the elements in the m_lu matrix 400 401 // sort the L-part of the row 402 sizel = len; 403 len = (std::min)(sizel, nnzL); 404 typename Vector::SegmentReturnType ul(u.segment(0, sizel)); 405 typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel)); 406 internal::QuickSplit(ul, jul, len); 407 408 // store the largest m_fill elements of the L part 409 m_lu.startVec(ii); 410 for(Index k = 0; k < len; k++) 411 m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k); 412 413 // store the diagonal element 414 // apply a shifting rule to avoid zero pivots (we are doing an incomplete factorization) 415 if (u(ii) == Scalar(0)) 416 u(ii) = sqrt(m_droptol) * rownorm; 417 m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii); 418 419 // sort the U-part of the row 420 // apply the dropping rule first 421 len = 0; 422 for(Index k = 1; k < sizeu; k++) 423 { 424 if(abs(u(ii+k)) > m_droptol * rownorm ) 425 { 426 ++len; 427 u(ii + len) = u(ii + k); 428 ju(ii + len) = ju(ii + k); 429 } 430 } 431 sizeu = len + 1; // +1 to take into account the diagonal element 432 len = (std::min)(sizeu, nnzU); 433 typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1)); 434 typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1)); 435 internal::QuickSplit(uu, juu, len); 436 437 // store the largest elements of the U part 438 for(Index k = ii + 1; k < ii + len; k++) 439 m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k); 440 } 441 442 m_lu.finalize(); 443 m_lu.makeCompressed(); 444 445 m_factorizationIsOk = true; 446 m_isInitialized = m_factorizationIsOk; 447 m_info = Success; 448 } 449 450 namespace internal { 451 452 template<typename _MatrixType, typename Rhs> 453 struct solve_retval<IncompleteLUT<_MatrixType>, Rhs> 454 : solve_retval_base<IncompleteLUT<_MatrixType>, Rhs> 455 { 456 typedef IncompleteLUT<_MatrixType> Dec; 457 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 458 459 template<typename Dest> void evalTo(Dest& dst) const 460 { 461 dec()._solve(rhs(),dst); 462 } 463 }; 464 465 } // end namespace internal 466 467 } // end namespace Eigen 468 469 #endif // EIGEN_INCOMPLETE_LUT_H 470