10 #ifndef EIGEN_INCOMPLETE_CHOlESKY_H 11 #define EIGEN_INCOMPLETE_CHOlESKY_H 13 #include <Eigen/OrderingMethods> 29 template <
typename Scalar,
int _UpLo = Lower,
typename _OrderingType = NaturalOrdering<
int> >
73 template<
typename MatrixType>
77 ord(mat.template selfadjointView<UpLo>(),
m_perm);
81 template<
typename MatrixType>
84 template<
typename MatrixType>
91 template<
typename Rhs,
typename Dest>
92 void _solve(
const Rhs& b, Dest& x)
const 99 x =
m_scal.asDiagonal() * x;
100 x =
m_L.template triangularView<UnitLower>().
solve(x);
104 x =
m_scal.asDiagonal() * x;
112 &&
"IncompleteLLT::solve(): invalid number of rows of the right hand side matrix b");
126 template <
typename IdxType,
typename SclType>
127 inline void updateList(
const IdxType& colPtr, IdxType& rowIdx, SclType& vals,
const Index&
col,
const Index& jk, IndexType& firstElt, VectorList& listCol);
130 template<
typename Scalar,
int _UpLo,
typename OrderingType>
131 template<
typename _MatrixType>
142 m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>().twistedBy(
m_perm);
144 m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
159 for (
int j = 0; j < n; j++)
165 Scalar mindiag = vals[0];
166 for (
int j = 0; j < n; j++){
167 for (
int k = colPtr[j]; k < colPtr[j+1]; k++)
169 mindiag = (min)(vals[colPtr[j]], mindiag);
174 for (
int j = 0; j < n; j++)
177 for (
int j=0; j < n; ++j)
181 Scalar diag = vals[colPtr[j]];
183 irow.setLinSpaced(n,0,n-1);
184 for (
int i = colPtr[j] + 1; i < colPtr[j+1]; i++)
186 curCol(rowIdx[i]) = vals[i];
187 irow(rowIdx[i]) = rowIdx[i];
189 std::list<int>::iterator k;
191 for(k = listCol[j].begin(); k != listCol[j].end(); k++)
193 int jk = firstElt(*k);
195 for (
int i = jk; i < colPtr[*k+1]; i++)
197 curCol(rowIdx[i]) -= vals[i] * vals[jk] ;
199 updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
205 std::cerr <<
"\nNegative diagonal during Incomplete factorization... "<< j <<
"\n";
210 vals[colPtr[j]] = rdiag;
211 for (
int i = j+1; i < n; i++)
216 vals[colPtr[i]] -= curCol(i) * curCol(i);
220 int p = colPtr[j+1] - colPtr[j] - 1 ;
224 for (
int i = colPtr[j]+1; i < colPtr[j+1]; i++)
226 vals[i] = curCol(cpt);
227 rowIdx[i] = irow(cpt);
231 Index jk = colPtr(j)+1;
232 updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
239 template<
typename Scalar,
int _UpLo,
typename OrderingType>
240 template <
typename IdxType,
typename SclType>
243 if (jk < colPtr(col+1) )
245 Index p = colPtr(col+1) - jk;
247 rowIdx.segment(jk,p).minCoeff(&minpos);
249 if (rowIdx(minpos) != rowIdx(jk))
252 std::swap(rowIdx(jk),rowIdx(minpos));
253 std::swap(vals(jk),vals(minpos));
256 listCol[rowIdx(jk)].push_back(col);
261 template<
typename _Scalar,
int _UpLo,
typename OrderingType,
typename Rhs>
268 template<typename Dest>
void evalTo(Dest& dst)
const 270 dec()._solve(
rhs(),dst);
void _solve(const Rhs &b, Dest &x) const
IntermediateState sqrt(const Expression &arg)
const AdjointReturnType adjoint() const
IncompleteCholesky< _Scalar, _UpLo, OrderingType > Dec
MatrixType::RealScalar RealScalar
SparseMatrix< Scalar, ColMajor > MatrixType
A matrix or vector expression mapping an existing array of data.
Modified Incomplete Cholesky with dual threshold.
void factorize(const MatrixType &amat)
iterative scaling algorithm to equilibrate rows and column norms in matrices
void setShift(Scalar shift)
Set the initial shift parameter.
Matrix< Scalar, Dynamic, 1 > ScalarType
const Index * outerIndexPtr() const
PermutationMatrix< Dynamic, Dynamic, Index > PermutationType
_OrderingType OrderingType
EIGEN_STRONG_INLINE void resize(Index nbRows, Index nbCols)
Derived & setZero(Index size)
Transpose< PermutationBase > inverse() const
Index QuickSplit(VectorV &row, VectorI &ind, Index ncut)
void rhs(const real_t *x, real_t *f)
void analyzePattern(const MatrixType &mat)
Computes the fill reducing permutation vector.
std::vector< std::list< Index > > VectorList
void compute(const MatrixType &matrix)
Matrix< Index, Dynamic, 1 > IndexType
IncompleteCholesky(const MatrixType &matrix)
const internal::solve_retval< IncompleteCholesky, Rhs > solve(const MatrixBase< Rhs > &b) const
const Index * innerIndexPtr() const
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
ComputationInfo info() const
Reports whether previous computation was successful.
NumTraits< Scalar >::Real RealScalar
void updateList(const IdxType &colPtr, IdxType &rowIdx, SclType &vals, const Index &col, const Index &jk, IndexType &firstElt, VectorList &listCol)
SparseMatrix< Scalar, ColMajor > m_L
const Scalar * valuePtr() const
Base class for all dense matrices, vectors, and expressions.
internal::traits< Derived >::Index Index