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