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
IncompleteCholesky.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) 2015 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_CHOlESKY_H
12#define EIGEN_INCOMPLETE_CHOlESKY_H
13
14#include <vector>
15#include <list>
16
17namespace Eigen {
44template <typename Scalar, int _UpLo = Lower, typename _OrderingType = AMDOrdering<int> >
45class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> >
46{
47 protected:
49 using Base::m_isInitialized;
50 public:
51 typedef typename NumTraits<Scalar>::Real RealScalar;
52 typedef _OrderingType OrderingType;
53 typedef typename OrderingType::PermutationType PermutationType;
54 typedef typename PermutationType::StorageIndex StorageIndex;
59 typedef std::vector<std::list<StorageIndex> > VectorList;
60 enum { UpLo = _UpLo };
61 enum {
62 ColsAtCompileTime = Dynamic,
63 MaxColsAtCompileTime = Dynamic
64 };
65 public:
66
73 IncompleteCholesky() : m_initialShift(1e-3),m_analysisIsOk(false),m_factorizationIsOk(false) {}
74
77 template<typename MatrixType>
78 IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_analysisIsOk(false),m_factorizationIsOk(false)
79 {
80 compute(matrix);
81 }
82
84 EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_L.rows(); }
85
87 EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_L.cols(); }
88
89
99 {
100 eigen_assert(m_isInitialized && "IncompleteCholesky is not initialized.");
101 return m_info;
102 }
103
106 void setInitialShift(RealScalar shift) { m_initialShift = shift; }
107
110 template<typename MatrixType>
111 void analyzePattern(const MatrixType& mat)
112 {
113 OrderingType ord;
114 PermutationType pinv;
115 ord(mat.template selfadjointView<UpLo>(), pinv);
116 if(pinv.size()>0) m_perm = pinv.inverse();
117 else m_perm.resize(0);
118 m_L.resize(mat.rows(), mat.cols());
119 m_analysisIsOk = true;
120 m_isInitialized = true;
121 m_info = Success;
122 }
123
131 template<typename MatrixType>
132 void factorize(const MatrixType& mat);
133
140 template<typename MatrixType>
141 void compute(const MatrixType& mat)
142 {
143 analyzePattern(mat);
144 factorize(mat);
145 }
146
147 // internal
148 template<typename Rhs, typename Dest>
149 void _solve_impl(const Rhs& b, Dest& x) const
150 {
151 eigen_assert(m_factorizationIsOk && "factorize() should be called first");
152 if (m_perm.rows() == b.rows()) x = m_perm * b;
153 else x = b;
154 x = m_scale.asDiagonal() * x;
155 x = m_L.template triangularView<Lower>().solve(x);
156 x = m_L.adjoint().template triangularView<Upper>().solve(x);
157 x = m_scale.asDiagonal() * x;
158 if (m_perm.rows() == b.rows())
159 x = m_perm.inverse() * x;
160 }
161
163 const FactorType& matrixL() const { eigen_assert("m_factorizationIsOk"); return m_L; }
164
166 const VectorRx& scalingS() const { eigen_assert("m_factorizationIsOk"); return m_scale; }
167
169 const PermutationType& permutationP() const { eigen_assert("m_analysisIsOk"); return m_perm; }
170
171 protected:
172 FactorType m_L; // The lower part stored in CSC
173 VectorRx m_scale; // The vector for scaling the matrix
174 RealScalar m_initialShift; // The initial shift parameter
175 bool m_analysisIsOk;
176 bool m_factorizationIsOk;
177 ComputationInfo m_info;
178 PermutationType m_perm;
179
180 private:
181 inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol);
182};
183
184// Based on the following paper:
185// C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with
186// Limited memory, SIAM J. Sci. Comput. 21(1), pp. 24-45, 1999
187// http://ftp.mcs.anl.gov/pub/tech_reports/reports/P682.pdf
188template<typename Scalar, int _UpLo, typename OrderingType>
189template<typename _MatrixType>
190void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat)
191{
192 using std::sqrt;
193 eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
194
195 // Dropping strategy : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added
196
197 // Apply the fill-reducing permutation computed in analyzePattern()
198 if (m_perm.rows() == mat.rows() ) // To detect the null permutation
199 {
200 // The temporary is needed to make sure that the diagonal entry is properly sorted
201 FactorType tmp(mat.rows(), mat.cols());
202 tmp = mat.template selfadjointView<_UpLo>().twistedBy(m_perm);
203 m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
204 }
205 else
206 {
207 m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
208 }
209
210 Index n = m_L.cols();
211 Index nnz = m_L.nonZeros();
212 Map<VectorSx> vals(m_L.valuePtr(), nnz); //values
213 Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz); //Row indices
214 Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row
215 VectorIx firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
216 VectorList listCol(n); // listCol(j) is a linked list of columns to update column j
217 VectorSx col_vals(n); // Store a nonzero values in each column
218 VectorIx col_irow(n); // Row indices of nonzero elements in each column
219 VectorIx col_pattern(n);
220 col_pattern.fill(-1);
221 StorageIndex col_nnz;
222
223
224 // Computes the scaling factors
225 m_scale.resize(n);
226 m_scale.setZero();
227 for (Index j = 0; j < n; j++)
228 for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
229 {
230 m_scale(j) += numext::abs2(vals(k));
231 if(rowIdx[k]!=j)
232 m_scale(rowIdx[k]) += numext::abs2(vals(k));
233 }
234
235 m_scale = m_scale.cwiseSqrt().cwiseSqrt();
236
237 for (Index j = 0; j < n; ++j)
238 if(m_scale(j)>(std::numeric_limits<RealScalar>::min)())
239 m_scale(j) = RealScalar(1)/m_scale(j);
240 else
241 m_scale(j) = 1;
242
243 // TODO disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster)
244
245 // Scale and compute the shift for the matrix
246 RealScalar mindiag = NumTraits<RealScalar>::highest();
247 for (Index j = 0; j < n; j++)
248 {
249 for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
250 vals[k] *= (m_scale(j)*m_scale(rowIdx[k]));
251 eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored");
252 mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
253 }
254
255 FactorType L_save = m_L;
256
257 RealScalar shift = 0;
258 if(mindiag <= RealScalar(0.))
259 shift = m_initialShift - mindiag;
260
261 m_info = NumericalIssue;
262
263 // Try to perform the incomplete factorization using the current shift
264 int iter = 0;
265 do
266 {
267 // Apply the shift to the diagonal elements of the matrix
268 for (Index j = 0; j < n; j++)
269 vals[colPtr[j]] += shift;
270
271 // jki version of the Cholesky factorization
272 Index j=0;
273 for (; j < n; ++j)
274 {
275 // Left-looking factorization of the j-th column
276 // First, load the j-th column into col_vals
277 Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored
278 col_nnz = 0;
279 for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
280 {
281 StorageIndex l = rowIdx[i];
282 col_vals(col_nnz) = vals[i];
283 col_irow(col_nnz) = l;
284 col_pattern(l) = col_nnz;
285 col_nnz++;
286 }
287 {
288 typename std::list<StorageIndex>::iterator k;
289 // Browse all previous columns that will update column j
290 for(k = listCol[j].begin(); k != listCol[j].end(); k++)
291 {
292 Index jk = firstElt(*k); // First element to use in the column
293 eigen_internal_assert(rowIdx[jk]==j);
294 Scalar v_j_jk = numext::conj(vals[jk]);
295
296 jk += 1;
297 for (Index i = jk; i < colPtr[*k+1]; i++)
298 {
299 StorageIndex l = rowIdx[i];
300 if(col_pattern[l]<0)
301 {
302 col_vals(col_nnz) = vals[i] * v_j_jk;
303 col_irow[col_nnz] = l;
304 col_pattern(l) = col_nnz;
305 col_nnz++;
306 }
307 else
308 col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
309 }
310 updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
311 }
312 }
313
314 // Scale the current column
315 if(numext::real(diag) <= 0)
316 {
317 if(++iter>=10)
318 return;
319
320 // increase shift
321 shift = numext::maxi(m_initialShift,RealScalar(2)*shift);
322 // restore m_L, col_pattern, and listCol
323 vals = Map<const VectorSx>(L_save.valuePtr(), nnz);
324 rowIdx = Map<const VectorIx>(L_save.innerIndexPtr(), nnz);
325 colPtr = Map<const VectorIx>(L_save.outerIndexPtr(), n+1);
326 col_pattern.fill(-1);
327 for(Index i=0; i<n; ++i)
328 listCol[i].clear();
329
330 break;
331 }
332
333 RealScalar rdiag = sqrt(numext::real(diag));
334 vals[colPtr[j]] = rdiag;
335 for (Index k = 0; k<col_nnz; ++k)
336 {
337 Index i = col_irow[k];
338 //Scale
339 col_vals(k) /= rdiag;
340 //Update the remaining diagonals with col_vals
341 vals[colPtr[i]] -= numext::abs2(col_vals(k));
342 }
343 // Select the largest p elements
344 // p is the original number of elements in the column (without the diagonal)
345 Index p = colPtr[j+1] - colPtr[j] - 1 ;
346 Ref<VectorSx> cvals = col_vals.head(col_nnz);
347 Ref<VectorIx> cirow = col_irow.head(col_nnz);
348 internal::QuickSplit(cvals,cirow, p);
349 // Insert the largest p elements in the matrix
350 Index cpt = 0;
351 for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
352 {
353 vals[i] = col_vals(cpt);
354 rowIdx[i] = col_irow(cpt);
355 // restore col_pattern:
356 col_pattern(col_irow(cpt)) = -1;
357 cpt++;
358 }
359 // Get the first smallest row index and put it after the diagonal element
360 Index jk = colPtr(j)+1;
361 updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
362 }
363
364 if(j==n)
365 {
366 m_factorizationIsOk = true;
367 m_info = Success;
368 }
369 } while(m_info!=Success);
370}
371
372template<typename Scalar, int _UpLo, typename OrderingType>
373inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol)
374{
375 if (jk < colPtr(col+1) )
376 {
377 Index p = colPtr(col+1) - jk;
378 Index minpos;
379 rowIdx.segment(jk,p).minCoeff(&minpos);
380 minpos += jk;
381 if (rowIdx(minpos) != rowIdx(jk))
382 {
383 //Swap
384 std::swap(rowIdx(jk),rowIdx(minpos));
385 std::swap(vals(jk),vals(minpos));
386 }
387 firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
388 listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
389 }
390}
391
392} // end namespace Eigen
393
394#endif
Modified Incomplete Cholesky with dual threshold.
Definition: IncompleteCholesky.h:46
void setInitialShift(RealScalar shift)
Set the initial shift parameter .
Definition: IncompleteCholesky.h:106
void analyzePattern(const MatrixType &mat)
Computes the fill reducing permutation vector using the sparsity pattern of mat.
Definition: IncompleteCholesky.h:111
IncompleteCholesky(const MatrixType &matrix)
Definition: IncompleteCholesky.h:78
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: IncompleteCholesky.h:84
const VectorRx & scalingS() const
Definition: IncompleteCholesky.h:166
void compute(const MatrixType &mat)
Definition: IncompleteCholesky.h:141
const PermutationType & permutationP() const
Definition: IncompleteCholesky.h:169
void factorize(const MatrixType &mat)
Performs the numerical factorization of the input matrix mat.
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: IncompleteCholesky.h:87
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: IncompleteCholesky.h:98
IncompleteCholesky()
Definition: IncompleteCholesky.h:73
const FactorType & matrixL() const
Definition: IncompleteCholesky.h:163
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:283
A versatible sparse matrix representation.
Definition: SparseMatrix.h:98
Index rows() const
Definition: SparseMatrix.h:138
Index cols() const
Definition: SparseMatrix.h:140
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:626
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
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