Premade SciMLOperators

IdentityOperator(len)

Matrix-free identity operator of size (len, len).

Arguments

• len::Integer: Number of rows and columns.

Fields

• len::Int: Stored operator dimension.

Interface Rules

IdentityOperator is constant, square, invertible, and supports the generic AbstractSciMLOperator application, caching, and trait interfaces. It is returned unchanged when composed with a compatible operator.

Examples

using SciMLOperators

IdentityOperator(2) * [3.0, 4.0] == [3.0, 4.0]

NullOperator(M, N)
NullOperator(N)

Matrix-free zero operator of size (M, N). The one-argument constructor creates a square zero operator.

Arguments

• M::Integer: Number of output rows.
• N::Integer: Number of input columns.

Fields

• M::Int: Stored output dimension.
• N::Int: Stored input dimension.

Interface Rules

The operator is constant and linear. It returns a zero result with the input container's element type and batch shape, and composition with a compatible operator remains a NullOperator with the composed dimensions.

Examples

using SciMLOperators

NullOperator(2, 3) * ones(3) == zeros(2)

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.

Arguments

• val::Number: Current scalar value.

Keyword Arguments

• update_func: Out-of-place update with signature update_func(oldval, u, p, t; kwargs...) -> newval.
• accepted_kwargs: Val tuple of keyword names forwarded to update_func.

Fields

• val: Current scalar value.
• update_func: Out-of-place scalar update function.

Interface Rules

Lazy scalar algebra is defined for AbstractSciMLScalarOperators. The interface supports lazy addition, subtraction, multiplication, and division. Updates must return a number with the intended scalar action; in-place scalar updates are not supported because numbers are immutable.

Examples

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.

Arguments

• A::AbstractMatrix: Matrix used for the current operator action.

Keyword Arguments

• update_func: Out-of-place update with signature update_func(A, u, p, t; kwargs...) -> new_A.
• update_func!: In-place update with signature update_func!(A, u, p, t; kwargs...).
• accepted_kwargs: Val tuple of keyword names forwarded to the update function.

Fields

• A: Current matrix state.
• update_func: Out-of-place matrix update function.
• update_func!: In-place matrix update function.

Interface Rules

Lazy matrix algebra is defined for AbstractSciMLOperators. The Interface supports lazy addition, subtraction, multiplication, inversion, adjoints, and transposes. An update function must preserve the dimensions and mathematical meaning of A; use FunctionOperator when the action itself is nonlinear in the input vector.

Examples

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.

Arguments

• diag::AbstractVector: Diagonal entries of the operator.

Keyword Arguments

• update_func: Out-of-place update with signature update_func(diag, u, p, t; kwargs...) -> new_diag.
• update_func!: In-place update with signature update_func!(diag, u, p, t; kwargs...).
• accepted_kwargs: Val tuple of keyword names forwarded to the update function.

Interface Rules

Updates must preserve the diagonal's leading dimension. For multidimensional diag, the operator acts elementwise while reporting a matrix size based on its leading dimension.

Examples

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.

Arguments

• A: Matrix or SciML operator applied to the action vector.
• B: Matrix or SciML operator applied to the additive term.
• b::AbstractArray: Additive input to B.

Keyword Arguments

• update_func: Out-of-place update for b.
• update_func!: In-place update for b.
• accepted_kwargs: Val tuple of keyword names forwarded to the update function.

Fields

• A: Linear action component.
• B: Additive-term action component.
• b: Current additive input.
• update_func: Out-of-place update for b.
• update_func!: In-place update for b.

Interface Rules

AffineOperator is not linear and cannot generally be converted to an AbstractMatrix. Updates must preserve compatible dimensions for A, B, and b; caching must preserve the affine action.

Examples

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(L, 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)

Construct the affine operation v + b as an AffineOperator.

Arguments

• b::AbstractVecOrMat: Additive term.

Keyword Arguments

• update_func: Out-of-place update for b with signature update_func(b, u, p, t; kwargs...) -> new_b.
• update_func!: In-place update for b with signature update_func!(b, u, p, t; kwargs...).
• accepted_kwargs: Val tuple of forwarded update keyword names.

Interface Rules

This is shorthand for AffineOperator(I, I, b; kwargs...). The term b must have a leading dimension compatible with the action vector. See AffineOperator for update and caching rules.

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

Construct the affine operation v + B * b as an AffineOperator.

Arguments

• B: Matrix or SciML operator acting on the additive term.
• b::AbstractVecOrMat: Additive input to B.

Keyword Arguments

