11 #ifndef EIGEN_INCOMPLETE_CHOlESKY_H
12 #define EIGEN_INCOMPLETE_CHOlESKY_H
17 #include "./InternalHeaderCheck.h"
46 template <
typename Scalar,
int UpLo_ = Lower,
typename OrderingType_ = AMDOrdering<
int> >
51 using Base::m_isInitialized;
54 typedef OrderingType_ OrderingType;
55 typedef typename OrderingType::PermutationType PermutationType;
56 typedef typename PermutationType::StorageIndex StorageIndex;
61 typedef std::vector<std::list<StorageIndex> > VectorList;
62 enum { UpLo = UpLo_ };
79 template<
typename MatrixType>
80 IncompleteCholesky(
const MatrixType& matrix) : m_initialShift(1e-3),m_analysisIsOk(false),m_factorizationIsOk(false)
86 EIGEN_CONSTEXPR
Index rows() const EIGEN_NOEXCEPT {
return m_L.
rows(); }
89 EIGEN_CONSTEXPR
Index cols() const EIGEN_NOEXCEPT {
return m_L.
cols(); }
102 eigen_assert(m_isInitialized &&
"IncompleteCholesky is not initialized.");
112 template<
typename MatrixType>
116 PermutationType pinv;
117 ord(mat.template selfadjointView<UpLo>(), pinv);
118 if(pinv.size()>0) m_perm = pinv.inverse();
119 else m_perm.resize(0);
120 m_L.
resize(mat.rows(), mat.cols());
121 m_analysisIsOk =
true;
122 m_isInitialized =
true;
133 template<
typename MatrixType>
142 template<
typename MatrixType>
150 template<
typename Rhs,
typename Dest>
151 void _solve_impl(
const Rhs& b, Dest& x)
const
153 eigen_assert(m_factorizationIsOk &&
"factorize() should be called first");
154 if (m_perm.rows() == b.rows()) x = m_perm * b;
156 x = m_scale.asDiagonal() * x;
157 x = m_L.template triangularView<Lower>().solve(x);
158 x = m_L.adjoint().template triangularView<Upper>().solve(x);
159 x = m_scale.asDiagonal() * x;
160 if (m_perm.rows() == b.rows())
161 x = m_perm.inverse() * x;
165 const FactorType&
matrixL()
const { eigen_assert(m_factorizationIsOk &&
"factorize() should be called first");
return m_L; }
168 const VectorRx&
scalingS()
const { eigen_assert(m_factorizationIsOk &&
"factorize() should be called first");
return m_scale; }
171 const PermutationType&
permutationP()
const { eigen_assert(m_analysisIsOk &&
"analyzePattern() should be called first");
return m_perm; }
176 RealScalar m_initialShift;
178 bool m_factorizationIsOk;
180 PermutationType m_perm;
190 template<
typename Scalar,
int UpLo_,
typename OrderingType>
191 template<
typename MatrixType_>
195 eigen_assert(m_analysisIsOk &&
"analyzePattern() should be called first");
200 if (m_perm.rows() == mat.rows() )
203 FactorType tmp(mat.rows(), mat.cols());
204 tmp = mat.template selfadjointView<UpLo_>().twistedBy(m_perm);
205 m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
209 m_L.template selfadjointView<Lower>() = mat.template selfadjointView<UpLo_>();
212 Index n = m_L.cols();
213 Index nnz = m_L.nonZeros();
214 Map<VectorSx> vals(m_L.valuePtr(), nnz);
215 Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz);
216 Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1);
217 VectorIx firstElt(n-1);
218 VectorList listCol(n);
219 VectorSx col_vals(n);
220 VectorIx col_irow(n);
221 VectorIx col_pattern(n);
222 col_pattern.fill(-1);
223 StorageIndex col_nnz;
229 for (
Index j = 0; j < n; j++)
230 for (
Index k = colPtr[j]; k < colPtr[j+1]; k++)
232 m_scale(j) += numext::abs2(vals(k));
234 m_scale(rowIdx[k]) += numext::abs2(vals(k));
237 m_scale = m_scale.cwiseSqrt().cwiseSqrt();
239 for (
Index j = 0; j < n; ++j)
240 if(m_scale(j)>(std::numeric_limits<RealScalar>::min)())
241 m_scale(j) = RealScalar(1)/m_scale(j);
248 RealScalar mindiag = NumTraits<RealScalar>::highest();
249 for (
Index j = 0; j < n; j++)
251 for (
Index k = colPtr[j]; k < colPtr[j+1]; k++)
252 vals[k] *= (m_scale(j)*m_scale(rowIdx[k]));
253 eigen_internal_assert(rowIdx[colPtr[j]]==j &&
"IncompleteCholesky: only the lower triangular part must be stored");
254 mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
257 FactorType L_save = m_L;
259 RealScalar shift = 0;
260 if(mindiag <= RealScalar(0.))
261 shift = m_initialShift - mindiag;
270 for (
Index j = 0; j < n; j++)
271 vals[colPtr[j]] += shift;
279 Scalar diag = vals[colPtr[j]];
281 for (
Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
283 StorageIndex l = rowIdx[i];
284 col_vals(col_nnz) = vals[i];
285 col_irow(col_nnz) = l;
286 col_pattern(l) = col_nnz;
290 typename std::list<StorageIndex>::iterator k;
292 for(k = listCol[j].begin(); k != listCol[j].end(); k++)
294 Index jk = firstElt(*k);
295 eigen_internal_assert(rowIdx[jk]==j);
296 Scalar v_j_jk = numext::conj(vals[jk]);
299 for (
Index i = jk; i < colPtr[*k+1]; i++)
301 StorageIndex l = rowIdx[i];
304 col_vals(col_nnz) = vals[i] * v_j_jk;
305 col_irow[col_nnz] = l;
306 col_pattern(l) = col_nnz;
310 col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
312 updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
317 if(numext::real(diag) <= 0)
323 shift = numext::maxi(m_initialShift,RealScalar(2)*shift);
325 vals = Map<const VectorSx>(L_save.valuePtr(), nnz);
326 rowIdx = Map<const VectorIx>(L_save.innerIndexPtr(), nnz);
327 colPtr = Map<const VectorIx>(L_save.outerIndexPtr(), n+1);
328 col_pattern.fill(-1);
329 for(
Index i=0; i<n; ++i)
335 RealScalar rdiag =
sqrt(numext::real(diag));
336 vals[colPtr[j]] = rdiag;
337 for (
Index k = 0; k<col_nnz; ++k)
339 Index i = col_irow[k];
341 col_vals(k) /= rdiag;
343 vals[colPtr[i]] -= numext::abs2(col_vals(k));
347 Index p = colPtr[j+1] - colPtr[j] - 1 ;
348 Ref<VectorSx> cvals = col_vals.head(col_nnz);
349 Ref<VectorIx> cirow = col_irow.head(col_nnz);
350 internal::QuickSplit(cvals,cirow, p);
353 for (
Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
355 vals[i] = col_vals(cpt);
356 rowIdx[i] = col_irow(cpt);
358 col_pattern(col_irow(cpt)) = -1;
362 Index jk = colPtr(j)+1;
363 updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
368 m_factorizationIsOk =
true;
374 template<
typename Scalar,
int UpLo_,
typename OrderingType>
375 inline 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)
377 if (jk < colPtr(col+1) )
379 Index p = colPtr(col+1) - jk;
381 rowIdx.segment(jk,p).minCoeff(&minpos);
383 if (rowIdx(minpos) != rowIdx(jk))
386 std::swap(rowIdx(jk),rowIdx(minpos));
387 std::swap(vals(jk),vals(minpos));
389 firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
390 listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
Modified Incomplete Cholesky with dual threshold.
Definition: IncompleteCholesky.h:48
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: IncompleteCholesky.h:100
IncompleteCholesky(const MatrixType &matrix)
Definition: IncompleteCholesky.h:80
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: IncompleteCholesky.h:86
void factorize(const MatrixType &mat)
Performs the numerical factorization of the input matrix mat.
const VectorRx & scalingS() const
Definition: IncompleteCholesky.h:168
IncompleteCholesky()
Definition: IncompleteCholesky.h:75
const PermutationType & permutationP() const
Definition: IncompleteCholesky.h:171
void compute(const MatrixType &mat)
Definition: IncompleteCholesky.h:143
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: IncompleteCholesky.h:89
const FactorType & matrixL() const
Definition: IncompleteCholesky.h:165
void analyzePattern(const MatrixType &mat)
Computes the fill reducing permutation vector using the sparsity pattern of mat.
Definition: IncompleteCholesky.h:113
void setInitialShift(RealScalar shift)
Set the initial shift parameter .
Definition: IncompleteCholesky.h:108
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:285
Index cols() const
Definition: SparseMatrix.h:142
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:626
Index rows() const
Definition: SparseMatrix.h:140
A base class for sparse solvers.
Definition: SparseSolverBase.h:70
ComputationInfo
Definition: Constants.h:442
@ NumericalIssue
Definition: Constants.h:446
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const int Dynamic
Definition: Constants.h:24
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231