Merge pull request #275 from SciML/operator_algebras

Fix operator algebras tutorial

Author: ChrisRackauckas

README.md

@@ -41,52 +41,61 @@ julia> Pkg.add("SciMLOperators")
 Let `M`, `D`, `F` be matrix-based, diagonal-matrix-based, and function-based
 `SciMLOperators` respectively.
 
-```julia
+Let `M`, `D`, `F` be matrix-based, diagonal-matrix-based, and function-based
+`SciMLOperators` respectively.
+
+```@example operator_algebra
+using SciMLOperators, LinearAlgebra
 N = 4
-f(u, p, t) = u .* u
-f(v, u, p, t) = v .= u .* u
+function f(v, u, p, t) 
+    u .* v
+end
+function f(w, v, u, p, t) 
+    w .= u .* v
+end
+
+u = rand(4)
+p = nothing # parameter struct
+t = 0.0     # time
 
 M = MatrixOperator(rand(N, N))
 D = DiagonalOperator(rand(N))
-F = FunctionOperator(f, zeros(N), zeros(N))
+F = FunctionOperator(f, zeros(N), zeros(N); u, p, t)
 ```
 
 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
-p = nothing # parameter struct
-t = 0.0     # time
-
-u = rand(N)
-v = L1(u, p, t) # == L1 * u
+```@example operator_algebra
+v = rand(N)
+w = L1(v, u, p, t) # == L1 * v
 
-u_kron = rand(N^3)
-v_kron = L3(u_kron, p, t) # == L3 * u_kron
+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
-in-place cache so the operation is nonallocating.
+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
 L2 = cache_operator(L2, u)
 L4 = cache_operator(L4, u)
 
 # allocation-free evaluation
-L2(v, u, p, t) # == mul!(v, L2, u)
-L4(v, u, p, t, α, β) # == mul!(v, L4, u, α, β)
+L2(w, v, u, p, t) # == mul!(w, L2, v)
+L4(w, v, u, p, t, α, β) # == mul!(w, L4, v, α, β)
 ```
 
 ## Roadmap

docs/src/tutorials/operator_algebras.md

@@ -1,34 +1,41 @@
 ## [Demonstration of Operator Algebras and Kron](@id operator_algebras)
 
 Let `M`, `D`, `F` be matrix-based, diagonal-matrix-based, and function-based
-`SciMLOperators` respectively.
+`SciMLOperators` respectively. Here are some examples of composing operators
+in order to build more complex objects and using their operations.
 
-```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
+
+u = rand(4)
+p = nothing # parameter struct
+t = 0.0     # time
 
 M = MatrixOperator(rand(N, N))
 D = DiagonalOperator(rand(N))
-F = FunctionOperator(f, zeros(N), zeros(N))
+F = FunctionOperator(f, zeros(N), zeros(N); u, p, t)
 ```
 
 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
-p = nothing # parameter struct
-t = 0.0     # time
-
-u = rand(N)
+```@example operator_algebra
 v = rand(N)
 w = L1(v, u, p, t) # == L1 * v
 
@@ -39,7 +46,7 @@ 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
@@ -49,15 +56,4 @@ 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`.
+```

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

src/tensor.jl

@@ -172,7 +172,7 @@ function Base.:*(L::TensorProductOperator, v::AbstractVecOrMat)
     k = size(v, 2)
 
     U = reshape(v, (ni, no * k))
-    C = inner * U
+    C = stack([inner * _v for _v in eachcol(U)])
 
     V = outer_mul(L, v, C)