Linear solvers (LS)

• Linear solvers (LS)
• • List of implemented linear solvers
  • Preconditioner

If you provide your own linear solver, it must be a subtype of AbstractLinearSolver otherwise BifurcationKit.jl will not recognize it. See example just below.

The linear solvers provide a way of inverting the Jacobian J, i.e. of solving $J\, x = \mathrm{rhs}$. Throughout the package, such a linear solver linsolve will be called like

sol, success, itnumber = linsolve(J, rhs; kwargs...)

The three returned quantities are the solution sol, a boolean success telling whether the solve succeeded (resp. converged for iterative solvers) and itnumber the number of iterations required by the computation (equal to 1 for a direct solver).

Here is an example of the simplest one (see src/LinearSolver.jl for the true implementation) to give you an idea, the backslash operator:

struct DefaultLS <: AbstractLinearSolver endfunction (l::DefaultLS)(J, rhs; k...)	return J \ rhs, true, 1end

You can then call it as follows (and it will be called like this in newton):

ls = DefaultLS()J = rand(2, 2) # example of linear operatorls(J, rand(2))

List of implemented linear solvers

The detailed information for each one of them is located in the API.

1. DefaultLS — direct solver based on the default \ operator, which dispatches to LU, Cholesky or a sparse factorization depending on the type of the Jacobian. It works for dense and sparse matrices. The option useFactorization = true (the default) caches the factorization so that a second linear solve with the same operator is very cheap; set it to false when your operator does not support factorize (e.g. some operators coming from ApproxFun.jl, or matrices such as StaticArrays.MMatrix). You can create one via linsolver = DefaultLS().

2. GMRESIterativeSolvers — iterative solver based on gmres from IterativeSolvers.jl. You can create one via linsolver = GMRESIterativeSolvers() and tune it with the keyword arguments abstol, reltol, restart, maxiter, N, verbose, ... It supports Matrix-Free Jacobians (with N the dimension of the state space) and left/right preconditioners (fields Pl, Pr). The struct is mutable so that you can change the preconditioners on the fly.

3. GMRESKrylovKit — iterative solver based on linsolve from KrylovKit.jl. You can create one via linsolver = GMRESKrylovKit() and tune it with the keyword arguments dim, atol, rtol, maxiter, verbose, ...

   !!! tip "Different linear solvers" By tuning the options of GMRESKrylovKit, you can select CG, GMRES... Indeed, KrylovKit.jl dispatches on the options issymmetric, ishermitian, isposdef, see KrylovKit.jl. A left preconditioner can be set through the field Pl.

4. KrylovLS / KrylovLSInplace — interfaces to the many solvers of Krylov.jl. You can create one via linsolver = KrylovLS() and pass the required Krylov method through the keyword KrylovAlg = :gmres (the default): you have access to cg, cr, gmres, symmlq, minres, cg_lanczos, cg_lanczos_shift_seq, ... The extra keyword arguments are forwarded to the chosen Krylov method, e.g. KrylovLS(atol = 1e-11, rtol = 1e-8). Left/right preconditioners are passed through the keywords Pl, Pr (see the Preconditioner section below). KrylovLSInplace pre-allocates the Krylov space (you then have to provide the workspace dimensions, see the docstring) which is very useful on the GPU or when many linear solves are performed. Both structs are mutable so that you can modify the preconditioners.

   !!! note "Other solvers" Thanks to the interface to Krylov.jl, you do not need to implement the Conjugate Gradients from IterativeSolvers.jl, minres, ... yourself: just select the algorithm with the keyword KrylovAlg as explained above.

Preconditioner

Preconditioners should be considered when using Matrix Free methods such as GMRES. GMRESIterativeSolvers provides a very simple interface for using them (fields Pl, Pr). For GMRESKrylovKit, we implemented a left preconditioner (field Pl). Note that, for GMRESKrylovKit, you are not restricted to use Vectors anymore.

A curated list is provided by the package LinearSolve.jl.

Finally, here are some packages to use preconditioners:

1. Krylov.jl provides some preconditioners.
2. IncompleteLU.jl an ILU like preconditioner
3. AlgebraicMultigrid.jl Algebraic Multigrid (AMG) preconditioners. This works especially well for symmetric positive definite matrices.
4. Preconditioners.jl A convenient interface to conveniently called most of the above preconditioners using a single syntax.
5. We provide a preconditioner based on deflation of eigenvalues (also called preconditioner based on Leading Invariant Subspaces) using a partial Schur decomposition. There are two ways to define one i.e. PrecPartialSchurKrylovKit and PrecPartialSchurArnoldiMethod.