Go to the documentation of this file.
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>
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>
156 x =
m_L.adjoint().template triangularView<Upper>().solve(
x);
188 template<
typename Scalar,
int _UpLo,
typename OrderingType>
189 template<
typename _MatrixType>
193 eigen_assert(m_analysisIsOk &&
"analyzePattern() should be called first");
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>();
211 Index nnz = m_L.nonZeros();
220 col_pattern.fill(-1);
228 for (
Index k = colPtr[
j]; k < colPtr[
j+1]; k++)
235 m_scale = m_scale.cwiseSqrt().cwiseSqrt();
249 for (
Index k = colPtr[
j]; k < colPtr[
j+1]; k++)
250 vals[k] *= (m_scale(
j)*m_scale(rowIdx[k]));
251 eigen_internal_assert(rowIdx[colPtr[
j]]==
j &&
"IncompleteCholesky: only the lower triangular part must be stored");
259 shift = m_initialShift - mindiag;
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;
361 updateList(colPtr,rowIdx,vals,
j,jk,firstElt,listCol);
366 m_factorizationIsOk =
true;
372 template<
typename Scalar,
int _UpLo,
typename OrderingType>
375 if (jk < colPtr(
col+1) )
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 resize(Index rows, Index cols)
Namespace containing all symbols from the Eigen library.
RealScalar m_initialShift
Array< double, 1, 3 > e(1./3., 0.5, 2.)
SparseMatrix< Scalar, ColMajor, StorageIndex > FactorType
void updateList(Ref< const VectorIx > colPtr, Ref< VectorIx > rowIdx, Ref< VectorSx > vals, const Index &col, const Index &jk, VectorIx &firstElt, VectorList &listCol)
void _solve_impl(const Rhs &b, Dest &x) const
void compute(const MatrixType &mat)
Matrix diag(const std::vector< Matrix > &Hs)
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
IncompleteCholesky(const MatrixType &matrix)
Modified Incomplete Cholesky with dual threshold.
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
#define eigen_internal_assert(x)
const Scalar * valuePtr() const
SparseSolverBase< IncompleteCholesky< Scalar, _UpLo, _OrderingType > > Base
OrderingType::PermutationType PermutationType
const StorageIndex * innerIndexPtr() const
std::vector< std::list< StorageIndex > > VectorList
void factorize(const MatrixType &mat)
Performs the numerical factorization of the input matrix mat.
Matrix< StorageIndex, Dynamic, 1 > VectorIx
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T mini(const T &x, const T &y)
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
static const Line3 l(Rot3(), 1, 1)
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
_OrderingType OrderingType
const VectorRx & scalingS() const
const PermutationType & permutationP() const
AnnoyingScalar conj(const AnnoyingScalar &x)
A matrix or vector expression mapping an existing array of data.
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
const StorageIndex * outerIndexPtr() const
NumTraits< Scalar >::Real RealScalar
NumTraits< Scalar >::Real RealScalar
const FactorType & matrixL() const
A base class for sparse solvers.
EIGEN_DEVICE_FUNC bool abs2(bool x)
A matrix or vector expression mapping an existing expression.
Matrix< RealScalar, Dynamic, 1 > VectorRx
PermutationType::StorageIndex StorageIndex
iterator iter(handle obj)
Index QuickSplit(VectorV &row, VectorI &ind, Index ncut)
ComputationInfo info() const
Reports whether previous computation was successful.
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Matrix< Scalar, Dynamic, 1 > VectorSx
A triangularView< Lower >().adjoint().solveInPlace(B)
void setInitialShift(RealScalar shift)
Set the initial shift parameter .
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Jet< T, N > sqrt(const Jet< T, N > &f)
void analyzePattern(const MatrixType &mat)
Computes the fill reducing permutation vector using the sparsity pattern of mat.
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:45