11 #ifndef EIGEN_INCOMPLETE_CHOlESKY_H 12 #define EIGEN_INCOMPLETE_CHOlESKY_H 44 template <
typename Scalar,
int _UpLo = Lower,
typename _OrderingType = AMDOrdering<
int> >
59 typedef std::vector<std::list<StorageIndex> >
VectorList;
77 template<
typename MatrixType>
110 template<
typename MatrixType>
114 PermutationType pinv;
115 ord(mat.template selfadjointView<UpLo>(), pinv);
116 if(pinv.size()>0)
m_perm = pinv.inverse();
131 template<
typename MatrixType>
140 template<
typename MatrixType>
148 template<
typename Rhs,
typename Dest>
158 if (
m_perm.rows() == b.rows())
188 template<
typename Scalar,
int _UpLo,
typename OrderingType>
189 template<
typename _MatrixType>
198 if (
m_perm.rows() == mat.rows() )
202 tmp = mat.template selfadjointView<_UpLo>().twistedBy(
m_perm);
203 m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
207 m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
220 col_pattern.fill(-1);
228 for (
Index k = colPtr[
j]; k < colPtr[
j+1]; k++)
249 for (
Index k = colPtr[
j]; k < colPtr[
j+1]; k++)
251 eigen_internal_assert(rowIdx[colPtr[
j]]==
j &&
"IncompleteCholesky: only the lower triangular part must be stored");
269 vals[colPtr[
j]] += shift;
279 for (
Index i = colPtr[j] + 1;
i < colPtr[j+1];
i++)
282 col_vals(col_nnz) = vals[
i];
283 col_irow(col_nnz) =
l;
284 col_pattern(l) = col_nnz;
288 typename std::list<StorageIndex>::iterator k;
290 for(k = listCol[j].begin(); k != listCol[
j].end(); k++)
292 Index jk = firstElt(*k);
297 for (
Index i = jk;
i < colPtr[*k+1];
i++)
302 col_vals(col_nnz) = vals[
i] * v_j_jk;
303 col_irow[col_nnz] =
l;
304 col_pattern(l) = col_nnz;
308 col_vals(col_pattern[l]) -= vals[
i] * v_j_jk;
310 updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
326 col_pattern.fill(-1);
334 vals[colPtr[
j]] = rdiag;
335 for (
Index k = 0; k<col_nnz; ++k)
339 col_vals(k) /= rdiag;
345 Index p = colPtr[j+1] - colPtr[
j] - 1 ;
351 for (
Index i = colPtr[j]+1;
i < colPtr[j+1];
i++)
353 vals[
i] = col_vals(cpt);
354 rowIdx[
i] = col_irow(cpt);
356 col_pattern(col_irow(cpt)) = -1;
360 Index jk = colPtr(j)+1;
361 updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
372 template<
typename Scalar,
int _UpLo,
typename OrderingType>
375 if (jk < colPtr(col+1) )
377 Index p = colPtr(col+1) - jk;
379 rowIdx.segment(jk,p).minCoeff(&minpos);
381 if (rowIdx(minpos) != rowIdx(jk))
387 firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
388 listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
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)
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
Matrix diag(const std::vector< Matrix > &Hs)
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)
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
Namespace containing all symbols from the Eigen library.
iterator iter(handle obj)
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)
AnnoyingScalar conj(const AnnoyingScalar &x)
_OrderingType OrderingType
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Matrix< Scalar, Dynamic, 1 > VectorSx
static const Line3 l(Rot3(), 1, 1)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Index QuickSplit(VectorV &row, VectorI &ind, Index ncut)
SparseMatrix< Scalar, ColMajor, StorageIndex > FactorType
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Array< double, 1, 3 > e(1./3., 0.5, 2.)
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)
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
const StorageIndex * innerIndexPtr() const
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.
A triangularView< Lower >().adjoint().solveInPlace(B)
Matrix< RealScalar, Dynamic, 1 > VectorRx
Jet< T, N > sqrt(const Jet< T, N > &f)
void compute(const MatrixType &mat)
NumTraits< Scalar >::Real RealScalar
#define eigen_internal_assert(x)
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
EIGEN_DEVICE_FUNC bool abs2(bool x)
std::vector< std::list< StorageIndex > > VectorList
const Scalar * valuePtr() const
void _solve_impl(const Rhs &b, Dest &x) const