11 #ifndef EIGEN_INCOMPLETE_CHOlESKY_H
12 #define EIGEN_INCOMPLETE_CHOlESKY_H
44 template <
typename Scalar,
int _UpLo =
Lower,
typename _OrderingType =
45 #ifndef EIGEN_MPL2_ONLY
55 using Base::m_isInitialized;
58 typedef _OrderingType OrderingType;
59 typedef typename OrderingType::PermutationType PermutationType;
60 typedef typename PermutationType::StorageIndex StorageIndex;
65 typedef std::vector<std::list<StorageIndex> > VectorList;
66 enum { UpLo = _UpLo };
83 template<
typename MatrixType>
106 eigen_assert(m_isInitialized &&
"IncompleteCholesky is not initialized.");
116 template<
typename MatrixType>
120 PermutationType pinv;
121 ord(mat.template selfadjointView<UpLo>(), pinv);
122 if(pinv.size()>0) m_perm = pinv.inverse();
123 else m_perm.resize(0);
124 m_L.
resize(mat.rows(), mat.cols());
125 m_analysisIsOk =
true;
126 m_isInitialized =
true;
137 template<
typename MatrixType>
146 template<
typename MatrixType>
154 template<
typename Rhs,
typename Dest>
155 void _solve_impl(
const Rhs& b, Dest& x)
const
157 eigen_assert(m_factorizationIsOk &&
"factorize() should be called first");
158 if (m_perm.rows() == b.rows()) x = m_perm * b;
160 x = m_scale.asDiagonal() * x;
161 x = m_L.template triangularView<Lower>().solve(x);
162 x = m_L.adjoint().template triangularView<Upper>().solve(x);
163 x = m_scale.asDiagonal() * x;
164 if (m_perm.rows() == b.rows())
165 x = m_perm.inverse() * x;
175 const PermutationType&
permutationP()
const { eigen_assert(
"m_analysisIsOk");
return m_perm; }
180 RealScalar m_initialShift;
182 bool m_factorizationIsOk;
184 PermutationType m_perm;
194 template<
typename Scalar,
int _UpLo,
typename OrderingType>
195 template<
typename _MatrixType>
199 eigen_assert(m_analysisIsOk &&
"analyzePattern() should be called first");
204 if (m_perm.rows() == mat.rows() )
207 FactorType tmp(mat.rows(), mat.cols());
208 tmp = mat.template selfadjointView<_UpLo>().twistedBy(m_perm);
209 m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
213 m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
216 Index n = m_L.cols();
217 Index nnz = m_L.nonZeros();
218 Map<VectorSx> vals(m_L.valuePtr(), nnz);
219 Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz);
220 Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1);
221 VectorIx firstElt(n-1);
222 VectorList listCol(n);
223 VectorSx col_vals(n);
224 VectorIx col_irow(n);
225 VectorIx col_pattern(n);
226 col_pattern.fill(-1);
227 StorageIndex col_nnz;
233 for (
Index j = 0; j < n; j++)
234 for (
Index k = colPtr[j]; k < colPtr[j+1]; k++)
236 m_scale(j) += numext::abs2(vals(k));
238 m_scale(rowIdx[k]) += numext::abs2(vals(k));
241 m_scale = m_scale.cwiseSqrt().cwiseSqrt();
243 for (
Index j = 0; j < n; ++j)
244 if(m_scale(j)>(std::numeric_limits<RealScalar>::min)())
245 m_scale(j) = RealScalar(1)/m_scale(j);
252 RealScalar mindiag = NumTraits<RealScalar>::highest();
253 for (
Index j = 0; j < n; j++)
255 for (
Index k = colPtr[j]; k < colPtr[j+1]; k++)
256 vals[k] *= (m_scale(j)*m_scale(rowIdx[k]));
257 eigen_internal_assert(rowIdx[colPtr[j]]==j &&
"IncompleteCholesky: only the lower triangular part must be stored");
258 mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
261 FactorType L_save = m_L;
263 RealScalar shift = 0;
264 if(mindiag <= RealScalar(0.))
265 shift = m_initialShift - mindiag;
274 for (
Index j = 0; j < n; j++)
275 vals[colPtr[j]] += shift;
283 Scalar diag = vals[colPtr[j]];
285 for (
Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
287 StorageIndex l = rowIdx[i];
288 col_vals(col_nnz) = vals[i];
289 col_irow(col_nnz) = l;
290 col_pattern(l) = col_nnz;
294 typename std::list<StorageIndex>::iterator k;
296 for(k = listCol[j].begin(); k != listCol[j].end(); k++)
298 Index jk = firstElt(*k);
299 eigen_internal_assert(rowIdx[jk]==j);
300 Scalar v_j_jk = numext::conj(vals[jk]);
303 for (
Index i = jk; i < colPtr[*k+1]; i++)
305 StorageIndex l = rowIdx[i];
308 col_vals(col_nnz) = vals[i] * v_j_jk;
309 col_irow[col_nnz] = l;
310 col_pattern(l) = col_nnz;
314 col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
316 updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
321 if(numext::real(diag) <= 0)
327 shift = numext::maxi(m_initialShift,RealScalar(2)*shift);
329 vals = Map<const VectorSx>(L_save.valuePtr(), nnz);
330 rowIdx = Map<const VectorIx>(L_save.innerIndexPtr(), nnz);
331 colPtr = Map<const VectorIx>(L_save.outerIndexPtr(), n+1);
332 col_pattern.fill(-1);
333 for(
Index i=0; i<n; ++i)
339 RealScalar rdiag =
sqrt(numext::real(diag));
340 vals[colPtr[j]] = rdiag;
341 for (
Index k = 0; k<col_nnz; ++k)
343 Index i = col_irow[k];
345 col_vals(k) /= rdiag;
347 vals[colPtr[i]] -= numext::abs2(col_vals(k));
351 Index p = colPtr[j+1] - colPtr[j] - 1 ;
352 Ref<VectorSx> cvals = col_vals.head(col_nnz);
353 Ref<VectorIx> cirow = col_irow.head(col_nnz);
354 internal::QuickSplit(cvals,cirow, p);
357 for (
Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
359 vals[i] = col_vals(cpt);
360 rowIdx[i] = col_irow(cpt);
362 col_pattern(col_irow(cpt)) = -1;
366 Index jk = colPtr(j)+1;
367 updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
372 m_factorizationIsOk =
true;
378 template<
typename Scalar,
int _UpLo,
typename OrderingType>
379 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)
381 if (jk < colPtr(col+1) )
383 Index p = colPtr(col+1) - jk;
385 rowIdx.segment(jk,p).minCoeff(&minpos);
387 if (rowIdx(minpos) != rowIdx(jk))
390 std::swap(rowIdx(jk),rowIdx(minpos));
391 std::swap(vals(jk),vals(minpos));
393 firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
394 listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));