Premade SciMLOperators

struct IdentityOperator <: SciMLOperators.AbstractSciMLOperator{Bool}

Operator representing the identity function id(v) = v

struct NullOperator <: SciMLOperators.AbstractSciMLOperator{Bool}

Operator representing the null function n(v) = 0 * v

ScalarOperator(val; update_func, accepted_kwargs)

Represents a linear scaling operator that may be applied to a Number, or an AbstractArray subtype. Its state is updated by the user-provided update_func during operator evaluation (L([w,] v, u, p, t)), or by calls to update_coefficients[!]. Both recursively call the update function, update_func which is assumed to have the signature:

update_func(oldval::Number, u, p, t; <accepted kwargs>) -> newval

The set of keyword-arguments accepted by update_func must be provided to ScalarOperator via the kwarg accepted_kwargs as a tuple of Symbols. kwargs cannot be passed down to update_func if accepted_kwargs are not provided.

Interface

Lazy scalar algebra is defined for AbstractSciMLScalarOperators. The interface supports lazy addition, subtraction, multiplication and division.

Example

v = rand(4)
u = rand(4)
w = zeros(4)
p = nothing
t = 0.0

val_update = (a, u, p, t; scale = 0.0) -> scale
α = ScalarOperator(0.0; update_func = val_update, accepted_kwargs = (:scale,))
β = 2 * α + 3 / α

# Update β and evaluate with the new interface
result = β(v, u, p, t; scale = 1.0)

# In-place application
β(w, v, u, p, t; scale = 1.0)

# In-place with scaling
w_orig = copy(w)
α_val = 2.0
β_val = 0.5
β(w, v, u, p, t, α_val, β_val; scale = 1.0) # w = α_val*(β*v) + β_val*w

Represents a linear operator given by an AbstractMatrix that may be applied to an AbstractVecOrMat. Its state is updated by the user-provided update_func during operator evaluation (L([w,], v, u, p, t)), or by calls to update_coefficients[!](L, u, p, t). Both recursively call the update_function, update_func which is assumed to have the signature

update_func(A::AbstractMatrix, u, p, t; <accepted kwargs>) -> newA

or

update_func!(A::AbstractMatrix, u, p, t; <accepted kwargs>) -> [modifies A]

The set of keyword-arguments accepted by update_func[!] should be provided to MatrixOperator via the kwarg accepted_kwargs as a Val of a tuple of Symbols for zero-allocation kwarg filtering. For example, accepted_kwargs = Val((:dtgamma,)). Plain tuples like (:dtgamma,) are deprecated but still supported. kwargs cannot be passed down to update_func[!] if accepted_kwargs are not provided.

Interface

Lazy matrix algebra is defined for AbstractSciMLOperators. The Interface supports lazy addition, subtraction, multiplication, inversion, adjoints, transposes.

Example

Out-of-place update and usage

v = rand(4)
u = rand(4)
p = rand(4, 4)
t = rand()

mat_update = (A, u, p, t; scale = 0.0) -> t * p
M = MatrixOperator(0.0; update_func = mat_update, accepted_kwargs = Val((:scale,)))

L = M * M + 3I
L = cache_operator(L, v)

# update and evaluate 
w = L(v, u, p, t; scale = 1.0)

# In-place evaluation
w = similar(v)
L(w, v, u, p, t; scale = 1.0)

# In-place with scaling
β = 0.5
L(w, v, u, p, t, 2.0, β; scale = 1.0) # w = 2.0*(L*v) + 0.5*w

In-place update and usage

w = zeros(4)
v = zeros(4)
u = rand(4)
p = rand(4) # Must be non-nothing
t = rand()

