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

Sequential IC(0) on CPU

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(...)
              .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<gko::solver::Bicgstab<double>>::build()
    .with_factorization(
        gko::factorization::ParIc<double, gko::int32>::build().on(exec))
    .with_l_solver(
        gko::solver::Bicgstab<double>::build()
            .with_criteria(
                gko::stop::Iteration::build()
                    .with_max_iters(3u).on(exec))
            .on(exec))
    .on(exec);

The first template parameter of preconditioner::Ic is the triangular solver type; the default is solver::LowerTrs. The \(L^H\) solve is constructed automatically as the conjugate-transpose of the \(L\) solve, so only one solver factory is configured.

Factory parameters#

Parameter

Type

Purpose

factorization

LinOpFactory

The factorisation that produces the \(L\) (and implicitly \(L^H\)) factor. Typically Ic, ParIc, or ParIct. Required.

l_solver

LinOpFactory

Factory for the \(L\) triangular solver. Defaults to a direct triangular solve.

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: