IC family (IC, ParIC, ParICT)#

The Incomplete Cholesky family preconditions a symmetric positive-definite system \(A\) by computing an approximation \(L L^H \approx A\) and using its action as the preconditioner: at apply time, the preconditioner solves \(L L^H z = r\) by chaining a forward triangular solve with its conjugate-transpose backward solve.

Three pieces work together:

The matching preconditioner LinOp, gko::preconditioner::Ic, wraps the factorisation result and configures the triangular-solve LinOp used at apply time. The default solver is a direct LowerTrs; a Krylov inner solver can be substituted for cases where the triangular solve is the bottleneck.

When to use IC#

  • SPD matrices where Jacobi is too weak but a full factorisation would not fit in memory. IC’s storage is bounded by the sparsity pattern of \(A\) (for IC(0)) or by a tunable budget (for ParICT).

  • As the preconditioner inside CG, FCG, or Minres: IC preserves symmetry, so the preconditioned operator stays SPD.

  • Pair with a fill-reducing reordering (AMD, nested dissection) for higher quality on mesh-based problems. The reordering improves the sparsity pattern that IC(0) is constrained to.

For non-symmetric matrices, use the ILU family instead.

Picking a variant#

Variant

Algorithm

Best when

factorization::Ic

Exact IC(0) on the input sparsity pattern: cuSPARSE/hipSPARSE or Ginkgo’s sync-free kernel, chosen by the algorithm parameter; runs on the Reference, OpenMP, CUDA and HIP executors

Small to moderate matrices; the canonical correctness baseline.

factorization::ParIc

Asynchronous fixed-point IC(0)

Larger matrices on GPU or many-core CPU; same sparsity pattern as Ic.

factorization::ParIct

Threshold-based IC with sparsity-pattern refinement

When IC(0) is too sparse to be effective; ParICT improves the pattern at the cost of more storage and longer generation.

Construction#

// Standard IC(0) preconditioner inside CG
auto cg = gko::solver::Cg<double>::build()
              .with_criteria(gko::stop::max_iters(1000u),
                             gko::stop::relative_residual_norm(1e-8))
              .with_preconditioner(
                  gko::preconditioner::Ic<>::build()
                      .with_factorization(
                          gko::factorization::ParIc<double, gko::int32>::build()
                              .on(exec))
                      .on(exec))
              .on(exec);

// With an iterative inner solve for the triangular system
auto preconditioner = gko::preconditioner::Ic<double, gko::int32>::build()
    .with_factorization(
        gko::factorization::ParIc<double, gko::int32>::build())
    .with_l_solver(
        gko::solver::Bicgstab<double>::build().with_criteria(
            gko::stop::max_iters(3u)))
    .on(exec);

preconditioner::Ic<ValueType, IndexType> defaults to <default_precision, int32>. The \(L\) solver is chosen at runtime with with_l_solver, which accepts any LinOpFactory whose generated solver is Transposable; when unset, Ic generates solver::LowerTrs<ValueType, IndexType>. The \(L^H\) solve is the conjugate transpose of the generated \(L\) solver, so only one solver factory is configured.

Factory parameters#

Parameter

Type

Purpose

factorization

LinOpFactory

Factorisation that produces \(L\): Ic, ParIc or ParIct. Defaults to ParIc computing only \(L\). Not used when generate() receives a Composition holding the factor(s) directly.

l_solver

LinOpFactory

Factory for the \(L\) triangular solver. Defaults to solver::LowerTrs.

Implementation notes#

The triangular-solve dispatcher caches the lower factor and its conjugate-transpose adapter, so each apply runs two triangular solves without rebuilding the operator chain. ParIC and ParICT generation are asynchronous fixed-point iterations of Chow & Patel; the number of sweeps needed for a usable preconditioner depends on the parallelism level and the problem (a single sweep is enough for sequential execution; 3 sweeps on OpenMP and 10+ on GPU are typical).

See also

References: