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