Merge pull request #91 from SciML/docs

Document developing algorithms and keyword arguments

Author: ChrisRackauckas

docs/make.jl

@@ -19,12 +19,16 @@ makedocs(
         ],
         "Basics" => Any[
             "basics/LinearProblem.md",
+            "basics/common_solver_opts.md",
             "basics/CachingAPI.md",
             "basics/Preconditioners.md",
             "basics/FAQ.md"
         ],
         "Solvers" => Any[
             "solvers/solvers.md"
+        ],
+        "Advanced" => Any[
+            "advanced/developing.md"
         ]
     ]
 )

docs/src/advanced/developing.md

@@ -0,0 +1,60 @@
+# Developing New Linear Solvers
+
+Developing new or custom linear solvers for the SciML interface can be done in
+one of two ways:
+
+1. You can either create a completely new set of dispatches for `init` and `solve`.
+2. You can extend LinearSolve.jl's internal mechanisms.
+
+For developer ease, we highly recommend (2) as that will automatically make the
+caching API work. Thus this is the documentation for how to do that.
+
+## Developing New Linear Solvers with LinearSolve.jl Primitives
+
+Let's create a new wrapper for a simple LU-factorization which uses only the
+basic machinery. A simplified version is:
+
+```julia
+struct MyLUFactorization{P} <: SciMLBase.AbstractLinearAlgorithm end
+
+init_cacheval(alg::MyLUFactorization, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose) = lu!(convert(AbstractMatrix,A))
+
+function SciMLBase.solve(cache::LinearCache, alg::MyLUFactorization; kwargs...)
+    if cache.isfresh
+        A = convert(AbstractMatrix,A)
+        fact = lu!(A)
+        cache = set_cacheval(cache, fact)
+    end
+    y = ldiv!(cache.u, cache.cacheval, cache.b)
+    SciMLBase.build_linear_solution(alg,y,nothing,cache)
+end
+```
+
+The way this works is as follows. LinearSolve.jl has a `LinearCache` that everything
+shares (this is what gives most of the ease of use). However, many algorithms
+need to cache their own things, and so there's one value `cacheval` that is
+for the algorithms to modify. The function:
+
+```julia
+init_cacheval(alg::MyLUFactorization, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
+```
+
+is what is called at `init` time to create the first `cacheval`. Note that this
+should match the type of the cache later used in `solve` as many algorithms, like
+those in OrdinaryDiffEq.jl, expect type-groundedness in the linear solver definitions.
+While there are cheaper ways to obtain this type for LU factorizations (specifically,
+`ArrayInterface.lu_instance(A)`), for a demonstration this just performs an
+LU-factorization to get an `LU{T, Matrix{T}}` which it puts into the `cacheval`
+so its typed for future use.
+
+After the `init_cacheval`, the only thing left to do is to define
+`SciMLBase.solve(cache::LinearCache, alg::MyLUFactorization)`. Many algorithms
+may use a lazy matrix-free representation of the operator `A`. Thus if the
+algorithm requires a concrete matrix, like LU-factorization does, the algorithm
+should `convert(AbstractMatrix,cache.A)`. The flag `cache.isfresh` states whether
+`A` has changed since the last `solve`. Since we only need to factorize when
+`A` is new, the factorization part of the algorithm is done in a `if cache.isfresh`.
+`cache = set_cacheval(cache, fact)` puts the new factorization into the cache
+so it's updated for future solves. Then `y = ldiv!(cache.u, cache.cacheval, cache.b)`
+performs the solve and a linear solution is returned via
+`SciMLBase.build_linear_solution(alg,y,nothing,cache)`.

docs/src/basics/Preconditioners.md

@@ -1,4 +1,4 @@
-# Preconditioners
+# [Preconditioners](@id prec)
 
 Many linear solvers can be accelerated by using what is known as a **preconditioner**,
 an approximation to the matrix inverse action which is cheap to evaluate. These

docs/src/basics/common_solver_opts.md

@@ -0,0 +1,27 @@
+# Common Solver Options (Keyword Arguments to Solve)
+
+While many algorithms have specific arguments within their constructor,
+the keyword arguments for `solve` are common across all of the algorithms
+in order to give composability. These are also the options taken at `init` time.
+The following are the options these algorithms take, along with their defaults.
+
+## General Controls
+
+- `alias_A`: Whether to alias the matrix `A` or use a copy by default. When true,
+  algorithms like LU-factorization can be faster by reusing the memory via `lu!`,
+  but care must be taken as the original input will be modified. Default is `false`.
+- `alias_b`: Whether to alias the matrix `b` or use a copy by default. When true,
+  algorithms can write and change `b` upon usage. Care must be taken as the
+  original input will be modified. Default is `false`.
+- `verbose`: Whether to print extra information. Defaults to `false`.
+
+## Iterative Solver Controls
+
+Error controls are not used by all algorithms. Specifically, direct solves always
+solve completely. Error controls only apply to iterative solvers.
+
+- `abstol`: The absolute tolerance. Defaults to `√(eps(eltype(A)))`
+- `reltol`: The relative tolerance. Defaults to `√(eps(eltype(A)))`
+- `maxiters`: The number of iterations allowed. Defaults to `length(prob.b)`
+- `Pl,Pr`: The left and right preconditioners respectively. For more information
+  see [the Preconditioners page](@ref prec).