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 65 typedef std::vector<std::list<StorageIndex> >
VectorList;
83 template<
typename MatrixType>
116 template<
typename MatrixType>
120 PermutationType pinv;
121 ord(mat.template selfadjointView<UpLo>(), pinv);
122 if(pinv.size()>0)
m_perm = pinv.inverse();
137 template<
typename MatrixType>
146 template<
typename MatrixType>
154 template<
typename Rhs,
typename Dest>
164 if (
m_perm.rows() == b.rows())
194 template<
typename Scalar,
int _UpLo,
typename OrderingType>
195 template<
typename _MatrixType>
204 if (
m_perm.rows() == mat.rows() )
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>();
226 col_pattern.fill(-1);
233 for (
Index j = 0; j <
n; j++)
234 for (
Index k = colPtr[j]; k < colPtr[j+1]; k++)
243 for (
Index j = 0; j <
n; ++j)
253 for (
Index j = 0; j <
n; j++)
255 for (
Index k = colPtr[j]; k < colPtr[j+1]; k++)
257 eigen_internal_assert(rowIdx[colPtr[j]]==j &&
"IncompleteCholesky: only the lower triangular part must be stored");
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++)
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);
303 for (
Index i = jk; i < colPtr[*k+1]; 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);
332 col_pattern.fill(-1);
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;
351 Index p = colPtr[j+1] - colPtr[j] - 1 ;
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);
378 template<
typename Scalar,
int _UpLo,
typename OrderingType>
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))
393 firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
394 listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
void updateList(Ref< const VectorIx > colPtr, Ref< VectorIx > rowIdx, Ref< VectorSx > vals, const Index &col, const Index &jk, VectorIx &firstElt, VectorList &listCol)
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
ComputationInfo info() const
Reports whether previous computation was successful.
const FactorType & matrixL() const
void factorize(const MatrixType &mat)
Performs the numerical factorization of the input matrix mat.
const StorageIndex * outerIndexPtr() const
PermutationType::StorageIndex StorageIndex
SparseSolverBase< IncompleteCholesky< Scalar, _UpLo, _OrderingType > > Base
A matrix or vector expression mapping an existing array of data.
Modified Incomplete Cholesky with dual threshold.
void resize(Index rows, Index cols)
EIGEN_DEVICE_FUNC ColXpr col(Index i)
This is the const version of col().
A base class for sparse solvers.
Matrix< StorageIndex, Dynamic, 1 > VectorIx
IncompleteCholesky(const MatrixType &matrix)
const Solve< IncompleteCholesky< Scalar, _UpLo, _OrderingType >, Rhs > solve(const MatrixBase< Rhs > &b) const
EIGEN_DEVICE_FUNC const SqrtReturnType sqrt() const
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
_OrderingType OrderingType
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Abs2ReturnType abs2() const
Matrix< Scalar, Dynamic, 1 > VectorSx
int EIGEN_BLAS_FUNC() swap(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Map< Matrix< Scalar, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > MatrixType
Index QuickSplit(VectorV &row, VectorI &ind, Index ncut)
SparseMatrix< Scalar, ColMajor, StorageIndex > FactorType
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
const VectorRx & scalingS() const
void setInitialShift(RealScalar shift)
Set the initial shift parameter .
OrderingType::PermutationType PermutationType
const AdjointReturnType adjoint() const
RealScalar m_initialShift
const PermutationType & permutationP() const
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T mini(const T &x, const T &y)
const StorageIndex * innerIndexPtr() const
A triangularView< Lower >().adjoint().solveInPlace(B)
A matrix or vector expression mapping an existing expression.
void analyzePattern(const MatrixType &mat)
Computes the fill reducing permutation vector using the sparsity pattern of mat.
Matrix< RealScalar, Dynamic, 1 > VectorRx
void compute(const MatrixType &mat)
NumTraits< Scalar >::Real RealScalar
#define eigen_internal_assert(x)
std::vector< std::list< StorageIndex > > VectorList
const Scalar * valuePtr() const
void _solve_impl(const Rhs &b, Dest &x) const