KSP

The KSP (Krylov Subspace Methods) module provides iterative linear solvers for systems of the form Ax = b. PETSc offers a wide variety of Krylov methods and preconditioners.

Overview

KSP provides:

Krylov methods: GMRES, CG, BiCGStab, and many more
Preconditioners: Jacobi, ILU, multigrid, direct solvers, etc.
Runtime configuration: Choose methods via command-line options
Convergence monitoring: Built-in residual tracking

Creating a KSP Solver

From a Matrix

# Basic creation with default options. The communicator comes from `A`
ksp = KSP(A)

# With the preconditioner construction matrix P
ksp = KSP(A, P)

# With options
ksp = KSP(A; 
    ksp_type = "gmres",
    pc_type = "ilu",
    ksp_rtol = 1e-8
)

From a DM

# Create KSP associated with a DM (for multigrid, etc.)
ksp = KSP(dm; 
    ksp_type = "cg",
    pc_type = "mg"
)

When the options apply

The keyword options of KSP, SNES and TS are applied through XSetFromOptions, and options_at says when:

`options_at`	Applied	Default for
`:construction`	in the constructor, once	`KSP`
`:each_solve`	at the start of every `solve!`	`SNES`, `TS`
`:first_solve`	once, at the first `setup!` or `solve!` (or `step!` for a TS)

The order matters because a setting made in code and an option for the same thing overwrite each other, and the later one wins. Under :each_solve, pc_type = "ilu" replaces a set_type!(pc(ksp), :jacobi) made before every solve. Under :construction, the options are read before anything set up after the constructor. :first_solve waits until the solver is complete and then reads them once:

ksp = KSP(dm; options_at = :first_solve, pc_type = "mg")
PETSc.set_dm_active!(ksp, false)   # the DM gives multigrid its hierarchy only
PETSc.set_operators!(ksp, J)
PETSc.solve!(x, ksp, b)            # options applied here, and not again

While a solver runs, its options stay pushed on PETSc's options stack, so sub-solvers PETSc creates during setup (the splits of a fieldsplit, the levels of a multigrid) read them too. PETSc.with_options does the same around any other code: with_options(ksp) do ... end.

From a Sparse Matrix

# Directly from Julia SparseMatrixCSC
using SparseArrays
S = sprand(100, 100, 0.1) + 10I
ksp = KSP(petsclib, MPI.COMM_SELF, S)

Solving

# Solve Ax = b, writing into x. The written vector comes first
PETSc.solve!(x, ksp, b)

# `ldiv!` and `\` are the LinearAlgebra spellings of the same call
using LinearAlgebra
ldiv!(x, ksp, b)
x = ksp \ b

# What PETSc did
PETSc.type_name(ksp)         # :gmres, a Symbol
PETSc.converged_reason(ksp)  # positive when it converged
PETSc.iteration_number(ksp)

PETSc.destroy!(ksp)

KSP is the type, not a factory function: ksp isa PETSc.KSP holds, and PETSc.dm(ksp) hands back the DM it was built on as a borrowed handle — destroy! on it is a no-op (naming conventions, §3.3).

Type names are Symbol at the Julia API: PETSc.set_type!(ksp, :cg) and PETSc.type_name(ksp) === :cg. The String spelling warns in v0.5 and is a MethodError in v0.6 (§3.1).

A nested solver

An inner solve, such as the Schur complement solve inside a preconditioner, starts from a bare KSP and an assembled operator. An options prefix keeps its options apart from the outer solver's:

inner = LibPETSc.KSPCreate(petsclib, comm)
PETSc.set_operators!(inner, S)              # S also builds the preconditioner; pass P to differ
PETSc.set_options_prefix!(inner, "schur_")  # reads -schur_ksp_type, -schur_pc_type, ...
PETSc.set_from_options!(inner)              # apply them now rather than at the first solve!
PETSc.solve!(y, inner, r)
PETSc.options_prefix(inner)                 # "schur_"

A KSP built on a DM computes its operator, right-hand side and initial guess from it. PETSc.set_dm_active!(ksp, false) keeps the DM for its geometry (multigrid needs it) and takes the rest from set_operators! and solve!; PETSc.set_dm_active!(ksp, :rhs, false) switches off one part.

