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

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

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

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}