Linear Solve

AbstractLinearSolverCache

Abstract Type for all Linear Solvers used in NonlinearSolveBase. Subtypes of these are meant to be constructured via construct_linear_solver.

construct_linear_solver(alg, linsolve, A, b, u, p; stats, kwargs...)

Construct a cache for solving linear systems of the form A * u = b. Following cases are handled:

1. A is Number, then we solve it with u = b / A
2. A is SMatrix, then we solve it with u = A \ b (using the defaults from base Julia) (unless a preconditioner is specified)
3. If linsolve is \, then we solve it with directly using ldiv!(u, A, b)
4. In all other cases, we use alg to solve the linear system using LinearSolve.jl

Solving the System

(cache::LinearSolverCache)(;
    A = nothing, b = nothing, linu = nothing, reuse_A_if_factorization = false, kwargs...
)

Returns the solution of the system u and stores the updated cache in cache.lincache.

Special Handling for Rank-deficient Matrix A

If we detect a failure in the linear solve (mostly due to using an algorithm that doesn't support rank-deficient matrices), we emit a warning and attempt to solve the problem using Pivoted QR factorization. This is quite efficient if there are only a few rank-deficient that originate in the problem. However, if these are quite frequent for the main nonlinear system, then it is recommended to use a different linear solver that supports rank-deficient matrices.

Keyword Arguments

• reuse_A_if_factorization: If true, then the factorization of A is reused if possible. This is useful when solving the same system with different b values. If the algorithm is an iterative solver, then we reset the internal linear solve cache.

• alias: A LinearAliasSpecifier forwarded to the LinearSolve cache construction, or nothing (the default) to let the alias_A_for_refactorization trait decide: unalias both A and b so the caller's arrays are used as caches safely, unless the LinearSolve extension opts the algorithm into an owned-copy refactorization buffer. Callers that own the passed A/b outright can pass a LinearAliasSpecifier explicitly so factorizations may work in place without any defensive copy.

One distinct feature of this compared to the cache from LinearSolve is that it respects the aliasing arguments even after cache construction, i.e., if we passed in an A that A is not mutated, we do this by copying over A to a preconstructed cache.

get_linear_cache(cache) -> Union{Nothing, LinearSolve.LinearCache}

Return the LinearSolve.jlLinearCache the solver holds its Jacobian (and, for a factorization, that factorization) in, or nothing when there is no single reusable linear cache.

This is a read-only accessor: it hands back the cache the solver already maintains so a caller can reuse the current Jacobian for an auxiliary linear solve J x = b (for example a smoothed error estimate) instead of building and factorizing a second copy. It does not refactorize; the returned cache is in exactly the state the last step!/solve! left it, so call it after the solve for the current point. Reuse it through the normal LinearSolve caching interface (set b/u and solve!, passing no new A), which reuses the existing factorization for a direct solver and re-runs the iterative solve for a Krylov cache.

Returns nothing for algorithms with no single descent linear solve (e.g. polyalgorithms) and for the native \/SMatrix/Number/Diagonal paths that hold no reusable LinearCache; callers should fall back accordingly.

needs_square_A(linsolve, u)::Bool

Return whether linsolve requires a square matrix for a state shaped like u.

NonlinearSolveBase uses this developer trait to decide whether least-squares nonlinear problems should build normal equations before calling a linear solver. LinearSolve extensions add methods for concrete LinearSolve algorithms.

Arguments

• linsolve: A linear solver algorithm, nothing, or \.
• u: Current nonlinear state, used to distinguish scalar and array cases.

Examples

using NonlinearSolveBase

NonlinearSolveBase.needs_square_A(nothing, [1.0, 2.0])

needs_concrete_A(linsolve)::Bool

Return whether linsolve requires NonlinearSolveBase to materialize a concrete Jacobian.

Matrix-free Jacobian operators can be used only when this trait is false.

Arguments

• linsolve: A linear solver algorithm, missing, nothing, or \.

Examples

using NonlinearSolveBase

NonlinearSolveBase.needs_concrete_A(\)