mat_update! = (A, u, p, t; scale = 0.0) -> (A .= t * p * u' * scale)
M = MatrixOperator(zeros(4, 4); update_func! = mat_update!, accepted_kwargs = Val((:scale,)))
L = M * M + 3I
L = cache_operator(L, v) 

# update L in-place and evaluate
update_coefficients!(L, u, p, t; scale = 1.0)
mul!(w, L, v)

# Or use the new interface that separates update and application
L(w, v, u, p, t; scale = 1.0)

DiagonalOperator(
    diag;
    update_func,
    update_func!,
    accepted_kwargs
)

Represents an elementwise scaling (diagonal-scaling) operation that may be applied to an AbstractVecOrMat. When diag is an AbstractVector of length N, L = DiagonalOperator(diag, ...) can be applied to AbstractArrays with size(u, 1) == N. Each column of the v will be scaled by diag, as in LinearAlgebra.Diagonal(diag) * v.

When diag is a multidimensional array, L = DiagonalOperator(diag, ...) forms an operator of size (N, N) where N = size(diag, 1) is the leading length of diag. L then is the elementwise-scaling operation on arrays of length(v) = length(diag) with leading length size(u, 1) = N.

Its state is updated by the user-provided update_func during operator evaluation (L([w,], v, u, p, t)), or by calls to update_coefficients[!](L, u, p, t). Both recursively call the update_function, update_func which is assumed to have the signature

update_func(diag::AbstractVecOrMat, u, p, t; <accepted kwargs>) -> new_diag

or

update_func!(diag::AbstractVecOrMat, u, p, t; <accepted kwargs>) -> [modifies diag]

The set of keyword-arguments accepted by update_func[!] should be provided to DiagonalOperator via the kwarg accepted_kwargs as a Val of a tuple of Symbols for zero-allocation kwarg filtering. For example, accepted_kwargs = Val((:dtgamma,)). Plain tuples like (:dtgamma,) are deprecated but still supported. kwargs cannot be passed down to update_func[!] if accepted_kwargs are not provided.

Example

BatchedDiagonalOperator(diag; update_func, update_func!, accepted_kwargs)

Represents a time-dependent elementwise scaling (diagonal-scaling) operation. Acts on AbstractArrays of the same size as diag. The update function is called by update_coefficients! and is assumed to have the following signature:

update_func(diag::AbstractArray, u, p, t; <accepted kwarg fields>) -> [modifies diag]

Represents a generalized affine operation (w = A * v + B * b) that may be applied to an AbstractVecOrMat. The user-provided update functions, update_func[!] update the AbstractVecOrMat b, and are called during operator evaluation (L([w,], v, u, p, t)), or by calls to update_coefficients[!](L, u, p, t). The update functions are assumed to have the syntax

update_func(b::AbstractVecOrMat, u, p, t; <accepted kwargs>) -> new_b

or

update_func!(b::AbstractVecOrMat, u ,p , t; <accepted kwargs>) -> [modifies b]

and B, b are expected to have an appropriate size so that A * v + B * b makes sense. Specifically, size(A, 1) == size(B, 1), and size(v, 2) == size(b, 2).

The set of keyword-arguments accepted by update_func[!] should be provided to AffineOperator via the kwarg accepted_kwargs as a Val of a tuple of Symbols for zero-allocation kwarg filtering. For example, accepted_kwargs = Val((:dtgamma,)). Plain tuples like (:dtgamma,) are deprecated but still supported. kwargs cannot be passed down to update_func[!] if accepted_kwargs are not provided.

Example

v = rand(4)
u = rand(4)
p = rand(4)
t = rand()

A = MatrixOperator(rand(4, 4))
B = MatrixOperator(rand(4, 4))

vec_update_func = (b, u, p, t) -> p .* u * t
L = AffineOperator(A, B, zeros(4); update_func = vec_update_func)
L = cache_operator(M, v)

# update L and evaluate
w = L(v, u, p, t) # == A * v + B * (p .* u * t)

AddVector(b; update_func, update_func!, accepted_kwargs)

Represents the affine operation w = I * v + I * b. The update functions, update_func[!] update the state of AbstractVecOrMat b. See documentation of AffineOperator for more details.

AddVector(B, b; update_func, update_func!, accepted_kwargs)

Represents the affine operation w = I * v + B * b. The update functions, update_func[!] update the state of AbstractVecOrMat b. See documentation of AffineOperator for more details.

Matrix free operator given by a function

• op: Function with signature op(v, u, p, t) and (if isinplace) op(w, v, u, p, t)

• op_adjoint: Adjoint operator

• op_inverse: Inverse operator

• op_adjoint_inverse: Adjoint inverse operator

• traits: Traits

• u: State

• p: Parameters

• t: Time

• cache: Cache

struct BlockDiagonalOperator{T, O<:Tuple{Vararg{SciMLOperators.AbstractSciMLOperator}}} <: SciMLOperators.AbstractSciMLOperator{T}

Lazy block diagonal operator built from AbstractSciMLOperator blocks.

Arguments

• ops: Operators or matrices to place on the block diagonal. Matrix arguments are wrapped in MatrixOperator.

Fields

• ops

Interface Rules

BlockDiagonalOperator applies each block to the corresponding slice of the input and concatenates the results. Its size is the sum of block row and column sizes. update_coefficients[!], caching, and trait queries are forwarded to each block, so a block diagonal operator has concretization, in-place multiplication, or adjoint support only when the required component operators do.

Examples

using LinearAlgebra, SciMLOperators

A = MatrixOperator([1.0 2.0; 3.0 4.0])
B = MatrixOperator(Diagonal([5.0, 6.0, 7.0]))
L = BlockDiagonalOperator(A, B)

v = ones(5)
L * v == Matrix(L) * v

Computes the lazy pairwise Kronecker product, or tensor product, operator of AbstractMatrix, and AbstractSciMLOperator subtypes. Calling ⊗(ops...) is equivalent to Base.kron(ops...). Fast operator evaluation is performed without forming the full tensor product operator.

TensorProductOperator(A, B) = A ⊗ B
TensorProductOperator(A, B, C) = A ⊗ B ⊗ C

(A ⊗ B)(v) = vec(B * reshape(v, M, N) * transpose(A))

where M = size(B, 2), and N = size(A, 2)

Example

using SciMLOperators, LinearAlgebra

# Create basic operators
A = rand(3, 3)
B = rand(4, 4)
A_op = MatrixOperator(A)
B_op = MatrixOperator(B)

# Create tensor product operator
T = A_op ⊗ B_op

# Apply to a vector using the new interface
v = rand(3*4)    # Action vector
u = rand(3*4)    # Update vector
p = nothing
t = 0.0

# Out-of-place application
result = T(v, u, p, t)

# For in-place operations, need to cache the operator first
T_cached = cache_operator(T, v)

# In-place application
w = zeros(size(T, 1))
T_cached(w, v, u, p, t)

# In-place with scaling
w_orig = copy(w)
α = 2.0
β = 0.5
T_cached(w, v, u, p, t, α, β) # w = α*(T*v) + β*w_orig

⊗(ops)

Computes the lazy pairwise Kronecker product, or tensor product, operator of AbstractMatrix, and AbstractSciMLOperator subtypes. Calling ⊗(ops...) is equivalent to Base.kron(ops...). Fast operator evaluation is performed without forming the full tensor product operator.

TensorProductOperator(A, B) = A ⊗ B
TensorProductOperator(A, B, C) = A ⊗ B ⊗ C

(A ⊗ B)(v) = vec(B * reshape(v, M, N) * transpose(A))

where M = size(B, 2), and N = size(A, 2)

Example

using SciMLOperators, LinearAlgebra

# Create basic operators
A = rand(3, 3)
B = rand(4, 4)
A_op = MatrixOperator(A)
B_op = MatrixOperator(B)

# Create tensor product operator
T = A_op ⊗ B_op

# Apply to a vector using the new interface
v = rand(3*4)    # Action vector
u = rand(3*4)    # Update vector
p = nothing
t = 0.0

# Out-of-place application
result = T(v, u, p, t)

# For in-place operations, need to cache the operator first
T_cached = cache_operator(T, v)

# In-place application
w = zeros(size(T, 1))
T_cached(w, v, u, p, t)

# In-place with scaling
w_orig = copy(w)
α = 2.0
β = 0.5
T_cached(w, v, u, p, t, α, β) # w = α*(T*v) + β*w_orig

kron(A, B)

Construct a lazy representation of the Kronecker product A ⊗ B. One of the two factors can be an AbstractMatrix, which is then promoted to a MatrixOperator automatically. To avoid fallback to the generic Base.kron, at least one of A and B must be an AbstractSciMLOperator.

struct TensorSumOperator{T, O, P} <: SciMLOperators.AbstractSciMLOperator{T}

Lazy Kronecker sum operator.

Arguments

• outer: A square matrix or AbstractSciMLOperator representing the first term in outer ⊗ I.
• inner: A square matrix or AbstractSciMLOperator representing the second term in I ⊗ inner.

Fields

• ops

• products

Interface Rules

TensorSumOperator(outer, inner) represents outer ⊗ I + I ⊗ inner without eagerly forming the Kronecker products. Both input operators must be square. The operator forwards state updates to outer and inner, and its cached application stores the two tensor-product terms needed by mul!.

isconvertible(::TensorSumOperator) is false because eager fusion is not the default algebra path, but has_concretization(L) is true when both operands can be materialized.

Examples

using LinearAlgebra, SciMLOperators

A = MatrixOperator([1.0 2.0; 3.0 4.0])
B = MatrixOperator(Diagonal([5.0, 6.0, 7.0]))
L = TensorSumOperator(A, B)

v = ones(6)
L * v == Matrix(L) * v

kronsum(A, B)

Construct the lazy Kronecker sum A ⊗ I + I ⊗ B.

Arguments

• A: A square matrix or AbstractSciMLOperator.
• B: A square matrix or AbstractSciMLOperator.

Returns

A TensorSumOperator whose action is equivalent to kron(A, I(size(B, 1))) + kron(I(size(A, 1)), B).

Interface Rules

Both inputs must be square. Matrix inputs are wrapped in MatrixOperator so the returned object participates in the AbstractSciMLOperator update, caching, multiplication, and trait interfaces.

Examples

using LinearAlgebra, SciMLOperators

A = [1.0 2.0; 3.0 4.0]
B = Diagonal([5.0, 6.0, 7.0])
L = kronsum(A, B)

v = ones(6)
L * v == Matrix(L) * v

mutable struct WOperator{IIP, T, MType, GType, JType, F, C, JV} <: SciMLOperators.AbstractWOperator{T}

WOperator{IIP}(mass_matrix, gamma, J, u[, jacvec])

A linear operator that represents the W matrix of an ODEProblem, defined as

\[W = \frac{1}{\gamma}MM - J\]

where MM is the mass matrix, γ is a scalar, and J is the Jacobian operator.

Arguments

• mass_matrix: A matrix-like object, UniformScaling, or MatrixOperator representing MM.
• gamma: Scalar coefficient in the W-operator definition.
• J: Jacobian represented as a number, matrix, or AbstractSciMLOperator.
• u: Prototype state used to allocate the internal multiplication cache.
• jacvec: Optional operator used for Jacobian-vector products in mul!.

Fields

• mass_matrix

• gamma

• J

• _func_cache

• _concrete_form

• jacvec

Interface Rules

WOperator is part of the public solver-developer interface used by implicit ODE solvers. It supports matrix-like *, \, mul!, indexing, sizing, and concretization. Calling update_coefficients!(W, u, p, t; gamma) updates the Jacobian, mass matrix, optional Jacobian-vector operator, and stored gamma. Omitting (u, p, t) leaves those operators unchanged and only updates gamma when it is supplied.

IIP controls whether conversion reuses the internally stored concrete form as an in-place operator. The public contract is the mathematical action of W; downstream code should not depend on _func_cache or _concrete_form.

Examples

using LinearAlgebra, SciMLOperators

J = MatrixOperator([1.0 2.0; 3.0 4.0])
W = WOperator{true}(I, 0.5, J, zeros(2))

v = [1.0, 2.0]
W * v == (Matrix(W) * v)

struct StaticWOperator{isinv, T, F} <: SciMLOperators.AbstractWOperator{T}

Small dense factorization helper for repeated solves with a fixed W matrix.

Arguments

• W: Concrete square matrix to solve against.
• callinv: Whether very small matrices may store inv(W) for direct multiplication during \.

Fields

• W

• F

Interface Rules

StaticWOperator is a specialized helper for solver internals that need a fixed W-operator solve. It supports Wstatic \ v; it does not participate in coefficient updates and should be reconstructed when the underlying matrix changes.

Examples

using SciMLOperators

W = StaticWOperator([2.0 0.0; 0.0 4.0])
W \ [2.0, 8.0]

struct AddedScalarOperator{T, O} <: SciMLOperators.AbstractSciMLScalarOperator{T}

Lazy addition of AbstractSciMLScalarOperators

struct ComposedScalarOperator{T, O} <: SciMLOperators.AbstractSciMLScalarOperator{T}

Lazy multiplication of AbstractSciMLScalarOperators

struct InvertedScalarOperator{T, λType} <: SciMLOperators.AbstractSciMLScalarOperator{T}

Lazy inverse of AbstractSciMLScalarOperators

struct ScaledOperator{T, λType, LType} <: SciMLOperators.AbstractSciMLOperator{T}

ScaledOperator

(λ L)*(v) = λ * L(v)

Lazy operator addition

(A1 + A2 + A3...)v = A1*v + A2*v + A3*v ....

Lazy operator composition

∘(A, B, C)(v) = A(B(C(v)))

ops = (A, B, C)
cache = (B*C*v , C*v)

Lazy Operator Inverse

Stores an operator and its factorization (or inverse operator). Supports left division and ldiv! via F, and operator evaluation via L.

struct AdjointOperator{T, LType} <: SciMLOperators.AbstractSciMLOperator{T}

struct TransposedOperator{T, LType} <: SciMLOperators.AbstractSciMLOperator{T}