Demonstration of Operator Algebras and Kron
Let M
, D
, F
be matrix-based, diagonal-matrix-based, and function-based SciMLOperators
respectively. Here are some examples of composing operators in order to build more complex objects and using their operations.
using SciMLOperators, LinearAlgebra
N = 4
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); u, p, t)
FunctionOperator(4 × 4)
Then, the following codes just work.
L1 = 2M + 3F + LinearAlgebra.I + rand(N, N)
L2 = D * F * M'
L3 = kron(M, D, F)
L4 = lu(M) \ D
L5 = [M; D]' * [M F; F D] * [F; D]
((((MatrixOperator(4 × 4) * MatrixOperator(4 × 4)) + (DiagonalOperator(4 × 4) * FunctionOperator(4 × 4))) * FunctionOperator(4 × 4)) + (((MatrixOperator(4 × 4) * FunctionOperator(4 × 4)) + (DiagonalOperator(4 × 4) * DiagonalOperator(4 × 4))) * DiagonalOperator(4 × 4)))
Each L#
can be applied to AbstractVector
s of appropriate sizes:
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
64-element reshape(transpose(::Matrix{Float64}), 64) with eltype Float64:
0.4028108527772814
0.3647533377065399
0.5218260337786796
0.2636762462821975
0.6530901434972598
0.5757734126287503
0.22665739417762296
0.1865263437871233
0.3473670511523424
0.43430094753753484
⋮
0.3074739197068582
0.415863021185039
0.516952544209508
0.4086682974715123
0.5280402537666745
0.4452563879940305
0.6630309309569733
0.9388202728434266
0.4379340659090096
For mutating operator evaluations, call cache_operator
to generate an in-place cache, so the operation is nonallocating.
α, β = 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, α, β)
4-element Vector{Float64}:
0.43962324365104816
0.5824609006732383
0.2716601765576226
0.5262672723133085