Add a new getting started tutorial and make operator updating safer

Fixes SciML/OrdinaryDiffEq.jl#2661 and uses that in a tutorial.

Author: ChrisRackauckas

docs/pages.jl

@@ -1,12 +1,11 @@
 pages = [
     "Home" => "index.md",
-    "sciml.md",
+    "tutorials/getting_started.md",
+    "Tutorials" => Any[
+        "operator_algebras.md",
+        "FFT Tutorial" => "tutorials/fftw.md"
+    ],
     "interface.md",
     "Premade Operators" => "premade_operators.md",
-    "Tutorials" => Any["FFT Tutorial" => "tutorials/fftw.md"
-    # "tutorials/linear.md",
-    # "tutorials/nonlin.md",
-    # "tutorials/ode.md",
-    # "tutorials/lux.md",
-    ]
+
 ]

docs/src/index.md

@@ -1,4 +1,4 @@
-# SciMLOperators.jl: Unified operator interface for `SciML.ai` and beyond
+# SciMLOperators.jl: Unified operator interface for Julia and SciML
 
 `SciMLOperators` is a package for managing linear, nonlinear,
 time-dependent, and parameter dependent operators acting on vectors,
@@ -25,69 +25,51 @@ using Pkg
 Pkg.add("SciMLOperators")
 ```
 
-## Examples
-
-Let `M`, `D`, `F` be matrix-based, diagonal-matrix-based, and function-based
-`SciMLOperators` respectively.
-
-```julia
-N = 4
-f = (v, u, p, t) -> u .* v
-
-M = MatrixOperator(rand(N, N))
-D = DiagonalOperator(rand(N))
-F = FunctionOperator(f, zeros(N), zeros(N))
-```
-
-Then, the following codes just work.
-
-```julia
-L1 = 2M + 3F + LinearAlgebra.I + rand(N, N)
-L2 = D * F * M'
-L3 = kron(M, D, F)
-L4 = M \ D
-L5 = [M; D]' * [M F; F D] * [F; D]
-```
-
-Each `L#` can be applied to `AbstractVector`s of appropriate sizes:
-
-```julia
-p = nothing # parameter struct
-t = 0.0     # time
-
-u = rand(N)
-v = rand(N)
-w = L1(v, u, p, t) # == L1 * v
-
-v_kron = rand(N^3)
-w_kron = L3(v_kron, u, p, t) # == L3 * v_kron
-```
-
-For mutating operator evaluations, call `cache_operator` to generate an
-in-place cache, so the operation is nonallocating.
-
-```julia
-α, β = rand(2)
-
-# allocate cache
-L2 = cache_operator(L2, u)
-L4 = cache_operator(L4, u)
-
-# allocation-free evaluation
-L2(w, v, u, p, t) # == mul!(w, L2, v)
-L4(w, v, u, p, t, α, β) # == mul!(w, L4, v, α, β)
-```
-
-The calling signature `L(v, u, p, t)`, for out-of-place evaluations, is
-equivalent to `L * v`, and the in-place evaluation `L(w, v, u, p, t, args...)`
-is equivalent to `LinearAlgebra.mul!(w, L, v, args...)`, where the arguments
-`u, p, t` are passed to `L` to update its state. More details are provided
-in the operator update section below.
-
-The `(v, u, p, t)` calling signature is standardized over the `SciML`
-ecosystem and is flexible enough to support use cases such as time-evolution
-in ODEs, as well as sensitivity computation with respect to the parameter
-object `p`.
+## Why `SciMLOperators`?
+
+Many functions, from linear solvers to differential equations, require
+the use of matrix-free operators to achieve maximum performance in
+many scenarios. `SciMLOperators.jl` defines the abstract interface for how
+operators in the SciML ecosystem are supposed to be defined. It gives the
+common set of functions and traits that solvers can rely on for properly
+performing their tasks. Along with that, `SciMLOperators.jl` provides
+definitions for the basic standard operators that are used as building
+blocks for most tasks, simplifying the use of operators while also
+demonstrating to users how such operators can be built and used in practice.
+
+`SciMLOperators.jl` has the design that is required to be used in
+all scenarios of equation solvers. For example, Magnus integrators for
+differential equations require defining an operator ``u' = A(t) u``, while
+Munthe-Kaas methods require defining operators of the form ``u' = A(u) u``.
+Thus, the operators need some form of time and state dependence, which the
+solvers can update and query when they are non-constant
+(`update_coefficients!`). Additionally, the operators need the ability to
+act like “normal” functions for equation solvers. For example, if `A(v,u,p,t)`
+has the same operation as `update_coefficients(A, u, p, t); A * v`, then `A`
+can be used in any place where a differential equation definition
+`(u,p,t) -> A(u, u, p, t)` is used without requiring the user or solver to do any extra
+work. 
+
+Another example is state-dependent mass matrices. `M(u,p,t)*u' = f(u,p,t)`.
+When solving such an equation, the solver must understand how to "update M"
+during operations, and thus the ability to update the state of `M` is a required
+function in the interface. This is also required for the definition of Jacobians
+`J(u,p,t)` in order to be properly used with Krylov methods inside of ODE solves
+without reconstructing the matrix-free operator at each step.
+
+Thus while previous good efforts for matrix-free operators have existed
+in the Julia ecosystem, such as
+[LinearMaps.jl](https://github.com/JuliaLinearAlgebra/LinearMaps.jl), those
+operator interfaces lack these aspects to actually be fully seamless
+with downstream equation solvers. This necessitates the definition and use of
+an extended operator interface with all of these properties, hence the
+`AbstractSciMLOperator` interface.
+
+!!! warn
+
+    This means that LinearMaps.jl is fundamentally lacking and is incompatible
+    with many of the tools in the SciML ecosystem, except for the specific cases
+    where the matrix-free operator is a constant!
 
 ## Features
 

docs/src/interface.md

@@ -1,4 +1,4 @@
-# The `AbstractSciMLOperator` Interface
+# [The `AbstractSciMLOperator` Interface](@id operator_interface)
 
 ```@docs
 SciMLOperators.AbstractSciMLOperator

docs/src/premade_operators.md

@@ -1,4 +1,4 @@
-# Premade SciMLOperators
+# [Premade SciMLOperators](@id premade_operators)
 
 ## Direct Operator Definitions
 

docs/src/sciml.md

@@ -1,4 +1,4 @@
-# Wrap a Fourier transform with SciMLOperators
+# [Wrap a Fourier transform with SciMLOperators](@id fft)
 
 In this tutorial, we will wrap a Fast Fourier Transform (FFT) in a SciMLOperator via the
 `FunctionOperator` interface. FFTs are commonly used algorithms for performing numerical