11 #ifndef EIGEN_INCOMPLETE_CHOlESKY_H 12 #define EIGEN_INCOMPLETE_CHOlESKY_H 42 template <
typename Scalar,
int _UpLo =
Lower,
typename _OrderingType =
43 #ifndef EIGEN_MPL2_ONLY 53 using Base::m_isInitialized;
56 typedef _OrderingType OrderingType;
57 typedef typename OrderingType::PermutationType PermutationType;
58 typedef typename PermutationType::StorageIndex StorageIndex;
63 typedef std::vector<std::list<StorageIndex> > VectorList;
64 enum { UpLo = _UpLo };
66 ColsAtCompileTime = Dynamic,
67 MaxColsAtCompileTime = Dynamic
81 template<
typename MatrixType>
104 eigen_assert(m_isInitialized &&
"IncompleteCholesky is not initialized.");
114 template<
typename MatrixType>
118 PermutationType pinv;
119 ord(mat.template selfadjointView<UpLo>(), pinv);
120 if(pinv.size()>0) m_perm = pinv.inverse();
121 else m_perm.resize(0);
122 m_L.
resize(mat.rows(), mat.cols());
123 m_analysisIsOk =
true;
124 m_isInitialized =
true;
135 template<
typename MatrixType>
144 template<
typename MatrixType>
152 template<
typename Rhs,
typename Dest>
153 void _solve_impl(
const Rhs& b, Dest& x)
const 155 eigen_assert(m_factorizationIsOk &&
"factorize() should be called first");
156 if (m_perm.rows() == b.rows()) x = m_perm * b;
158 x = m_scale.asDiagonal() * x;
159 x = m_L.template triangularView<Lower>().
solve(x);
160 x = m_L.adjoint().template triangularView<Upper>().
solve(x);
161 x = m_scale.asDiagonal() * x;
162 if (m_perm.rows() == b.rows())
163 x = m_perm.inverse() * x;
167 const FactorType&
matrixL()
const { eigen_assert(
"m_factorizationIsOk");
return m_L; }
170 const VectorRx&
scalingS()
const { eigen_assert(
"m_factorizationIsOk");
return m_scale; }
173 const PermutationType&
permutationP()
const { eigen_assert(
"m_analysisIsOk");
return m_perm; }
178 RealScalar m_initialShift;
180 bool m_factorizationIsOk;
182 PermutationType m_perm;
188 template<
typename Scalar,
int _UpLo,
typename OrderingType>
189 template<
typename _MatrixType>
193 eigen_assert(m_analysisIsOk &&
"analyzePattern() should be called first");
198 if (m_perm.rows() == mat.rows() )
202 tmp = mat.template selfadjointView<_UpLo>().twistedBy(m_perm);
203 m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
207 m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
210 Index n = m_L.
cols();
211 Index nnz = m_L.nonZeros();
216 VectorList listCol(n);
220 col_pattern.fill(-1);
221 StorageIndex col_nnz;
227 for (Index j = 0; j < n; j++)
228 for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
230 m_scale(j) += numext::abs2(vals(k));
232 m_scale(rowIdx[k]) += numext::abs2(vals(k));
235 m_scale = m_scale.cwiseSqrt().cwiseSqrt();
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);
247 for (Index j = 0; j < n; j++)
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);
255 RealScalar shift = 0;
256 if(mindiag <= RealScalar(0.))
257 shift = m_initialShift - mindiag;
260 for (Index j = 0; j < n; j++)
261 vals[colPtr[j]] += shift;
264 for (Index j=0; j < n; ++j)
268 Scalar diag = vals[colPtr[j]];
270 for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
272 StorageIndex l = rowIdx[i];
273 col_vals(col_nnz) = vals[i];
274 col_irow(col_nnz) = l;
275 col_pattern(l) = col_nnz;
279 typename std::list<StorageIndex>::iterator k;
281 for(k = listCol[j].begin(); k != listCol[j].end(); k++)
283 Index jk = firstElt(*k);
284 eigen_internal_assert(rowIdx[jk]==j);
285 Scalar v_j_jk = numext::conj(vals[jk]);
288 for (Index i = jk; i < colPtr[*k+1]; i++)
290 StorageIndex l = rowIdx[i];
293 col_vals(col_nnz) = vals[i] * v_j_jk;
294 col_irow[col_nnz] = l;
295 col_pattern(l) = col_nnz;
299 col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
301 updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
306 if(numext::real(diag) <= 0)
312 RealScalar rdiag = sqrt(numext::real(diag));
313 vals[colPtr[j]] = rdiag;
314 for (Index k = 0; k<col_nnz; ++k)
316 Index i = col_irow[k];
318 col_vals(k) /= rdiag;
320 vals[colPtr[i]] -= numext::abs2(col_vals(k));
324 Index p = colPtr[j+1] - colPtr[j] - 1 ;
327 internal::QuickSplit(cvals,cirow, p);
330 for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
332 vals[i] = col_vals(cpt);
333 rowIdx[i] = col_irow(cpt);
335 col_pattern(col_irow(cpt)) = -1;
339 Index jk = colPtr(j)+1;
340 updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
342 m_factorizationIsOk =
true;
346 template<
typename Scalar,
int _UpLo,
typename OrderingType>
349 if (jk < colPtr(col+1) )
351 Index p = colPtr(col+1) - jk;
353 rowIdx.segment(jk,p).minCoeff(&minpos);
355 if (rowIdx(minpos) != rowIdx(jk))
358 std::swap(rowIdx(jk),rowIdx(minpos));
359 std::swap(vals(jk),vals(minpos));
361 firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
362 listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
Index cols() const
Definition: IncompleteCholesky.h:91
void factorize(const MatrixType &mat)
Performs the numerical factorization of the input matrix mat.
Index cols() const
Definition: SparseMatrix.h:133
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:89
Modified Incomplete Cholesky with dual threshold.
Definition: IncompleteCholesky.h:49
const PermutationType & permutationP() const
Definition: IncompleteCholesky.h:173
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
IncompleteCholesky(const MatrixType &matrix)
Definition: IncompleteCholesky.h:82
const Solve< IncompleteCholesky< Scalar, _UpLo, _OrderingType >, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SparseSolverBase.h:74
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:161
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:252
IncompleteCholesky()
Definition: IncompleteCholesky.h:77
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:520
Definition: Constants.h:204
Definition: Constants.h:434
const VectorRx & scalingS() const
Definition: IncompleteCholesky.h:170
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:611
void setInitialShift(RealScalar shift)
Set the initial shift parameter .
Definition: IncompleteCholesky.h:110
Definition: Constants.h:432
Index rows() const
Definition: IncompleteCholesky.h:88
const StorageIndex * innerIndexPtr() const
Definition: SparseMatrix.h:152
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:186
void analyzePattern(const MatrixType &mat)
Computes the fill reducing permutation vector using the sparsity pattern of mat.
Definition: IncompleteCholesky.h:115
const FactorType & matrixL() const
Definition: IncompleteCholesky.h:167
void compute(const MatrixType &mat)
Definition: IncompleteCholesky.h:145
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: IncompleteCholesky.h:102
const Scalar * valuePtr() const
Definition: SparseMatrix.h:143
ComputationInfo
Definition: Constants.h:430
Index rows() const
Definition: SparseMatrix.h:131