• update_func: Out-of-place update for b.
• update_func!: In-place update for b.
• accepted_kwargs: Val tuple of forwarded update keyword names.

Interface Rules

B * b must have the same leading dimension as the action vector. This is shorthand for AffineOperator(I, B, b; kwargs...); see AffineOperator for the update and caching contract.

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}

AddedScalarOperator(α, β, ...)

Lazy sum of scalar operators.

Arguments

• α, β, ...: One or more AbstractSciMLScalarOperators.

Fields

• ops: Tuple of component scalar operators.

Interface Rules

Construct through scalar addition. The current scalar value is the sum of the updated component values, so updates and division traits are evaluated componentwise.

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

ComposedScalarOperator(α, β, ...)

Lazy product of scalar operators.

Arguments

• α, β, ...: One or more AbstractSciMLScalarOperators.

Fields

• ops: Tuple of component scalar operators.

Interface Rules

Construct through scalar multiplication or composition. Updates are forwarded to every component and the converted scalar is their product. Division is available only when every component supports it in the current state.

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

InvertedScalarOperator(α)

Lazy reciprocal of a scalar operator.

Arguments

• α::AbstractSciMLScalarOperator: Scalar operator to invert.

Fields

• λ: Wrapped scalar operator.

Interface Rules

Construct through inv(α). The current scalar value must be nonzero whenever the reciprocal is evaluated; updates are forwarded to λ before conversion.

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

ScaledOperator(λ, L)

Lazy scalar multiple of an AbstractSciMLOperator, representing λ * L.

Arguments

• λ: A number, UniformScaling, or AbstractSciMLScalarOperator.
• L::AbstractSciMLOperator: Operator being scaled.

Fields

• λ: Lazy scalar factor.
• L: Wrapped operator.

Interface Rules

Updates, caching, and trait queries are forwarded to both fields. The result is linear exactly when L is linear; division traits additionally require a nonzero current scalar value. Prefer λ * L to constructing this type directly.

AddedOperator(A, B, ...)

Lazy sum of compatible AbstractSciMLOperators.

Arguments

• A, B, ...: One or more operators with identical sizes.

Fields

• ops: Tuple of component operators. Nested sums are flattened.

Interface Rules

The action is the sum of each component action. Updates, caching, and traits are forwarded componentwise; a trait is true only when every required component supports it. Use A + B rather than constructing this type directly.

ComposedOperator(A, B, ...)

Lazy composition representing A * B * ..., with the rightmost operator applied first.

Arguments

• A, B, ...: One or more dimension-compatible operators.

Fields

• ops: Tuple of component operators in multiplication order.
• cache: Intermediate arrays used by in-place multiplication, or nothing.

Interface Rules

Updates are forwarded to every component. Call cache_operator before repeated in-place application when intermediate storage is required. The composition is linear, convertible, or supports a trait only when all components satisfy the corresponding contract.

InvertedOperator(L)

Lazy inverse of an AbstractSciMLOperator.

Arguments

• L::AbstractSciMLOperator: Operator whose current action is inverted.

Fields

• L: Wrapped operator.
• cache: Optional work arrays used by in-place division.

Interface Rules

Construct this type through inv(L). The wrapped operator must support the division operations required by the chosen application form. Updates and caching are forwarded to L; inversion does not materialize a matrix.

InvertibleOperator(L, F)

Pair an operator with a factorization or inverse object used for division.

Arguments

• L: Original operator used for multiplication and updates.
• F: Factorization or inverse supporting \ or ldiv!.

Fields

• L: Wrapped operator.
• F: Factorization or inverse representation.

Interface Rules

Constructors such as factorize, lu, and qr create this wrapper for concretizable operators. L and F must represent the same current action; after a state update, extension code is responsible for keeping both current. Use ldiv! or \ for the inverse action rather than assuming F is a public field contract.

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

AdjointOperator(L)

Lazy adjoint wrapper for an AbstractSciMLOperator.

Arguments

• L::AbstractSciMLOperator: Operator to adjoint.

Fields

• L: Wrapped operator.

Interface Rules

Construct through adjoint(L) or L'. The wrapper delegates updates and matrix-like application to the adjoint action of L; it is valid only when has_adjoint(L) is true.

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

TransposedOperator(L)

Lazy transpose wrapper for an AbstractSciMLOperator.

Arguments

• L::AbstractSciMLOperator: Operator to transpose.

Fields

• L: Wrapped operator.

Interface Rules

Construct through transpose(L). The wrapper delegates updates and matrix-like application to the transpose action of L; complex operators must distinguish this from adjoint(L).