blob: e697f32530cbedb125221f37bee90963b3775e7a [file] [log] [blame]
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
// Copyright (C) 2012 Désiré Nuentsa-Wakam <>
// Copyright (C) 2015 Gael Guennebaud <>
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at
#include <vector>
#include <list>
#include "./InternalHeaderCheck.h"
namespace Eigen {
* \brief Modified Incomplete Cholesky with dual threshold
* References : C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with
* Limited memory, SIAM J. Sci. Comput. 21(1), pp. 24-45, 1999
* \tparam Scalar the scalar type of the input matrices
* \tparam UpLo_ The triangular part that will be used for the computations. It can be Lower
* or Upper. Default is Lower.
* \tparam OrderingType_ The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<int>,
* unless EIGEN_MPL2_ONLY is defined, in which case the default is NaturalOrdering<int>.
* \implsparsesolverconcept
* It performs the following incomplete factorization: \f$ S P A P' S \approx L L' \f$
* where L is a lower triangular factor, S is a diagonal scaling matrix, and P is a
* fill-in reducing permutation as computed by the ordering method.
* \b Shifting \b strategy: Let \f$ B = S P A P' S \f$ be the scaled matrix on which the factorization is carried out,
* and \f$ \beta \f$ be the minimum value of the diagonal. If \f$ \beta > 0 \f$ then, the factorization is directly performed
* on the matrix B. Otherwise, the factorization is performed on the shifted matrix \f$ B + (\sigma+|\beta| I \f$ where
* \f$ \sigma \f$ is the initial shift value as returned and set by setInitialShift() method. The default value is \f$ \sigma = 10^{-3} \f$.
* If the factorization fails, then the shift in doubled until it succeed or a maximum of ten attempts. If it still fails, as returned by
* the info() method, then you can either increase the initial shift, or better use another preconditioning technique.
template <typename Scalar, int UpLo_ = Lower, typename OrderingType_ = AMDOrdering<int> >
class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,UpLo_,OrderingType_> >
typedef SparseSolverBase<IncompleteCholesky<Scalar,UpLo_,OrderingType_> > Base;
using Base::m_isInitialized;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef OrderingType_ OrderingType;
typedef typename OrderingType::PermutationType PermutationType;
typedef typename PermutationType::StorageIndex StorageIndex;
typedef SparseMatrix<Scalar,ColMajor,StorageIndex> FactorType;
typedef Matrix<Scalar,Dynamic,1> VectorSx;
typedef Matrix<RealScalar,Dynamic,1> VectorRx;
typedef Matrix<StorageIndex,Dynamic, 1> VectorIx;
typedef std::vector<std::list<StorageIndex> > VectorList;
enum { UpLo = UpLo_ };
enum {
ColsAtCompileTime = Dynamic,
MaxColsAtCompileTime = Dynamic
/** Default constructor leaving the object in a partly non-initialized stage.
* You must call compute() or the pair analyzePattern()/factorize() to make it valid.
* \sa IncompleteCholesky(const MatrixType&)
IncompleteCholesky() : m_initialShift(1e-3),m_analysisIsOk(false),m_factorizationIsOk(false) {}
/** Constructor computing the incomplete factorization for the given matrix \a matrix.
template<typename MatrixType>
IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_analysisIsOk(false),m_factorizationIsOk(false)
/** \returns number of rows of the factored matrix */
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_L.rows(); }
/** \returns number of columns of the factored matrix */
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_L.cols(); }
/** \brief Reports whether previous computation was successful.
* It triggers an assertion if \c *this has not been initialized through the respective constructor,
* or a call to compute() or analyzePattern().
* \returns \c Success if computation was successful,
* \c NumericalIssue if the matrix appears to be negative.
ComputationInfo info() const
eigen_assert(m_isInitialized && "IncompleteCholesky is not initialized.");
return m_info;
/** \brief Set the initial shift parameter \f$ \sigma \f$.
void setInitialShift(RealScalar shift) { m_initialShift = shift; }
/** \brief Computes the fill reducing permutation vector using the sparsity pattern of \a mat
template<typename MatrixType>
void analyzePattern(const MatrixType& mat)
OrderingType ord;
PermutationType pinv;
ord(mat.template selfadjointView<UpLo>(), pinv);
if(pinv.size()>0) m_perm = pinv.inverse();
else m_perm.resize(0);
m_L.resize(mat.rows(), mat.cols());
m_analysisIsOk = true;
m_isInitialized = true;
m_info = Success;
/** \brief Performs the numerical factorization of the input matrix \a mat
* The method analyzePattern() or compute() must have been called beforehand
* with a matrix having the same pattern.
* \sa compute(), analyzePattern()
template<typename MatrixType>
void factorize(const MatrixType& mat);
/** Computes or re-computes the incomplete Cholesky factorization of the input matrix \a mat
* It is a shortcut for a sequential call to the analyzePattern() and factorize() methods.
* \sa analyzePattern(), factorize()
template<typename MatrixType>
void compute(const MatrixType& mat)
// internal
template<typename Rhs, typename Dest>
void _solve_impl(const Rhs& b, Dest& x) const
eigen_assert(m_factorizationIsOk && "factorize() should be called first");
if (m_perm.rows() == b.rows()) x = m_perm * b;
else x = b;
x = m_scale.asDiagonal() * x;
x = m_L.template triangularView<Lower>().solve(x);
x = m_L.adjoint().template triangularView<Upper>().solve(x);
x = m_scale.asDiagonal() * x;
if (m_perm.rows() == b.rows())
x = m_perm.inverse() * x;
/** \returns the sparse lower triangular factor L */
const FactorType& matrixL() const { eigen_assert(m_factorizationIsOk && "factorize() should be called first"); return m_L; }
/** \returns a vector representing the scaling factor S */
const VectorRx& scalingS() const { eigen_assert(m_factorizationIsOk && "factorize() should be called first"); return m_scale; }
/** \returns the fill-in reducing permutation P (can be empty for a natural ordering) */
const PermutationType& permutationP() const { eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); return m_perm; }
FactorType m_L; // The lower part stored in CSC
VectorRx m_scale; // The vector for scaling the matrix
RealScalar m_initialShift; // The initial shift parameter
bool m_analysisIsOk;
bool m_factorizationIsOk;
ComputationInfo m_info;
PermutationType m_perm;
inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol);
// Based on the following paper:
// C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with
// Limited memory, SIAM J. Sci. Comput. 21(1), pp. 24-45, 1999
template<typename Scalar, int UpLo_, typename OrderingType>
template<typename MatrixType_>
void IncompleteCholesky<Scalar,UpLo_, OrderingType>::factorize(const MatrixType_& mat)
using std::sqrt;
eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
// Dropping strategy : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added
// Apply the fill-reducing permutation computed in analyzePattern()
if (m_perm.rows() == mat.rows() ) // To detect the null permutation
// The temporary is needed to make sure that the diagonal entry is properly sorted
FactorType tmp(mat.rows(), mat.cols());
tmp = mat.template selfadjointView<UpLo_>().twistedBy(m_perm);
m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
m_L.template selfadjointView<Lower>() = mat.template selfadjointView<UpLo_>();
Index n = m_L.cols();
Index nnz = m_L.nonZeros();
Map<VectorSx> vals(m_L.valuePtr(), nnz); //values
Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz); //Row indices
Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row
VectorIx firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
VectorList listCol(n); // listCol(j) is a linked list of columns to update column j
VectorSx col_vals(n); // Store a nonzero values in each column
VectorIx col_irow(n); // Row indices of nonzero elements in each column
VectorIx col_pattern(n);
StorageIndex col_nnz;
// Computes the scaling factors
for (Index j = 0; j < n; j++)
for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
m_scale(j) += numext::abs2(vals(k));
m_scale(rowIdx[k]) += numext::abs2(vals(k));
m_scale = m_scale.cwiseSqrt().cwiseSqrt();
for (Index j = 0; j < n; ++j)
m_scale(j) = RealScalar(1)/m_scale(j);
m_scale(j) = 1;
// TODO disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster)
// Scale and compute the shift for the matrix
RealScalar mindiag = NumTraits<RealScalar>::highest();
for (Index j = 0; j < n; j++)
for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
vals[k] *= (m_scale(j)*m_scale(rowIdx[k]));
eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored");
mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
FactorType L_save = m_L;
RealScalar shift = 0;
if(mindiag <= RealScalar(0.))
shift = m_initialShift - mindiag;
m_info = NumericalIssue;
// Try to perform the incomplete factorization using the current shift
int iter = 0;
// Apply the shift to the diagonal elements of the matrix
for (Index j = 0; j < n; j++)
vals[colPtr[j]] += shift;
// jki version of the Cholesky factorization
Index j=0;
for (; j < n; ++j)
// Left-looking factorization of the j-th column
// First, load the j-th column into col_vals
Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored
col_nnz = 0;
for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
StorageIndex l = rowIdx[i];
col_vals(col_nnz) = vals[i];
col_irow(col_nnz) = l;
col_pattern(l) = col_nnz;
typename std::list<StorageIndex>::iterator k;
// Browse all previous columns that will update column j
for(k = listCol[j].begin(); k != listCol[j].end(); k++)
Index jk = firstElt(*k); // First element to use in the column
Scalar v_j_jk = numext::conj(vals[jk]);
jk += 1;
for (Index i = jk; i < colPtr[*k+1]; i++)
StorageIndex l = rowIdx[i];
col_vals(col_nnz) = vals[i] * v_j_jk;
col_irow[col_nnz] = l;
col_pattern(l) = col_nnz;
col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
// Scale the current column
if(numext::real(diag) <= 0)
// increase shift
shift = numext::maxi(m_initialShift,RealScalar(2)*shift);
// restore m_L, col_pattern, and listCol
vals = Map<const VectorSx>(L_save.valuePtr(), nnz);
rowIdx = Map<const VectorIx>(L_save.innerIndexPtr(), nnz);
colPtr = Map<const VectorIx>(L_save.outerIndexPtr(), n+1);
for(Index i=0; i<n; ++i)
RealScalar rdiag = sqrt(numext::real(diag));
vals[colPtr[j]] = rdiag;
for (Index k = 0; k<col_nnz; ++k)
Index i = col_irow[k];
col_vals(k) /= rdiag;
//Update the remaining diagonals with col_vals
vals[colPtr[i]] -= numext::abs2(col_vals(k));
// Select the largest p elements
// p is the original number of elements in the column (without the diagonal)
Index p = colPtr[j+1] - colPtr[j] - 1 ;
Ref<VectorSx> cvals = col_vals.head(col_nnz);
Ref<VectorIx> cirow = col_irow.head(col_nnz);
internal::QuickSplit(cvals,cirow, p);
// Insert the largest p elements in the matrix
Index cpt = 0;
for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
vals[i] = col_vals(cpt);
rowIdx[i] = col_irow(cpt);
// restore col_pattern:
col_pattern(col_irow(cpt)) = -1;
// Get the first smallest row index and put it after the diagonal element
Index jk = colPtr(j)+1;
m_factorizationIsOk = true;
m_info = Success;
} while(m_info!=Success);
template<typename Scalar, int UpLo_, typename OrderingType>
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)
if (jk < colPtr(col+1) )
Index p = colPtr(col+1) - jk;
Index minpos;
minpos += jk;
if (rowIdx(minpos) != rowIdx(jk))
firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
} // end namespace Eigen