Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
IncompleteLUT.h
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 <desire.nuentsa_wakam@inria.fr>
5// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@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
15namespace Eigen {
16
17namespace internal {
18
28template <typename VectorV, typename VectorI>
29Index 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
98template <typename _Scalar, typename _StorageIndex = int>
99class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar, _StorageIndex> >
100{
101 protected:
103 using Base::m_isInitialized;
104 public:
105 typedef _Scalar Scalar;
106 typedef _StorageIndex StorageIndex;
107 typedef typename NumTraits<Scalar>::Real RealScalar;
111
112 enum {
113 ColsAtCompileTime = Dynamic,
114 MaxColsAtCompileTime = Dynamic
115 };
116
117 public:
118
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 EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_lu.rows(); }
134
135 EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_lu.cols(); }
136
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
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
179protected:
180
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
189protected:
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
205template<typename Scalar, typename StorageIndex>
207{
208 this->m_droptol = droptol;
209}
210
215template<typename Scalar, typename StorageIndex>
217{
218 this->m_fillfactor = fillfactor;
219}
220
221template <typename Scalar, typename StorageIndex>
222template<typename _MatrixType>
223void 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.
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...
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
242template <typename Scalar, typename StorageIndex>
243template<typename _MatrixType>
244void 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
Definition: Ordering.h:51
Incomplete LU factorization with dual-threshold strategy.
Definition: IncompleteLUT.h:100
void setFillfactor(int fillfactor)
Definition: IncompleteLUT.h:216
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: IncompleteLUT.h:142
void setDroptol(const RealScalar &droptol)
Definition: IncompleteLUT.h:206
IncompleteLUT & compute(const MatrixType &amat)
Definition: IncompleteLUT.h:160
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Permutation matrix.
Definition: PermutationMatrix.h:298
A versatible sparse matrix representation.
Definition: SparseMatrix.h:98
Index rows() const
Definition: SparseMatrix.h:138
Index cols() const
Definition: SparseMatrix.h:140
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
static const symbolic::SymbolExpr< internal::symbolic_last_tag > last
Definition: IndexedViewHelper.h:38
ComputationInfo
Definition: Constants.h:440
@ NumericalIssue
Definition: Constants.h:444
@ Success
Definition: Constants.h:442
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const int Dynamic
Definition: Constants.h:22
Definition: IncompleteLUT.h:182
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:233