Fix operator algebras tutorial

Author: ChrisRackauckas

docs/src/tutorials/operator_algebras.md

@@ -3,9 +3,15 @@
 Let `M`, `D`, `F` be matrix-based, diagonal-matrix-based, and function-based
 `SciMLOperators` respectively.
 
-```julia
+```@example operator_algebra
+using SciMLOperators, LinearAlgebra
 N = 4
-f = (v, u, p, t) -> u .* v
+function f(v, u, p, t) 
+    u .* v
+end
+function f(w, v, u, p, t) 
+    w .= u .* v
+end
 
 M = MatrixOperator(rand(N, N))
 D = DiagonalOperator(rand(N))
@@ -14,32 +20,34 @@ F = FunctionOperator(f, zeros(N), zeros(N))
 
 Then, the following codes just work.
 
-```julia
+```@example operator_algebra
 L1 = 2M + 3F + LinearAlgebra.I + rand(N, N)
 L2 = D * F * M'
 L3 = kron(M, D, F)
-L4 = M \ D
+L4 = lu(M) \ D
 L5 = [M; D]' * [M F; F D] * [F; D]
 ```
 
 Each `L#` can be applied to `AbstractVector`s of appropriate sizes:
 
-```julia
+```@example operator_algebra
 p = nothing # parameter struct
 t = 0.0     # time
 
 u = rand(N)
 v = rand(N)
 w = L1(v, u, p, t) # == L1 * v
 
+L3 = cache_operator(L3, 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
+```@example operator_algebra
 α, β = rand(2)
 
 # allocate cache

src/func.jl

@@ -734,6 +734,17 @@ function Base.:*(L::FunctionOperator{iip, true}, v::AbstractArray) where {iip}
     vec_output ? vec(W) : W
 end
 
+function Base.:*(L::FunctionOperator{iip, false}, v::AbstractArray) where {iip}
+    _sizecheck(L, v, nothing)
+    V, _, vec_output = _unvec(L, v, nothing)
+
+    W, _ = L.cache
+    W = copy(W)
+    L.op(W, V, L.u, L.p, L.t; L.traits.kwargs...)
+
+    vec_output ? vec(W) : W
+end
+
 function Base.:\(L::FunctionOperator{iip, true}, v::AbstractArray) where {iip}
     _sizecheck(L, nothing, v)
     _, V, vec_output = _unvec(L, nothing, v)
@@ -743,16 +754,30 @@ function Base.:\(L::FunctionOperator{iip, true}, v::AbstractArray) where {iip}
     vec_output ? vec(W) : W
 end
 
+function Base.:\(L::FunctionOperator{iip, false}, v::AbstractArray) where {iip}
+    _sizecheck(L, nothing, v)
+    _, V, vec_output = _unvec(L, nothing, v)
+
+    W, _ = L.cache
+    W = copy(W)
+    L.op_inverse(W, V, L.u, L.p, L.t; L.traits.kwargs...)
+
+    vec_output ? vec(W) : W
+end
+
 function LinearAlgebra.mul!(w::AbstractArray, L::FunctionOperator{true}, v::AbstractArray)
     _sizecheck(L, v, w)
     V, W, vec_output = _unvec(L, v, w)
     L.op(W, V, L.u, L.p, L.t; L.traits.kwargs...)
     vec_output ? vec(W) : W
 end
 
-function LinearAlgebra.mul!(::AbstractArray, L::FunctionOperator{false}, ::AbstractArray,
+function LinearAlgebra.mul!(w::AbstractArray, L::FunctionOperator{false}, ::AbstractArray,
         args...)
-    @error "LinearAlgebra.mul! not defined for out-of-place operator $L"
+    _sizecheck(L, v, w)
+    V, W, vec_output = _unvec(L, v, w)
+    W .= L.op(V, L.u, L.p, L.t; L.traits.kwargs...)
+    vec_output ? vec(W) : W
 end
 
 function LinearAlgebra.mul!(w::AbstractArray, L::FunctionOperator{true, oop, false},
@@ -797,15 +822,26 @@ function LinearAlgebra.ldiv!(L::FunctionOperator{true}, v::AbstractArray)
     copy!(W, V)
     L.op_inverse(W, V, L.u, L.p, L.t; L.traits.kwargs...) # ldiv!(U, L, V)
 
-    vec_output ? vec(W) : W
+    V .= W
+    vec_output ? vec(V) : V
 end
 
-function LinearAlgebra.ldiv!(v::AbstractArray, L::FunctionOperator{false}, u::AbstractArray)
-    @error "LinearAlgebra.ldiv! not defined for out-of-place $L"
+function LinearAlgebra.ldiv!(w::AbstractArray, L::FunctionOperator{false}, v::AbstractArray)
+    _sizecheck(L, v, w)
+    W, V, _ = _unvec(L, w, v)
+
+    W .= L.op_inverse(V, L.u, L.p, L.t; L.traits.kwargs...)
+
+    w
 end
 
-function LinearAlgebra.ldiv!(L::FunctionOperator{false}, u::AbstractArray)
-    @error "LinearAlgebra.ldiv! not defined for out-of-place $L"
+function LinearAlgebra.ldiv!(L::FunctionOperator{false}, v::AbstractArray)
+    _sizecheck(L, nothing, v)
+    V, _, vec_output = _unvec(L, v, nothing)
+
+    V .= L.op_inverse(V, L.u, L.p, L.t; L.traits.kwargs...) # ldiv!(W, L, V)
+
+    vec_output ? vec(V) : V
 end
 
 # Out-of-place: v is action vector, u is update vector