10 #ifndef EIGEN_INCOMPLETE_CHOlESKY_H
11 #define EIGEN_INCOMPLETE_CHOlESKY_H
12 #include "Eigen/src/IterativeLinearSolvers/IncompleteLUT.h"
13 #include <Eigen/OrderingMethods>
29 template <
typename Scalar,
int _UpLo = Lower,
typename _OrderingType = AMDOrdering<
int> >
30 class IncompleteCholesky :
public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> >
33 typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base;
34 using Base::m_isInitialized;
36 typedef typename NumTraits<Scalar>::Real RealScalar;
37 typedef _OrderingType OrderingType;
38 typedef typename OrderingType::PermutationType PermutationType;
39 typedef typename PermutationType::StorageIndex StorageIndex;
40 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> FactorType;
41 typedef FactorType MatrixType;
42 typedef Matrix<Scalar,Dynamic,1> VectorSx;
43 typedef Matrix<RealScalar,Dynamic,1> VectorRx;
44 typedef Matrix<StorageIndex,Dynamic, 1> VectorIx;
45 typedef std::vector<std::list<StorageIndex> > VectorList;
46 enum { UpLo = _UpLo };
50 template<
typename MatrixType>
51 IncompleteCholesky(
const MatrixType& matrix) : m_initialShift(1e-3),m_factorizationIsOk(
false)
56 Index rows()
const {
return m_L.rows(); }
58 Index cols()
const {
return m_L.cols(); }
66 ComputationInfo
info()
const
68 eigen_assert(m_isInitialized &&
"IncompleteLLT is not initialized.");
80 template<
typename MatrixType>
85 ord(mat.template selfadjointView<UpLo>(), pinv);
86 if(pinv.size()>0) m_perm = pinv.inverse();
87 else m_perm.resize(0);
88 m_analysisIsOk =
true;
91 template<
typename MatrixType>
92 void factorize(
const MatrixType& amat);
94 template<
typename MatrixType>
95 void compute(
const MatrixType& matrix)
101 template<
typename Rhs,
typename Dest>
102 void _solve_impl(
const Rhs& b, Dest& x)
const
104 eigen_assert(m_factorizationIsOk &&
"factorize() should be called first");
105 if (m_perm.rows() == b.rows()) x = m_perm * b;
107 x = m_scale.asDiagonal() * x;
108 x = m_L.template triangularView<Lower>().solve(x);
109 x = m_L.adjoint().template triangularView<Upper>().solve(x);
110 x = m_scale.asDiagonal() * x;
111 if (m_perm.rows() == b.rows())
112 x = m_perm.inverse() * x;
119 RealScalar m_initialShift;
121 bool m_factorizationIsOk;
122 ComputationInfo m_info;
123 PermutationType m_perm;
126 inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals,
const Index& col,
const Index& jk, VectorIx& firstElt, VectorList& listCol);
129 template<
typename Scalar,
int _UpLo,
typename OrderingType>
130 template<
typename _MatrixType>
131 void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(
const _MatrixType& mat)
134 eigen_assert(m_analysisIsOk &&
"analyzePattern() should be called first");
138 m_L.resize(mat.rows(), mat.cols());
141 if (m_perm.rows() == mat.rows() )
144 FactorType tmp(mat.rows(), mat.cols());
145 tmp = mat.template selfadjointView<_UpLo>().twistedBy(m_perm);
146 m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
150 m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
153 Index n = m_L.cols();
154 Index nnz = m_L.nonZeros();
155 Map<VectorSx> vals(m_L.valuePtr(), nnz);
156 Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz);
157 Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1);
158 VectorIx firstElt(n-1);
159 VectorList listCol(n);
160 VectorSx col_vals(n);
161 VectorIx col_irow(n);
162 VectorIx col_pattern(n);
163 col_pattern.fill(-1);
164 StorageIndex col_nnz;
170 for (Index j = 0; j < n; j++)
171 for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
173 m_scale(j) += numext::abs2(vals(k));
175 m_scale(rowIdx[k]) += numext::abs2(vals(k));
178 m_scale = m_scale.cwiseSqrt().cwiseSqrt();
181 RealScalar mindiag = NumTraits<RealScalar>::highest();
182 for (Index j = 0; j < n; j++)
184 for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
185 vals[k] /= (m_scale(j)*m_scale(rowIdx[k]));
186 eigen_internal_assert(rowIdx[colPtr[j]]==j &&
"IncompleteCholesky: only the lower triangular part must be stored");
187 mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
190 RealScalar shift = 0;
191 if(mindiag <= RealScalar(0.))
192 shift = m_initialShift - mindiag;
195 for (Index j = 0; j < n; j++)
196 vals[colPtr[j]] += shift;
199 for (Index j=0; j < n; ++j)
203 Scalar diag = vals[colPtr[j]];
205 for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
207 StorageIndex l = rowIdx[i];
208 col_vals(col_nnz) = vals[i];
209 col_irow(col_nnz) = l;
210 col_pattern(l) = col_nnz;
214 typename std::list<StorageIndex>::iterator k;
216 for(k = listCol[j].begin(); k != listCol[j].end(); k++)
218 Index jk = firstElt(*k);
219 eigen_internal_assert(rowIdx[jk]==j);
220 Scalar v_j_jk = numext::conj(vals[jk]);
223 for (Index i = jk; i < colPtr[*k+1]; i++)
225 StorageIndex l = rowIdx[i];
228 col_vals(col_nnz) = vals[i] * v_j_jk;
229 col_irow[col_nnz] = l;
230 col_pattern(l) = col_nnz;
234 col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
236 updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
241 if(numext::real(diag) <= 0)
243 std::cerr <<
"\nNegative diagonal during Incomplete factorization at position " << j <<
" (value = " << diag <<
")\n";
244 m_info = NumericalIssue;
248 RealScalar rdiag = sqrt(numext::real(diag));
249 vals[colPtr[j]] = rdiag;
250 for (Index k = 0; k<col_nnz; ++k)
252 Index i = col_irow[k];
254 col_vals(k) /= rdiag;
256 vals[colPtr[i]] -= numext::abs2(col_vals(k));
260 Index p = colPtr[j+1] - colPtr[j] - 1 ;
261 Ref<VectorSx> cvals = col_vals.head(col_nnz);
262 Ref<VectorIx> cirow = col_irow.head(col_nnz);
263 internal::QuickSplit(cvals,cirow, p);
266 for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
268 vals[i] = col_vals(cpt);
269 rowIdx[i] = col_irow(cpt);
271 col_pattern(col_irow(cpt)) = -1;
275 Index jk = colPtr(j)+1;
276 updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
278 m_factorizationIsOk =
true;
279 m_isInitialized =
true;
283 template<
typename Scalar,
int _UpLo,
typename OrderingType>
284 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)
286 if (jk < colPtr(col+1) )
288 Index p = colPtr(col+1) - jk;
290 rowIdx.segment(jk,p).minCoeff(&minpos);
292 if (rowIdx(minpos) != rowIdx(jk))
295 std::swap(rowIdx(jk),rowIdx(minpos));
296 std::swap(vals(jk),vals(minpos));
298 firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
299 listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
Modified Incomplete Cholesky with dual threshold.
Definition: IncompleteCholesky.h:30
Namespace containing all symbols from the Eigen library.
Definition: CXX11Meta.h:13
void setInitialShift(RealScalar shift)
Set the initial shift parameter.
Definition: IncompleteCholesky.h:75
void analyzePattern(const MatrixType &mat)
Computes the fill reducing permutation vector.
Definition: IncompleteCholesky.h:81
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: IncompleteCholesky.h:66