1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012 Désiré Nuentsa-Wakam <[email protected]>
5 // Copyright (C) 2014 Gael Guennebaud <[email protected]>
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>
QuickSplit(VectorV & row,VectorI & ind,Index ncut)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
IncompleteLUT()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)
m_droptol(droptol)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
rows()133 EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_lu.rows(); }
134
cols()135 EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_lu.cols(); }
136
137 /** \brief Reports whether previous computation was successful.
138 *
139 * \returns \c Success if computation was successful,
140 * \c NumericalIssue if the matrix.appears to be negative.
141 */
info()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>
compute(const MatrixType & amat)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>
_solve_impl(const Rhs & b,Dest & x)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 {
operatorkeep_diag183 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>
setDroptol(const RealScalar & droptol)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>
setFillfactor(int fillfactor)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>
analyzePattern(const _MatrixType & amat)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 // To this end, let's symmetrize the pattern and perform AMD on it.
229 SparseMatrix<Scalar,ColMajor, StorageIndex> mat1 = amat;
230 SparseMatrix<Scalar,ColMajor, StorageIndex> mat2 = amat.transpose();
231 // FIXME for a matrix with nearly symmetric pattern, mat2+mat1 is the appropriate choice.
232 // on the other hand for a really non-symmetric pattern, mat2*mat1 should be preferred...
233 SparseMatrix<Scalar,ColMajor, StorageIndex> AtA = mat2 + mat1;
234 AMDOrdering<StorageIndex> ordering;
235 ordering(AtA,m_P);
236 m_Pinv = m_P.inverse(); // cache the inverse permutation
237 m_analysisIsOk = true;
238 m_factorizationIsOk = false;
239 m_isInitialized = true;
240 }
241
242 template <typename Scalar, typename StorageIndex>
243 template<typename _MatrixType>
factorize(const _MatrixType & amat)244 void IncompleteLUT<Scalar,StorageIndex>::factorize(const _MatrixType& amat)
245 {
246 using std::sqrt;
247 using std::swap;
248 using std::abs;
249 using internal::convert_index;
250
251 eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix");
252 Index n = amat.cols(); // Size of the matrix
253 m_lu.resize(n,n);
254 // Declare Working vectors and variables
255 Vector u(n) ; // real values of the row -- maximum size is n --
256 VectorI ju(n); // column position of the values in u -- maximum size is n
257 VectorI jr(n); // Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1
258
259 // Apply the fill-reducing permutation
260 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
261 SparseMatrix<Scalar,RowMajor, StorageIndex> mat;
262 mat = amat.twistedBy(m_Pinv);
263
264 // Initialization
265 jr.fill(-1);
266 ju.fill(0);
267 u.fill(0);
268
269 // number of largest elements to keep in each row:
270 Index fill_in = (amat.nonZeros()*m_fillfactor)/n + 1;
271 if (fill_in > n) fill_in = n;
272
273 // number of largest nonzero elements to keep in the L and the U part of the current row:
274 Index nnzL = fill_in/2;
275 Index nnzU = nnzL;
276 m_lu.reserve(n * (nnzL + nnzU + 1));
277
278 // global loop over the rows of the sparse matrix
279 for (Index ii = 0; ii < n; ii++)
280 {
281 // 1 - copy the lower and the upper part of the row i of mat in the working vector u
282
283 Index sizeu = 1; // number of nonzero elements in the upper part of the current row
284 Index sizel = 0; // number of nonzero elements in the lower part of the current row
285 ju(ii) = convert_index<StorageIndex>(ii);
286 u(ii) = 0;
287 jr(ii) = convert_index<StorageIndex>(ii);
288 RealScalar rownorm = 0;
289
290 typename FactorType::InnerIterator j_it(mat, ii); // Iterate through the current row ii
291 for (; j_it; ++j_it)
292 {
293 Index k = j_it.index();
294 if (k < ii)
295 {
296 // copy the lower part
297 ju(sizel) = convert_index<StorageIndex>(k);
298 u(sizel) = j_it.value();
299 jr(k) = convert_index<StorageIndex>(sizel);
300 ++sizel;
301 }
302 else if (k == ii)
303 {
304 u(ii) = j_it.value();
305 }
306 else
307 {
308 // copy the upper part
309 Index jpos = ii + sizeu;
310 ju(jpos) = convert_index<StorageIndex>(k);
311 u(jpos) = j_it.value();
312 jr(k) = convert_index<StorageIndex>(jpos);
313 ++sizeu;
314 }
315 rownorm += numext::abs2(j_it.value());
316 }
317
318 // 2 - detect possible zero row
319 if(rownorm==0)
320 {
321 m_info = NumericalIssue;
322 return;
323 }
324 // Take the 2-norm of the current row as a relative tolerance
325 rownorm = sqrt(rownorm);
326
327 // 3 - eliminate the previous nonzero rows
328 Index jj = 0;
329 Index len = 0;
330 while (jj < sizel)
331 {
332 // In order to eliminate in the correct order,
333 // we must select first the smallest column index among ju(jj:sizel)
334 Index k;
335 Index minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment
336 k += jj;
337 if (minrow != ju(jj))
338 {
339 // swap the two locations
340 Index j = ju(jj);
341 swap(ju(jj), ju(k));
342 jr(minrow) = convert_index<StorageIndex>(jj);
343 jr(j) = convert_index<StorageIndex>(k);
344 swap(u(jj), u(k));
345 }
346 // Reset this location
347 jr(minrow) = -1;
348
349 // Start elimination
350 typename FactorType::InnerIterator ki_it(m_lu, minrow);
351 while (ki_it && ki_it.index() < minrow) ++ki_it;
352 eigen_internal_assert(ki_it && ki_it.col()==minrow);
353 Scalar fact = u(jj) / ki_it.value();
354
355 // drop too small elements
356 if(abs(fact) <= m_droptol)
357 {
358 jj++;
359 continue;
360 }
361
362 // linear combination of the current row ii and the row minrow
363 ++ki_it;
364 for (; ki_it; ++ki_it)
365 {
366 Scalar prod = fact * ki_it.value();
367 Index j = ki_it.index();
368 Index jpos = jr(j);
369 if (jpos == -1) // fill-in element
370 {
371 Index newpos;
372 if (j >= ii) // dealing with the upper part
373 {
374 newpos = ii + sizeu;
375 sizeu++;
376 eigen_internal_assert(sizeu<=n);
377 }
378 else // dealing with the lower part
379 {
380 newpos = sizel;
381 sizel++;
382 eigen_internal_assert(sizel<=ii);
383 }
384 ju(newpos) = convert_index<StorageIndex>(j);
385 u(newpos) = -prod;
386 jr(j) = convert_index<StorageIndex>(newpos);
387 }
388 else
389 u(jpos) -= prod;
390 }
391 // store the pivot element
392 u(len) = fact;
393 ju(len) = convert_index<StorageIndex>(minrow);
394 ++len;
395
396 jj++;
397 } // end of the elimination on the row ii
398
399 // reset the upper part of the pointer jr to zero
400 for(Index k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1;
401
402 // 4 - partially sort and insert the elements in the m_lu matrix
403
404 // sort the L-part of the row
405 sizel = len;
406 len = (std::min)(sizel, nnzL);
407 typename Vector::SegmentReturnType ul(u.segment(0, sizel));
408 typename VectorI::SegmentReturnType jul(ju.segment(0, sizel));
409 internal::QuickSplit(ul, jul, len);
410
411 // store the largest m_fill elements of the L part
412 m_lu.startVec(ii);
413 for(Index k = 0; k < len; k++)
414 m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
415
416 // store the diagonal element
417 // apply a shifting rule to avoid zero pivots (we are doing an incomplete factorization)
418 if (u(ii) == Scalar(0))
419 u(ii) = sqrt(m_droptol) * rownorm;
420 m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii);
421
422 // sort the U-part of the row
423 // apply the dropping rule first
424 len = 0;
425 for(Index k = 1; k < sizeu; k++)
426 {
427 if(abs(u(ii+k)) > m_droptol * rownorm )
428 {
429 ++len;
430 u(ii + len) = u(ii + k);
431 ju(ii + len) = ju(ii + k);
432 }
433 }
434 sizeu = len + 1; // +1 to take into account the diagonal element
435 len = (std::min)(sizeu, nnzU);
436 typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
437 typename VectorI::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
438 internal::QuickSplit(uu, juu, len);
439
440 // store the largest elements of the U part
441 for(Index k = ii + 1; k < ii + len; k++)
442 m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
443 }
444 m_lu.finalize();
445 m_lu.makeCompressed();
446
447 m_factorizationIsOk = true;
448 m_info = Success;
449 }
450
451 } // end namespace Eigen
452
453 #endif // EIGEN_INCOMPLETE_LUT_H
454