Common Solver/Preconditioner Options

Krylov Methods (`ksp_type`)

cg - Conjugate Gradient (symmetric positive definite)
gmres - Generalized Minimum Residual
bicgstab - BiConjugate Gradient Stabilized
richardson - Richardson iteration
preonly - Apply preconditioner only (for direct solvers)

Preconditioners (`pc_type`)

jacobi - Diagonal scaling
ilu - Incomplete LU factorization
lu - Direct LU factorization
mg - Geometric multigrid
gamg - Algebraic multigrid
none - No preconditioning

Convergence Options

ksp_rtol - Relative tolerance (default: 1e-5)
ksp_atol - Absolute tolerance
ksp_max_it - Maximum iterations
ksp_monitor - Print residual each iteration

Example: Multigrid Solver

ksp = KSP(dm;
    ksp_type = "cg",
    pc_type = "mg",
    pc_mg_levels = 4,
    pc_mg_galerkin = true,
    mg_levels_ksp_type = "richardson",
    mg_levels_pc_type = "jacobi",
    mg_coarse_pc_type = "lu"
)

The preconditioner

PETSc.pc(ksp) hands back the preconditioner as a borrowed LibPETSc.PC. It takes the same set_type!/type_name pair as the solver, and is how a split preconditioner gets its index sets, which options alone cannot supply:

p = PETSc.pc(ksp)                         # not pc = pc(ksp), see naming.md §3.2
PETSc.set_type!(p, :fieldsplit)
PETSc.set_fieldsplit_is!(p, "u", is_u)    # rows of the first split, 0-based
PETSc.set_fieldsplit_is!(p, "p", is_p)    # options prefix -fieldsplit_p_
PETSc.type_name(p)                        # :fieldsplit

A :shell preconditioner is written in Julia. The action writes into its first argument, like every other callback here:

p = PETSc.pc(ksp)
PETSc.set_type!(p, :shell)
PETSc.set_shell_setup!(p) do p
    # rebuild whatever apply! needs, e.g. after the operator changed
end
PETSc.set_shell_apply!(p) do y, p, x
    PETSc.with_local_array!(y, x; read = (false, true), write = (true, false)) do ya, xa
        ya .= xa ./ diagonal                 # Jacobi, by hand
    end
end

The closures are kept with the PETSc preconditioner, not with the wrapper p, so p can be dropped. An exception thrown inside either one comes out of the solve that ran it, as itself: a DomainError in apply! makes ksp \ b throw that DomainError (naming conventions, §18.4).

The solvers of a :fieldsplit's splits exist only once the preconditioner has been set up, inside the first solve. set_post_setup! runs a closure at the end of every setup, which is where to configure them:

PETSc.set_post_setup!(p) do p
    _, subksps = LibPETSc.PCFieldSplitGetSubKSP(petsclib, p)
    PETSc.set_type!(PETSc.pc(subksps[1]), :lu)
end

It runs again whenever the operator changes, and the options given to a SNES are applied again on every solve!, so write it so that a second run is harmless. With a Schur split, set the factorization type in code (LibPETSc.PCFieldSplitSetType) before LibPETSc.PCFieldSplitSetSchurPre with a user matrix: when the type comes only from the options, PETSc drops the user matrix without saying so. set_pre_solve! is the matching hook on the KSP; it runs at the start of every solve, before the preconditioner is set up.

A LibPETSc.PC is also what every low-level PC* function takes, so anything without a high-level verb is one call away: LibPETSc.PCFieldSplitSetType(petsclib, p, LibPETSc.PC_COMPOSITE_SCHUR).

Functions

PETSc.LibPETSc.KSP — Method

KSP(petsclib, comm::MPI.Comm, A::SparseMatrixCSC; options...)

Create a KSP for the sparse matrix A, on comm, using the library petsclib. A is copied into a new PETSc matrix that the solver owns, so changing A afterwards does not change the solver. The keyword arguments are PETSc options, as for KSP(A::AbstractPetscMat; options...).