diff --git a/README.md b/README.md index cf6a6245..a48dbcd8 100644 --- a/README.md +++ b/README.md @@ -47,10 +47,10 @@ Let `M`, `D`, `F` be matrix-based, diagonal-matrix-based, and function-based ```@example operator_algebra using SciMLOperators, LinearAlgebra N = 4 -function f(v, u, p, t) +function f(v, u, p, t) u .* v end -function f(w, v, u, p, t) +function f(w, v, u, p, t) w .= u .* v end diff --git a/docs/pages.jl b/docs/pages.jl index 7bea2ba9..e3bdbf54 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -6,6 +6,4 @@ pages = [ "FFT Tutorial" => "tutorials/fftw.md" ], "interface.md", - "Premade Operators" => "premade_operators.md", - -] + "Premade Operators" => "premade_operators.md"] diff --git a/docs/src/index.md b/docs/src/index.md index 7e171d21..bacbdcf5 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -48,7 +48,7 @@ act like “normal” functions for equation solvers. For example, if `A(v,u,p,t has the same operation as `update_coefficients(A, u, p, t); A * v`, then `A` can be used in any place where a differential equation definition `(u,p,t) -> A(u, u, p, t)` is used without requiring the user or solver to do any extra -work. +work. Another example is state-dependent mass matrices. `M(u,p,t)*u' = f(u,p,t)`. When solving such an equation, the solver must understand how to "update M" @@ -66,7 +66,7 @@ an extended operator interface with all of these properties, hence the `AbstractSciMLOperator` interface. !!! warn - + This means that LinearMaps.jl is fundamentally lacking and is incompatible with many of the tools in the SciML ecosystem, except for the specific cases where the matrix-free operator is a constant! diff --git a/docs/src/tutorials/getting_started.md b/docs/src/tutorials/getting_started.md index 38d0ca92..bb91c835 100644 --- a/docs/src/tutorials/getting_started.md +++ b/docs/src/tutorials/getting_started.md @@ -11,16 +11,16 @@ and updating operators. Before we get into the deeper operators, let's show the simplest SciMLOperator: `MatrixOperator`. `MatrixOperator` just turns a matrix into an `AbstractSciMLOperator`, so it's not really a matrix-free operator but it's a starting point that is good for -understanding the interface and testing. To create a `MatrixOperator`, simply call the +understanding the interface and testing. To create a `MatrixOperator`, simply call the constructor on a matrix: ```@example getting_started using SciMLOperators, LinearAlgebra -A = [-2.0 1 0 0 0 - 1 -2 1 0 0 - 0 1 -2 1 0 - 0 0 1 -2 1 - 0 0 0 1 -2] +A = [-2.0 1 0 0 0 + 1 -2 1 0 0 + 0 1 -2 1 0 + 0 0 1 -2 1 + 0 0 0 1 -2] opA = MatrixOperator(A) ``` @@ -29,7 +29,7 @@ The operators can do [operations as defined in the operator interface](@ref oper matrix multiplication as the core action: ```@example getting_started -v = [3.0,2.0,1.0,2.0,3.0] +v = [3.0, 2.0, 1.0, 2.0, 3.0] opA*v ``` @@ -43,7 +43,8 @@ mul!(w, opA, v) ``` ```@example getting_started -α = 1.0; β = 1.0 +α = 1.0; +β = 1.0 mul!(w, opA, v, α, β) # α*opA*v + β*w ``` @@ -64,18 +65,20 @@ For example, let's make the operator `A .* u + dt*I` where `dt` is a parameter and `u` is a state vector: ```@example getting_started -A = [-2.0 1 0 0 0 - 1 -2 1 0 0 - 0 1 -2 1 0 - 0 0 1 -2 1 - 0 0 0 1 -2] +A = [-2.0 1 0 0 0 + 1 -2 1 0 0 + 0 1 -2 1 0 + 0 0 1 -2 1 + 0 0 0 1 -2] function update_function!(B, u, p, t) dt = p B .= A .* u + dt*I end -u = Array(1:1.0:5); p = 0.1; t = 0.0 +u = Array(1:1.0:5); +p = 0.1; +t = 0.0 opB = MatrixOperator(copy(A); update_func! = update_function!) ``` @@ -123,18 +126,18 @@ With `FunctionOperator`, we directly define the operator application function `o which means `w = opA(u,p,t)*v`. For example we can do the following: ```@example getting_started -function Afunc!(w,v,u,p,t) +function Afunc!(w, v, u, p, t) w[1] = -2v[1] + v[2] for i in 2:4 - w[i] = v[i-1] - 2v[i] + v[i+1] + w[i] = v[i - 1] - 2v[i] + v[i + 1] end w[5] = v[4] - 2v[5] nothing end -function Afunc!(v,u,p,t) +function Afunc!(v, u, p, t) w = zeros(5) - Afunc!(w,v,u,p,t) + Afunc!(w, v, u, p, t) w end @@ -148,29 +151,29 @@ mfopA*v - opA*v ``` ```@example getting_started -mfopA(v,u,p,t) - opA(v,u,p,t) +mfopA(v, u, p, t) - opA(v, u, p, t) ``` We can also create the state-dependent operator as well: ```@example getting_started -function Bfunc!(w,v,u,p,t) +function Bfunc!(w, v, u, p, t) dt = p w[1] = -(2*u[1]-dt)*v[1] + v[2]*u[1] for i in 2:4 - w[i] = v[i-1]*u[i] - (2*u[i]-dt)*v[i] + v[i+1]*u[i] + w[i] = v[i - 1]*u[i] - (2*u[i]-dt)*v[i] + v[i + 1]*u[i] end w[5] = v[4]*u[5] - (2*u[5]-dt)*v[5] nothing end -function Bfunc!(v,u,p,t) +function Bfunc!(v, u, p, t) w = zeros(5) - Bfunc!(w,v,u,p,t) + Bfunc!(w, v, u, p, t) w end -mfopB = FunctionOperator(Bfunc!, zeros(5), zeros(5); u, p, t, isconstant=false) +mfopB = FunctionOperator(Bfunc!, zeros(5), zeros(5); u, p, t, isconstant = false) ``` ```@example getting_started @@ -187,19 +190,19 @@ operator for `A .* u` (since right now there is not a built in operator for vect but that would be a fantastic thing to add!): ```@example getting_started -function Cfunc!(w,v,u,p,t) +function Cfunc!(w, v, u, p, t) w[1] = -2v[1] + v[2] for i in 2:4 - w[i] = v[i-1] - 2v[i] + v[i+1] + w[i] = v[i - 1] - 2v[i] + v[i + 1] end w[5] = v[4] - 2v[5] w .= w .* u nothing end -function Cfunc!(v,u,p,t) +function Cfunc!(v, u, p, t) w = zeros(5) - Cfunc!(w,v,u,p,t) + Cfunc!(w, v, u, p, t) w end @@ -227,11 +230,11 @@ adjoints, inverses, and more. For more information, see the [operator algebras t Great! You now know how to be state/parameter/time-dependent operators and make them matrix-free, along with doing algebras on operators. What's next? -* Interested in more examples of building operators? See the example of [making a fast fourier transform linear operator](@ref fft) -* Interested in more operators ready to go? See the [Premade Operators page](@ref premade_operators) for all of the operators included with SciMLOperators. Note that there are also downstream packages that make new operators. -* Want to make your own SciMLOperator? See the [AbstractSciMLOperator interface page](@ref operator_interface) which describes the full interface. + - Interested in more examples of building operators? See the example of [making a fast fourier transform linear operator](@ref fft) + - Interested in more operators ready to go? See the [Premade Operators page](@ref premade_operators) for all of the operators included with SciMLOperators. Note that there are also downstream packages that make new operators. + - Want to make your own SciMLOperator? See the [AbstractSciMLOperator interface page](@ref operator_interface) which describes the full interface. How do you use SciMLOperators? Check out the following downstream pages: -* [Using SciMLOperators in LinearSolve.jl for matrix-free Krylov methods](https://docs.sciml.ai/LinearSolve/stable/tutorials/linear/) -* [Using SciMLOperators in OrdinaryDiffEq.jl for semi-linear ODE solvers](https://docs.sciml.ai/DiffEqDocs/stable/solvers/nonautonomous_linear_ode/) + - [Using SciMLOperators in LinearSolve.jl for matrix-free Krylov methods](https://docs.sciml.ai/LinearSolve/stable/tutorials/linear/) + - [Using SciMLOperators in OrdinaryDiffEq.jl for semi-linear ODE solvers](https://docs.sciml.ai/DiffEqDocs/stable/solvers/nonautonomous_linear_ode/) diff --git a/docs/src/tutorials/operator_algebras.md b/docs/src/tutorials/operator_algebras.md index af29c904..c5762cdd 100644 --- a/docs/src/tutorials/operator_algebras.md +++ b/docs/src/tutorials/operator_algebras.md @@ -7,10 +7,10 @@ in order to build more complex objects and using their operations. ```@example operator_algebra using SciMLOperators, LinearAlgebra N = 4 -function f(v, u, p, t) +function f(v, u, p, t) u .* v end -function f(w, v, u, p, t) +function f(w, v, u, p, t) w .= u .* v end @@ -56,4 +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, α, β) -``` \ No newline at end of file +``` diff --git a/src/SciMLOperators.jl b/src/SciMLOperators.jl index 362d9f00..af5b0559 100644 --- a/src/SciMLOperators.jl +++ b/src/SciMLOperators.jl @@ -35,20 +35,20 @@ the following type of equation: w = L(u,p,t)[v] ``` -where `L[v]` is the operator application of ``L`` on the vector ``v``. +where `L[v]` is the operator application of ``L`` on the vector ``v``. ## Interface An `AbstractSciMLOperator` can be called like a function in the following ways: -- `L(v, u, p, t)` - Out-of-place application where `v` is the action vector and `u` is the update vector -- `L(w, v, u, p, t)` - In-place application where `w` is the destination, `v` is the action vector, and `u` is the update vector -- `L(w, v, u, p, t, α, β)` - In-place application with scaling: `w = α*(L*v) + β*w` + - `L(v, u, p, t)` - Out-of-place application where `v` is the action vector and `u` is the update vector + - `L(w, v, u, p, t)` - In-place application where `w` is the destination, `v` is the action vector, and `u` is the update vector + - `L(w, v, u, p, t, α, β)` - In-place application with scaling: `w = α*(L*v) + β*w` Operator state can be updated separately from application: -- `update_coefficients!(L, u, p, t)` for in-place operator update -- `L = update_coefficients(L, u, p, t)` for out-of-place operator update + - `update_coefficients!(L, u, p, t)` for in-place operator update + - `L = update_coefficients(L, u, p, t)` for out-of-place operator update SciMLOperators also overloads `Base.*`, `LinearAlgebra.mul!`, `LinearAlgebra.ldiv!` for operator evaluation without updating operator state. @@ -183,7 +183,6 @@ M = MatrixOperator(zero(N, N); update_func = mat_update_func, M(v, u, p, t) == zeros(N) # true M(v, u, p, t; scale = 1.0) != zero(N) ``` - """ abstract type AbstractSciMLOperator{T} end diff --git a/src/basic.jl b/src/basic.jl index 8236ed41..6689c9fd 100644 --- a/src/basic.jl +++ b/src/basic.jl @@ -77,14 +77,16 @@ function (ii::IdentityOperator)(v::AbstractVecOrMat, u, p, t; kwargs...) end # In-place: w is destination, v is action vector, u is update vector -function (ii::IdentityOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; kwargs...) +function (ii::IdentityOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; kwargs...) @assert size(v, 1) == ii.len update_coefficients!(ii, u, p, t; kwargs...) copy!(w, v) end # In-place with scaling: w = α*(ii*v) + β*w -function (ii::IdentityOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (ii::IdentityOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) @assert size(v, 1) == ii.len update_coefficients!(ii, u, p, t; kwargs...) mul!(w, I, v, α, β) @@ -181,7 +183,8 @@ function (nn::NullOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; k end # In-place with scaling: w = α*(nn*v) + β*w -function (nn::NullOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (nn::NullOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) @assert size(v, 1) == nn.len lmul!(β, w) w @@ -236,7 +239,8 @@ end # constructors for T in SCALINGNUMBERTYPES[2:end] - @eval ScaledOperator(λ::$T, L::AbstractSciMLOperator) = ScaledOperator( + @eval ScaledOperator( + λ::$T, L::AbstractSciMLOperator) = ScaledOperator( ScalarOperator(λ), L) end @@ -401,7 +405,8 @@ function (L::ScaledOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; end # In-place with scaling: w = α*(L*v) + β*w -function (L::ScaledOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (L::ScaledOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) update_coefficients!(L, u, p, t; kwargs...) if iszero(L.λ) lmul!(β, w) @@ -413,7 +418,6 @@ function (L::ScaledOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, end end - """ Lazy operator addition @@ -634,7 +638,8 @@ function (L::AddedOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; k end # In-place with scaling: w = α*(L*v) + β*w -function (L::AddedOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (L::AddedOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) update_coefficients!(L, u, p, t; kwargs...) lmul!(β, w) for op in L.ops @@ -680,15 +685,19 @@ end # constructors for op in (:*, :∘) @eval Base.$op(ops::AbstractSciMLOperator...) = reduce($op, ops) - @eval Base.$op(A::AbstractSciMLOperator, B::AbstractSciMLOperator) = ComposedOperator( + @eval Base.$op( + A::AbstractSciMLOperator, B::AbstractSciMLOperator) = ComposedOperator( A, B) - @eval Base.$op(A::ComposedOperator, B::AbstractSciMLOperator) = ComposedOperator( + @eval Base.$op( + A::ComposedOperator, B::AbstractSciMLOperator) = ComposedOperator( A.ops..., B) - @eval Base.$op(A::AbstractSciMLOperator, B::ComposedOperator) = ComposedOperator(A, + @eval Base.$op( + A::AbstractSciMLOperator, B::ComposedOperator) = ComposedOperator(A, B.ops...) - @eval Base.$op(A::ComposedOperator, B::ComposedOperator) = ComposedOperator(A.ops..., + @eval Base.$op( + A::ComposedOperator, B::ComposedOperator) = ComposedOperator(A.ops..., B.ops...) end @@ -786,7 +795,8 @@ for fact in (:lu, :lu!, :bunchkaufman, :bunchkaufman!, :lq, :lq!, :svd, :svd!) - @eval LinearAlgebra.$fact(L::ComposedOperator, args...) = prod( + @eval LinearAlgebra.$fact( + L::ComposedOperator, args...) = prod( op -> $fact(op, args...), reverse(L.ops)) end @@ -910,22 +920,23 @@ end function (L::ComposedOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; kwargs...) update_coefficients!(L, u, p, t; kwargs...) @assert iscached(L) "Cache needs to be set up for ComposedOperator. Call cache_operator(L, u) first." - - vecs = (w, L.cache[1:(end-1)]..., v) + + vecs = (w, L.cache[1:(end - 1)]..., v) for i in reverse(1:length(L.ops)) - L.ops[i](vecs[i], vecs[i+1], u, p, t; kwargs...) + L.ops[i](vecs[i], vecs[i + 1], u, p, t; kwargs...) end w end # In-place with scaling: w = α*(L*v) + β*w -function (L::ComposedOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (L::ComposedOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) update_coefficients!(L, u, p, t; kwargs...) @assert iscached(L) "Cache needs to be set up for ComposedOperator. Call cache_operator(L, u) first." - + cache = L.cache[end] copy!(cache, w) - + L(w, v, u, p, t; kwargs...) lmul!(α, w) axpy!(β, cache, w) @@ -1064,10 +1075,11 @@ function (L::InvertedOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t end # In-place with scaling: w = α*(L*v) + β*w -function (L::InvertedOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (L::InvertedOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) update_coefficients!(L, u, p, t; kwargs...) @assert iscached(L) "Cache needs to be set up for InvertedOperator. Call cache_operator(L, u) first." - + copy!(L.cache, w) ldiv!(w, L.L, v) lmul!(α, w) diff --git a/src/batch.jl b/src/batch.jl index bb2efbe5..9f037644 100644 --- a/src/batch.jl +++ b/src/batch.jl @@ -45,15 +45,20 @@ Base.iszero(L::BatchedDiagonalOperator) = iszero(L.diag) Base.transpose(L::BatchedDiagonalOperator) = L Base.adjoint(L::BatchedDiagonalOperator) = conj(L) function Base.conj(L::BatchedDiagonalOperator) # TODO - test this thoroughly - update_func, update_func! = if isreal(L) + update_func, + update_func! = if isreal(L) L.update_func, L.update_func! else - uf = L.update_func === nothing ? nothing : (L, u, p, t; kwargs...) -> conj(L.update_func(conj(L.diag), + uf = L.update_func === nothing ? nothing : + ( + L, u, p, t; kwargs...) -> conj(L.update_func(conj(L.diag), u, p, t; kwargs...)) - uf! = L.update_func! === nothing ? nothing : (L, u, p, t; kwargs...) -> begin + uf! = L.update_func! === nothing ? nothing : + ( + L, u, p, t; kwargs...) -> begin L.update_func!(conj!(L.diag), u, p, t; kwargs...) conj!(L.diag) end @@ -168,19 +173,21 @@ function (L::BatchedDiagonalOperator)(v::AbstractVecOrMat, u, p, t; kwargs...) L.diag .* v end -function (L::BatchedDiagonalOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; kwargs...) +function (L::BatchedDiagonalOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; kwargs...) update_coefficients!(L, u, p, t; kwargs...) w .= L.diag .* v return w end -function (L::BatchedDiagonalOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (L::BatchedDiagonalOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) update_coefficients!(L, u, p, t; kwargs...) if β == 0 w .= α .* (L.diag .* v) else w .= α .* (L.diag .* v) .+ β .* w - end + end return w end # diff --git a/src/func.jl b/src/func.jl index c9e7008b..ecccdaf6 100644 --- a/src/func.jl +++ b/src/func.jl @@ -45,7 +45,8 @@ function set_op( op) where {iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, iType, oType} return FunctionOperator{ iip, oop, mul5, T, typeof(op), Fa, Fi, Fai, Tr, U, P, Tt, C, iType, - oType}(op, f.op_adjoint, f.op_inverse, f.op_adjoint_inverse, f.traits, f.u, f.p, f.t, + oType}( + op, f.op_adjoint, f.op_inverse, f.op_adjoint_inverse, f.traits, f.u, f.p, f.t, f.cache) end @@ -54,7 +55,8 @@ function set_op_adjoint( iType, oType}, op_adjoint) where {iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, iType, oType} - return FunctionOperator{iip, oop, mul5, T, F, typeof(op_adjoint), Fi, Fai, Tr, U, P, Tt, + return FunctionOperator{ + iip, oop, mul5, T, F, typeof(op_adjoint), Fi, Fai, Tr, U, P, Tt, C, iType, oType}(f.op, op_adjoint, f.op_inverse, f.op_adjoint_inverse, f.traits, f.u, f.p, f.t, f.cache) end @@ -64,7 +66,8 @@ function set_op_inverse( iType, oType}, op_inverse) where {iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, iType, oType} - return FunctionOperator{iip, oop, mul5, T, F, Fa, typeof(op_inverse), Fai, Tr, U, P, Tt, + return FunctionOperator{ + iip, oop, mul5, T, F, Fa, typeof(op_inverse), Fai, Tr, U, P, Tt, C, iType, oType}(f.op, f.op_adjoint, op_inverse, f.op_adjoint_inverse, f.traits, f.u, f.p, f.t, f.cache) end @@ -94,8 +97,10 @@ function set_u( iType, oType}, u) where {iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, iType, oType} - return FunctionOperator{iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, typeof(u), P, Tt, C, iType, - oType}(f.op, f.op_adjoint, f.op_inverse, f.op_adjoint_inverse, f.traits, u, f.p, f.t, + return FunctionOperator{ + iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, typeof(u), P, Tt, C, iType, + oType}( + f.op, f.op_adjoint, f.op_inverse, f.op_adjoint_inverse, f.traits, u, f.p, f.t, f.cache) end @@ -104,8 +109,10 @@ function set_p( iType, oType}, p) where {iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, iType, oType} - return FunctionOperator{iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, typeof(p), Tt, C, iType, - oType}(f.op, f.op_adjoint, f.op_inverse, f.op_adjoint_inverse, f.traits, f.u, p, f.t, + return FunctionOperator{ + iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, typeof(p), Tt, C, iType, + oType}( + f.op, f.op_adjoint, f.op_inverse, f.op_adjoint_inverse, f.traits, f.u, p, f.t, f.cache) end @@ -113,8 +120,10 @@ function set_t( f::FunctionOperator{iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, iType, oType}, t) where {iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, iType, oType} - return FunctionOperator{iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, typeof(t), C, iType, - oType}(f.op, f.op_adjoint, f.op_inverse, f.op_adjoint_inverse, f.traits, f.u, f.p, t, + return FunctionOperator{ + iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, typeof(t), C, iType, + oType}( + f.op, f.op_adjoint, f.op_inverse, f.op_adjoint_inverse, f.traits, f.u, f.p, t, f.cache) end @@ -128,13 +137,15 @@ function set_cache( f.u, f.p, f.t, cache) end -function input_eltype(::FunctionOperator{iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, +function input_eltype(::FunctionOperator{ + iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, iType, oType }) where {iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, iType, oType} return iType end -function output_eltype(::FunctionOperator{iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, +function output_eltype(::FunctionOperator{ + iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, iType, oType }) where {iip, oop, mul5, T, F, Fa, Fi, Fai, Tr, U, P, Tt, C, iType, oType} return oType @@ -143,7 +154,7 @@ end """ $(SIGNATURES) -Wrap callable object `op` within an `AbstractSciMLOperator`. +Wrap callable object `op` within an `AbstractSciMLOperator`. ## Mathematical Description @@ -190,25 +201,25 @@ below traits. Keyword arguments are used to set operator traits, which are assumed to be uniform across `op`, `op_adjoint`, `op_inverse`, `op_adjoint_inverse`. -* `u` - Prototype of the state struct passed to the operator during evaluation, i.e. `L(u, p, t)`. `u` is set to `nothing` if no value is provided. -* `p` - Prototype of parameter struct passed to the operator during evaluation, i.e. `L(u, p, t)`. `p` is set to `nothing` if no value is provided. -* `t` - Protype of scalar time variable passed to the operator during evaluation. `t` is set to `zero(T)` if no value is provided. -* `accepted_kwargs` - `Tuple` of `Symbol`s corresponding to the keyword arguments accepted by `op*`, and `update_coefficients[!]`. For example, if `op` accepts kwarg `scale`, as in `op(u, p, t; scale)`, then `accepted_kwargs = (:scale,)`. -* `T` - `eltype` of the operator. If no value is provided, the constructor inferrs the value from types of `input`, and `output` -* `isinplace` - `true` if the operator can be used is a mutating way with in-place allocations. This trait is inferred if no value is provided. -* `outofplace` - `true` if the operator can be used is a non-mutating way with in-place allocations. This trait is inferred if no value is provided. -* `has_mul5` - `true` if the operator provides a five-argument `mul!` via the signature `op(v, u, p, t, α, β; )`. This trait is inferred if no value is provided. -* `isconstant` - `true` if the operator is constant, and doesn't need to be updated via `update_coefficients[!]` during operator evaluation. Defaults to false. -* `islinear` - `true` if the operator is linear. Defaults to `false`. -* `isconvertible` - `true` a cheap `convert(AbstractMatrix, L.op)` method is available. Defaults to `false`. -* `batch` - Boolean indicating if the input/output arrays comprise of batched column-vectors stacked in a matrix. If `true`, the input/output arrays must be `AbstractVecOrMat`s, and the length of the second dimension (the batch dimension) must be the same. The batch dimension is not involved in size computation. For example, with `batch = true`, and `size(output), size(input) = (M, K), (N, K)`, the `FunctionOperator` size is set to `(M, N)`. If `batch = false`, which is the default, the `input`/`output` arrays may of any size so long as `ndims(input) == ndims(output)`, and the `size` of `FunctionOperator` is set to `(length(input), length(output))`. -* `ifcache` - Allocate cache arrays in constructor. Defaults to `true`. Cache can be generated afterwards by calling `cache_operator(L, input, output)` -* `cache` - Pregenerated cache arrays for in-place evaluations. Expected to be of type and shape `(similar(input), similar(output),)`. The constructor generates cache if no values are provided. Cache generation by the constructor can be disabled by setting the kwarg `ifcache = false`. -* `opnorm` - The norm of `op`. Can be a `Number`, or function `opnorm(p::Integer)`. Defaults to `nothing`. -* `issymmetric` - `true` if the operator is linear and symmetric. Defaults to `false`. -* `ishermitian` - `true` if the operator is linear and hermitian. Defaults to `false`. -* `isposdef` - `true` if the operator is linear and positive-definite. Defaults to `false`. -* `kwargs` - Keyword arguments for cache initialization. If `accepted_kwargs` is provided, the corresponding keyword arguments must be passed. + - `u` - Prototype of the state struct passed to the operator during evaluation, i.e. `L(u, p, t)`. `u` is set to `nothing` if no value is provided. + - `p` - Prototype of parameter struct passed to the operator during evaluation, i.e. `L(u, p, t)`. `p` is set to `nothing` if no value is provided. + - `t` - Protype of scalar time variable passed to the operator during evaluation. `t` is set to `zero(T)` if no value is provided. + - `accepted_kwargs` - `Tuple` of `Symbol`s corresponding to the keyword arguments accepted by `op*`, and `update_coefficients[!]`. For example, if `op` accepts kwarg `scale`, as in `op(u, p, t; scale)`, then `accepted_kwargs = (:scale,)`. + - `T` - `eltype` of the operator. If no value is provided, the constructor inferrs the value from types of `input`, and `output` + - `isinplace` - `true` if the operator can be used is a mutating way with in-place allocations. This trait is inferred if no value is provided. + - `outofplace` - `true` if the operator can be used is a non-mutating way with in-place allocations. This trait is inferred if no value is provided. + - `has_mul5` - `true` if the operator provides a five-argument `mul!` via the signature `op(v, u, p, t, α, β; )`. This trait is inferred if no value is provided. + - `isconstant` - `true` if the operator is constant, and doesn't need to be updated via `update_coefficients[!]` during operator evaluation. Defaults to false. + - `islinear` - `true` if the operator is linear. Defaults to `false`. + - `isconvertible` - `true` a cheap `convert(AbstractMatrix, L.op)` method is available. Defaults to `false`. + - `batch` - Boolean indicating if the input/output arrays comprise of batched column-vectors stacked in a matrix. If `true`, the input/output arrays must be `AbstractVecOrMat`s, and the length of the second dimension (the batch dimension) must be the same. The batch dimension is not involved in size computation. For example, with `batch = true`, and `size(output), size(input) = (M, K), (N, K)`, the `FunctionOperator` size is set to `(M, N)`. If `batch = false`, which is the default, the `input`/`output` arrays may of any size so long as `ndims(input) == ndims(output)`, and the `size` of `FunctionOperator` is set to `(length(input), length(output))`. + - `ifcache` - Allocate cache arrays in constructor. Defaults to `true`. Cache can be generated afterwards by calling `cache_operator(L, input, output)` + - `cache` - Pregenerated cache arrays for in-place evaluations. Expected to be of type and shape `(similar(input), similar(output),)`. The constructor generates cache if no values are provided. Cache generation by the constructor can be disabled by setting the kwarg `ifcache = false`. + - `opnorm` - The norm of `op`. Can be a `Number`, or function `opnorm(p::Integer)`. Defaults to `nothing`. + - `issymmetric` - `true` if the operator is linear and symmetric. Defaults to `false`. + - `ishermitian` - `true` if the operator is linear and hermitian. Defaults to `false`. + - `isposdef` - `true` if the operator is linear and positive-definite. Defaults to `false`. + - `kwargs` - Keyword arguments for cache initialization. If `accepted_kwargs` is provided, the corresponding keyword arguments must be passed. """ function FunctionOperator(op, input::AbstractArray, @@ -849,7 +860,7 @@ function (L::FunctionOperator)(v::AbstractArray, u, p, t; kwargs...) L = update_coefficients(L, u, p, t; kwargs...) _sizecheck(L, v, nothing) V, _, vec_output = _unvec(L, v, nothing) - + # Apply the operator to action vector v after updating with u if L.traits.outofplace result = L.op(V, L.u, L.p, L.t; L.traits.kwargs...) @@ -867,11 +878,11 @@ end # In-place: w is destination, v is action vector, u is update vector function (L::FunctionOperator)(w::AbstractArray, v::AbstractArray, u, p, t; kwargs...) update_coefficients!(L, u, p, t; kwargs...) - + # Check dimensions _sizecheck(L, v, w) V, W, _ = _unvec(L, v, w) - + # Apply the operator in-place to action vector v after updating with u if L.traits.isinplace L.op(W, V, L.u, L.p, L.t; L.traits.kwargs...) @@ -880,18 +891,18 @@ function (L::FunctionOperator)(w::AbstractArray, v::AbstractArray, u, p, t; kwar result = L.op(V, L.u, L.p, L.t; L.traits.kwargs...) copyto!(W, result) end - + return w end # In-place with scaling: w = α*(L*v) + β*w function (L::FunctionOperator)(w::AbstractArray, v::AbstractArray, u, p, t, α, β; kwargs...) update_coefficients!(L, u, p, t; kwargs...) - + # Check dimensions _sizecheck(L, v, w) V, W, _ = _unvec(L, v, w) - + # Apply the operator in-place to action vector v with scaling if L.traits.isinplace && L.traits.has_mul5 # Direct 5-arg mul! if supported @@ -906,7 +917,7 @@ function (L::FunctionOperator)(w::AbstractArray, v::AbstractArray, u, p, t, α, result = L.op(V, L.u, L.p, L.t; L.traits.kwargs...) axpby!(β, W, α, result) end - + return w end # diff --git a/src/interface.jl b/src/interface.jl index 2aebd262..dfc7b839 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -63,7 +63,6 @@ result = L * v # Or use the interface which separrates the update from the application result = L(v, u, p, t; scale = 2.0) ``` - """ update_coefficients(L, u, p, t; kwargs...) = L @@ -98,7 +97,6 @@ t = 1.0 update_coefficients!(L, u, p, t) L * v ``` - """ update_coefficients!(L, u, p, t; kwargs...) = nothing @@ -157,10 +155,10 @@ computation. iscached(L) = true iscached(::Union{# LinearAlgebra -AbstractMatrix, -UniformScaling, -Factorization, # Base -Number}) = true + AbstractMatrix, + UniformScaling, + Factorization, # Base + Number}) = true """ $SIGNATURES @@ -253,10 +251,10 @@ Checks if an `L`'s state is constant or needs to be updated by calling `update_coefficients`. """ isconstant(::Union{# LinearAlgebra -AbstractMatrix, -UniformScaling, -Factorization, # Base -Number}) = true + AbstractMatrix, + UniformScaling, + Factorization, # Base + Number}) = true isconstant(L::AbstractSciMLOperator) = all(isconstant, getops(L)) isconstant(L) = false @@ -290,11 +288,13 @@ end Convert `SciMLOperator` to a concrete type via eager fusion. This method is a no-op for types that are already concrete. """ -concretize(L::Union{# LinearAlgebra -AbstractMatrix, -Factorization, # SciMLOperators -AbstractSciMLOperator -}) = convert(AbstractMatrix, L) +function concretize(L::Union{# LinearAlgebra + AbstractMatrix, + Factorization, # SciMLOperators + AbstractSciMLOperator +}) + convert(AbstractMatrix, L) +end function concretize(L::Union{ # LinearAlgebra @@ -318,48 +318,48 @@ islinear(::AbstractSciMLOperator) = false islinear(L) = false islinear(::Union{# LinearAlgebra -AbstractMatrix, -UniformScaling, -Factorization, # Base -Number + AbstractMatrix, + UniformScaling, + Factorization, # Base + Number }) = true has_mul(L) = false has_mul(::Union{# LinearAlgebra -AbstractVecOrMat, -AbstractMatrix, -UniformScaling, # Base -Number + AbstractVecOrMat, + AbstractMatrix, + UniformScaling, # Base + Number }) = true has_mul!(L) = false has_mul!(::Union{# LinearAlgebra -AbstractVecOrMat, -AbstractMatrix, -UniformScaling, # Base -Number + AbstractVecOrMat, + AbstractMatrix, + UniformScaling, # Base + Number }) = true has_ldiv(L) = false has_ldiv(::Union{ -AbstractMatrix, -Factorization, -Number + AbstractMatrix, + Factorization, + Number }) = true has_ldiv!(L) = false has_ldiv!(::Union{ -Diagonal, -Bidiagonal, -Factorization + Diagonal, + Bidiagonal, + Factorization }) = true has_adjoint(L) = islinear(L) has_adjoint(::Union{# LinearAlgebra -AbstractMatrix, -UniformScaling, -Factorization, # Base -Number + AbstractMatrix, + UniformScaling, + Factorization, # Base + Number }) = true """ @@ -368,9 +368,9 @@ Checks if `size(L, 1) == size(L, 2)`. issquare(L) = ndims(L) >= 2 && size(L, 1) == size(L, 2) issquare(::AbstractVector) = false issquare(::Union{# LinearAlgebra -UniformScaling, # SciMLOperators -AbstractSciMLScalarOperator, # Base -Number + UniformScaling, # SciMLOperators + AbstractSciMLScalarOperator, # Base + Number }) = true issquare(A...) = @. (&)(issquare(A)...) diff --git a/src/left.jl b/src/left.jl index 90055fa4..99719e0c 100644 --- a/src/left.jl +++ b/src/left.jl @@ -96,7 +96,8 @@ for (op, LType, VType) in ((:adjoint, :AdjointOperator, :AbstractAdjointVecOrMat # constructor @eval Base.$op(L::AbstractSciMLOperator) = $LType(L) - @eval Base.convert(::Type{AbstractMatrix}, L::$LType) = $op(convert(AbstractMatrix, + @eval Base.convert( + ::Type{AbstractMatrix}, L::$LType) = $op(convert(AbstractMatrix, L.L)) # traits @@ -158,7 +159,6 @@ for (op, LType, VType) in ((:adjoint, :AdjointOperator, :AbstractAdjointVecOrMat end end - # For AdjointOperator # Out-of-place: v is action vector, u is update vector function (L::AdjointOperator)(v::AbstractVecOrMat, u, p, t; kwargs...) @@ -181,9 +181,10 @@ function (L::AdjointOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; end # In-place with scaling: w = α*(L*v) + β*w -function (L::AdjointOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (L::AdjointOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) # Update the operator in-place - update_coefficients!(L.L, u, p ,t; kwargs...) + update_coefficients!(L.L, u, p, t; kwargs...) mul!(w', v', L.L, α, β) return w end @@ -191,20 +192,22 @@ end # For TransposedOperator # Out-of-place function (L::TransposedOperator)(v::AbstractVecOrMat, u, p, t; kwargs...) - L_updated = update_coefficients(L.L, u, p, t; kwargs...) - # (A^T)v = (v'A)' where v'A is computed by A'*v' - return (L_updated' * v')' + L_updated = update_coefficients(L.L, u, p, t; kwargs...) + # (A^T)v = (v'A)' where v'A is computed by A'*v' + return (L_updated' * v')' end # In-place -function (L::TransposedOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; kwargs...) +function (L::TransposedOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; kwargs...) update_coefficients!(L.L, u, p, t; kwargs...) mul!(w', v', L.L) return w end # In-place with scaling -function (L::TransposedOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (L::TransposedOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) update_coefficients!(L.L, u, p, t; kwargs...) mul!(w', v', L.L, α, β) return w diff --git a/src/matrix.jl b/src/matrix.jl index 4772e2c9..177af48a 100644 --- a/src/matrix.jl +++ b/src/matrix.jl @@ -9,6 +9,7 @@ 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; ) -> newA + or update_func!(A::AbstractMatrix, u, p, t; ) -> [modifies A] @@ -29,6 +30,7 @@ adjoints, transposes. # Example Out-of-place update and usage + ``` v = rand(4) u = rand(4) @@ -54,6 +56,7 @@ 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) @@ -127,12 +130,16 @@ for op in (:adjoint, @eval function Base.$op(L::MatrixOperator) isconstant(L) && return MatrixOperator($op(L.A)) - update_func = L.update_func === nothing ? nothing : (A, u, p, t; kwargs...) -> $op(L.update_func($op(L.A), + update_func = L.update_func === nothing ? nothing : + ( + A, u, p, t; kwargs...) -> $op(L.update_func($op(L.A), u, p, t; kwargs...)) - update_func! = L.update_func! === nothing ? nothing : (A, u, p, t; kwargs...) -> $op(L.update_func!($op(L.A), + update_func! = L.update_func! === nothing ? nothing : + ( + A, u, p, t; kwargs...) -> $op(L.update_func!($op(L.A), u, p, t; @@ -148,12 +155,16 @@ end function Base.conj(L::MatrixOperator) isconstant(L) && return MatrixOperator(conj(L.A)) - update_func = L.update_func === nothing ? nothing : (A, u, p, t; kwargs...) -> conj(L.update_func(conj(L.A), + update_func = L.update_func === nothing ? nothing : + ( + A, u, p, t; kwargs...) -> conj(L.update_func(conj(L.A), u, p, t; kwargs...)) - update_func! = L.update_func! === nothing ? nothing : (A, u, p, t; kwargs...) -> begin + update_func! = L.update_func! === nothing ? nothing : + ( + A, u, p, t; kwargs...) -> begin L.update_func!(conj!(L.A), u, p, t; kwargs...) conj!(L.A) end @@ -201,19 +212,21 @@ function (L::MatrixOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; end # In-place with scaling: w = α*(L*v) + β*w -function (L::MatrixOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (L::MatrixOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) update_coefficients!(L, u, p, t; kwargs...) mul!(w, L.A, v, α, β) end - # TODO - add tests for MatrixOperator indexing # propagate_inbounds here for the getindex fallback -Base.@propagate_inbounds Base.convert(::Type{AbstractMatrix}, L::MatrixOperator) = convert( +Base.@propagate_inbounds Base.convert( + ::Type{AbstractMatrix}, L::MatrixOperator) = convert( AbstractMatrix, L.A) Base.@propagate_inbounds Base.setindex!(L::MatrixOperator, v, i::Int) = (L.A[i] = v) -Base.@propagate_inbounds Base.setindex!(L::MatrixOperator, v, I::Vararg{Int, N}) where {N} = (L.A[I...] = v) +Base.@propagate_inbounds Base.setindex!( + L::MatrixOperator, v, I::Vararg{Int, N}) where {N} = (L.A[I...] = v) Base.eachcol(L::MatrixOperator) = eachcol(L.A) Base.eachrow(L::MatrixOperator) = eachrow(L.A) @@ -274,6 +287,7 @@ evaluation (`L([w,], v, u, p, t)`), or by calls to `update_function`, `update_func` which is assumed to have the signature update_func(diag::AbstractVecOrMat, u, p, t; ) -> new_diag + or update_func!(diag::AbstractVecOrMat, u, p, t; ) -> [modifies diag] @@ -286,15 +300,15 @@ are not provided. $(UPDATE_COEFFS_WARNING) # Example - """ function DiagonalOperator(diag::AbstractVector; update_func = nothing, update_func! = nothing, accepted_kwargs = nothing) diag_update_func = update_func_isconstant(update_func) ? update_func : - (A, u, p, t; kwargs...) -> update_func(A.diag, u, p, t; kwargs...) |> - Diagonal + ( + A, u, p, t; kwargs...) -> update_func(A.diag, u, p, t; kwargs...) |> + Diagonal diag_update_func! = update_func_isconstant(update_func!) ? update_func! : (A, u, p, t; kwargs...) -> update_func!(A.diag, u, p, t; kwargs...) @@ -348,9 +362,11 @@ for fact in (:lu, :lu!, :bunchkaufman, :bunchkaufman!, :lq, :lq!, :svd, :svd!) - @eval LinearAlgebra.$fact(L::AbstractSciMLOperator, args...) = InvertibleOperator(L, + @eval LinearAlgebra.$fact(L::AbstractSciMLOperator, + args...) = InvertibleOperator(L, $fact(convert(AbstractMatrix, L), args...)) - @eval LinearAlgebra.$fact(L::AbstractSciMLOperator; kwargs...) = InvertibleOperator(L, + @eval LinearAlgebra.$fact(L::AbstractSciMLOperator; + kwargs...) = InvertibleOperator(L, $fact(convert(AbstractMatrix, L); kwargs...)) end @@ -444,13 +460,15 @@ function (L::InvertibleOperator)(v::AbstractVecOrMat, u, p, t; kwargs...) end # In-place: w is destination, v is action vector, u is update vector -function (L::InvertibleOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; kwargs...) +function (L::InvertibleOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; kwargs...) update_coefficients!(L, u, p, t; kwargs...) mul!(w, L.L, v) end # In-place with scaling: w = α*(L*v) + β*w -function (L::InvertibleOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (L::InvertibleOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) update_coefficients!(L, u, p, t; kwargs...) mul!(w, L.L, v, α, β) end @@ -466,6 +484,7 @@ to `update_coefficients[!](L, u, p, t)`. The update functions are assumed to have the syntax update_func(b::AbstractVecOrMat, u, p, t; ) -> new_b + or update_func!(b::AbstractVecOrMat, u ,p , t; ) -> [modifies b] @@ -497,7 +516,6 @@ L = cache_operator(M, v) # update L and evaluate w = L(v, u, p, t) # == A * v + B * (p .* u * t) ``` - """ struct AffineOperator{T, AT, BT, bT, F, F!} <: AbstractSciMLOperator{T} A::AT @@ -705,7 +723,8 @@ function (L::AffineOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; end # In-place with scaling: w = α*(L*v) + β*w -function (L::AffineOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (L::AffineOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) update_coefficients!(L, u, p, t; kwargs...) # Scale the existing w by β lmul!(β, w) diff --git a/src/scalar.jl b/src/scalar.jl index b453e501..c5e4b788 100644 --- a/src/scalar.jl +++ b/src/scalar.jl @@ -174,7 +174,8 @@ ScalarOperator(λ::UniformScaling) = ScalarOperator(λ.λ) Base.show(io::IO, α::ScalarOperator) = print(io, "ScalarOperator($(α.val))") function Base.conj(α::ScalarOperator) # TODO - test val = conj(α.val) - update_func = (oldval, u, p, t; kwargs...) -> α.update_func(oldval |> conj, + update_func = (oldval, u, p, t; + kwargs...) -> α.update_func(oldval |> conj, u, p, t; @@ -205,7 +206,6 @@ function SciMLOperators.update_coefficients(L::ScalarOperator, u, p, t; kwargs.. return ScalarOperator(L.update_func(L.val, u, p, t; kwargs...), L.update_func) end - # Add ScalarOperator specific implementations for the new interface function (α::ScalarOperator)(v::AbstractArray, u, p, t; kwargs...) α = update_coefficients(α, u, p, t; kwargs...) @@ -222,7 +222,6 @@ function (α::ScalarOperator)(w::AbstractArray, v::AbstractArray, u, p, t, a, b; mul!(w, α, v, a, b) end - """ $TYPEDEF @@ -259,9 +258,11 @@ end for op in (:-, :+) for T in SCALINGNUMBERTYPES[2:end] - @eval Base.$op(α::AbstractSciMLScalarOperator, x::$T) = AddedScalarOperator(α, + @eval Base.$op(α::AbstractSciMLScalarOperator, + x::$T) = AddedScalarOperator(α, ScalarOperator($op(x))) - @eval Base.$op(x::$T, α::AbstractSciMLScalarOperator) = AddedScalarOperator( + @eval Base.$op(x::$T, + α::AbstractSciMLScalarOperator) = AddedScalarOperator( ScalarOperator(x), $op(α)) end @@ -307,12 +308,12 @@ function (α::AddedScalarOperator)(w::AbstractArray, v::AbstractArray, u, p, t; mul!(w, α, v) end -function (α::AddedScalarOperator)(w::AbstractArray, v::AbstractArray, u, p, t, a, b; kwargs...) +function (α::AddedScalarOperator)( + w::AbstractArray, v::AbstractArray, u, p, t, a, b; kwargs...) update_coefficients!(α, u, p, t; kwargs...) mul!(w, α, v, a, b) end - getops(α::AddedScalarOperator) = α.ops has_ldiv(α::AddedScalarOperator) = !iszero(convert(Number, α)) has_ldiv!(α::AddedScalarOperator) = has_ldiv(α) @@ -339,23 +340,29 @@ end for op in (:*, :∘) @eval Base.$op(ops::AbstractSciMLScalarOperator...) = reduce($op, ops) - @eval Base.$op(A::AbstractSciMLScalarOperator, B::AbstractSciMLScalarOperator) = ComposedScalarOperator( + @eval Base.$op(A::AbstractSciMLScalarOperator, + B::AbstractSciMLScalarOperator) = ComposedScalarOperator( A, B) - @eval Base.$op(A::ComposedScalarOperator, B::AbstractSciMLScalarOperator) = ComposedScalarOperator( + @eval Base.$op(A::ComposedScalarOperator, + B::AbstractSciMLScalarOperator) = ComposedScalarOperator( A.ops..., B) - @eval Base.$op(A::AbstractSciMLScalarOperator, B::ComposedScalarOperator) = ComposedScalarOperator( + @eval Base.$op(A::AbstractSciMLScalarOperator, + B::ComposedScalarOperator) = ComposedScalarOperator( A, B.ops...) - @eval Base.$op(A::ComposedScalarOperator, B::ComposedScalarOperator) = ComposedScalarOperator( + @eval Base.$op(A::ComposedScalarOperator, + B::ComposedScalarOperator) = ComposedScalarOperator( A.ops..., B.ops...) for T in SCALINGNUMBERTYPES[2:end] - @eval Base.$op(α::AbstractSciMLScalarOperator, x::$T) = ComposedScalarOperator(α, + @eval Base.$op(α::AbstractSciMLScalarOperator, + x::$T) = ComposedScalarOperator(α, ScalarOperator(x)) - @eval Base.$op(x::$T, α::AbstractSciMLScalarOperator) = ComposedScalarOperator( + @eval Base.$op(x::$T, + α::AbstractSciMLScalarOperator) = ComposedScalarOperator( ScalarOperator(x), α) end @@ -422,7 +429,8 @@ function (α::ComposedScalarOperator)(w::AbstractArray, v::AbstractArray, u, p, mul!(w, α, v) end -function (α::ComposedScalarOperator)(w::AbstractArray, v::AbstractArray, u, p, t, a, b; kwargs...) +function (α::ComposedScalarOperator)( + w::AbstractArray, v::AbstractArray, u, p, t, a, b; kwargs...) update_coefficients!(α, u, p, t; kwargs...) mul!(w, α, v, a, b) end @@ -435,7 +443,6 @@ has_ldiv!(α::ComposedScalarOperator) = all(has_ldiv!, α.ops) $TYPEDEF Lazy inverse of `AbstractSciMLScalarOperator`s - """ struct InvertedScalarOperator{T, λType} <: AbstractSciMLScalarOperator{T} λ::λType @@ -456,8 +463,9 @@ for op in (:/,) @eval Base.$op(x::$T, α::AbstractSciMLScalarOperator) = ScalarOperator(x) * inv(α) end - @eval Base.$op(α::AbstractSciMLScalarOperator, β::AbstractSciMLScalarOperator) = α * - inv(β) + @eval Base.$op( + α::AbstractSciMLScalarOperator, β::AbstractSciMLScalarOperator) = α * + inv(β) end for op in (:\,) @@ -466,8 +474,9 @@ for op in (:\,) @eval Base.$op(x::$T, α::AbstractSciMLScalarOperator) = inv(ScalarOperator(x)) * α end - @eval Base.$op(α::AbstractSciMLScalarOperator, β::AbstractSciMLScalarOperator) = inv(α) * - β + @eval Base.$op( + α::AbstractSciMLScalarOperator, β::AbstractSciMLScalarOperator) = inv(α) * + β end function Base.convert(T::Type{<:Number}, α::InvertedScalarOperator) @@ -497,7 +506,8 @@ function (α::InvertedScalarOperator)(w::AbstractArray, v::AbstractArray, u, p, mul!(w, α, v) end -function (α::InvertedScalarOperator)(w::AbstractArray, v::AbstractArray, u, p, t, a, b; kwargs...) +function (α::InvertedScalarOperator)( + w::AbstractArray, v::AbstractArray, u, p, t, a, b; kwargs...) update_coefficients!(α, u, p, t; kwargs...) mul!(w, α, v, a, b) end diff --git a/src/tensor.jl b/src/tensor.jl index 2f150135..6203f934 100644 --- a/src/tensor.jl +++ b/src/tensor.jl @@ -583,14 +583,16 @@ function (L::TensorProductOperator)(v::AbstractVecOrMat, u, p, t; kwargs...) end # In-place: w is destination, v is action vector, u is update vector -function (L::TensorProductOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; kwargs...) +function (L::TensorProductOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t; kwargs...) update_coefficients!(L, u, p, t; kwargs...) mul!(w, L, v) return w end # In-place with scaling: w = α*(L*v) + β*w -function (L::TensorProductOperator)(w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) +function (L::TensorProductOperator)( + w::AbstractVecOrMat, v::AbstractVecOrMat, u, p, t, α, β; kwargs...) update_coefficients!(L, u, p, t; kwargs...) mul!(w, L, v, α, β) return w diff --git a/src/utils.jl b/src/utils.jl index bcc5f9ec..632792a1 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -54,7 +54,7 @@ function (f::FilterKwargs)(args...; kwargs...) f.f(args...; filtered_kwargs...) end -isnothingfunc(f::FilterKwargs) = isnothingfunc(f.f) +isnothingfunc(f::FilterKwargs) = isnothingfunc(f.f) isnothingfunc(f::Nothing) = true isnothingfunc(f) = false # diff --git a/test/basic.jl b/test/basic.jl index 06dfd7bf..cf76ddd3 100644 --- a/test/basic.jl +++ b/test/basic.jl @@ -19,8 +19,8 @@ K = 12 @testset "IdentityOperator" begin A = rand(N, N) |> MatrixOperator u = rand(N, K) - v = rand(N, K) - w = zeros(N, K) + v = rand(N, K) + w = zeros(N, K) p = nothing t = 0 α = rand() @@ -47,15 +47,15 @@ K = 12 # Test with new interface - same update and action vector @test Id(u, u, p, t) ≈ u - + # Test with different vectors for update and action @test Id(v, u, p, t) ≈ v - + # Test in-place operation copy!(w, zeros(N, K)) Id(w, v, u, p, t) @test w ≈ v - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) @@ -82,8 +82,8 @@ end @testset "NullOperator" begin A = rand(N, N) |> MatrixOperator - u = rand(N, K) - v = rand(N, K) + u = rand(N, K) + v = rand(N, K) w = zeros(N, K) p = nothing t = 0 @@ -109,15 +109,15 @@ end # Test with new interface - same update and action vector @test Z(u, u, p, t) ≈ zero(u) - + # Test with different vectors for update and action @test Z(v, u, p, t) ≈ zero(v) - + # Test in-place operation copy!(w, ones(N, K)) Z(w, v, u, p, t) @test w ≈ zero(v) - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) @@ -167,15 +167,15 @@ end # Test with new interface - same vector for update and action @test op(u, u, p, t) ≈ α * A * u - + # Test with different vectors for update and action @test op(v, u, p, t) ≈ α * A * v - + # Test in-place operation copy!(w, zeros(N, K)) op(w, v, u, p, t) @test w ≈ α * A * v - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) @@ -235,24 +235,24 @@ end @test op2 * u ≈ op(α * A * u, B * u) @test op3 * u ≈ op(A * u, β * B * u) @test op4 * u ≈ op(α * A * u, β * B * u) - + # Test new interface - combined case @test op1(u, u, p, t) ≈ op(A * u, B * u) @test op2(u, u, p, t) ≈ op(α * A * u, B * u) @test op3(u, u, p, t) ≈ op(A * u, β * B * u) @test op4(u, u, p, t) ≈ op(α * A * u, β * B * u) - + # Test new interface - separate vectors @test op1(v, u, p, t) ≈ op(A * v, B * v) @test op2(v, u, p, t) ≈ op(α * A * v, B * v) @test op3(v, u, p, t) ≈ op(A * v, β * B * v) @test op4(v, u, p, t) ≈ op(α * A * v, β * B * v) - + # Test in-place operation copy!(w, zeros(N, K)) op1(w, v, u, p, t) @test w ≈ op(A * v, B * v) - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) @@ -301,7 +301,7 @@ end H_dense = c1 * A1_dense + c2 * A2_dense + c3 * A3_dense u = rand(T, N) - v = rand(T, N) + v = rand(T, N) du = similar(u) p = (ω = 0.1,) t = 0.1 @@ -341,22 +341,22 @@ end @test ABCmulu ≈ op * u @test ABCdivu ≈ op \ u ≈ opF \ u - + # Test new interface - combined case @test op(u, u, p, t) ≈ ABCmulu - + # Test new interface - separate vectors @test op(v, u, p, t) ≈ (A * B * C) * v @test !iscached(op) op = cache_operator(op, u) @test iscached(op) - + # Test in-place operation with new interface copy!(w, zeros(N, K)) op(w, v, u, p, t) @test w ≈ (A * B * C) * v - + # Test in-place with scaling with new interface copy!(w, rand(N, K)) orig_w = copy(w) @@ -404,8 +404,8 @@ end @testset "Adjoint, Transpose" begin for (op, - LType, - VType) in ((adjoint, AdjointOperator, AbstractAdjointVecOrMat), + LType, + VType) in ((adjoint, AdjointOperator, AbstractAdjointVecOrMat), (transpose, TransposedOperator, AbstractTransposedVecOrMat)) A = rand(N, N) D = Bidiagonal(rand(N, N), :L) @@ -480,18 +480,18 @@ end @test islinear(Di) @test Di * u ≈ u ./ s - + # Test new interface - same vectors @test Di(u, u, p, t) ≈ u ./ s - + # Test new interface - separate vectors @test Di(v, u, p, t) ≈ v ./ s - + # Test in-place operation copy!(w, zeros(N)) Di(w, v, u, p, t) @test w ≈ v ./ s - + # Test in-place with scaling copy!(w, rand(N)) orig_w = copy(w) @@ -510,4 +510,4 @@ end @test ldiv!(v, Di, u) ≈ u .* s v = copy(u) @test ldiv!(Di, u) ≈ v .* s -end \ No newline at end of file +end diff --git a/test/downstream/alloccheck.jl b/test/downstream/alloccheck.jl index 5a6b125c..47f9b062 100644 --- a/test/downstream/alloccheck.jl +++ b/test/downstream/alloccheck.jl @@ -52,7 +52,7 @@ for T in (Float32, Float64, ComplexF32, ComplexF64) H_dense = c1 * A1_dense + c2 * A2_dense + c3 * A3_dense u = rand(T, N) - v = rand(T, N) + v = rand(T, N) w = similar(u) p = (ω = 0.1,) t = 0.1 diff --git a/test/func.jl b/test/func.jl index 807e078b..c672afea 100644 --- a/test/func.jl +++ b/test/func.jl @@ -50,17 +50,18 @@ NK = N * K # test with ND-arrays and new interface @test _mul(A, v) ≈ L(v, u, p, t) ≈ L * v ≈ mul!(zero(w), L, v) @test α * _mul(A, v) + β * w ≈ mul!(copy(w), L, v, α, β) - + # Test with different update and action vectors action_vec = rand(sz_in...) @test _mul(A, action_vec) ≈ L(action_vec, u, p, t) - + if sz_in == sz_out @test _div(A, w) ≈ L \ w ≈ ldiv!(zero(v), L, w) ≈ ldiv!(L, copy(w)) end - + # test with vec(Array) - @test vec(_mul(A, v)) ≈ L(vec(v), u, p, t) ≈ L * vec(v) ≈ mul!(vec(zero(w)), L, vec(v)) + @test vec(_mul(A, v)) ≈ L(vec(v), u, p, t) ≈ L * vec(v) ≈ + mul!(vec(zero(w)), L, vec(v)) @test vec(α * _mul(A, v) + β * w) ≈ mul!(vec(copy(w)), L, vec(v), α, β) if sz_in == sz_out @@ -72,7 +73,7 @@ NK = N * K output_vec = zeros(sz_out...) L(output_vec, action_vec, u, p, t) @test output_vec ≈ _mul(A, action_vec) - + # Test in-place with scaling output_vec = rand(sz_out...) orig_output = copy(output_vec) @@ -128,7 +129,7 @@ end @test islinear(op1) @test islinear(op2) @test op1' === op1 - + # Test operator properties @test size(op1) == (NK, NK) @test has_adjoint(op1) @@ -181,11 +182,11 @@ end # Test with new interface - out of place @test _mul(A, action_vec) ≈ op1(action_vec, u, p, t) - + # Test with new interface - in place op2(result_vec, action_vec, u, p, t) @test result_vec ≈ _mul(A, action_vec) - + # Test in-place with scaling result_vec = rand(N, K) orig_result = copy(result_vec) @@ -278,11 +279,11 @@ end v = rand(N, K) @test *(A, v) ≈ op1 * v ≈ mul!(w, op2, v) - + # Test with new interface v = rand(N, K) @test *(A, v) ≈ op1(w, v, u, p, t) ≈ op2(w, v, u, p, t) - + v = rand(N, K) w = copy(v) @test α * *(A, v) + β * w ≈ mul!(w, op2, v, α, β) @@ -292,7 +293,7 @@ end @test \(A, w) ≈ op1 \ w ≈ ldiv!(v, op2, w) w = copy(v) @test \(A, w) ≈ ldiv!(op2, w) - + # Test new interface ldiv w = rand(N, K) ldiv_result = zeros(N, K) @@ -316,12 +317,12 @@ end for acc_kw in ((:scale,), Val((:scale,))) # Function operator with keyword arguments L = FunctionOperator(f, v, w; - u = u, - p = zero(p), - t = zero(t), - batch = true, - accepted_kwargs = acc_kw, - scale = 1.0) + u = u, + p = zero(p), + t = zero(t), + batch = true, + accepted_kwargs = acc_kw, + scale = 1.0) @test_throws ArgumentError FunctionOperator( f, v, w; u = u, p = zero(p), t = zero(t), batch = true, @@ -333,15 +334,15 @@ end A = Diagonal(u * p * t * scale) expected = A * v ans = u * p .* t .* scale - + # Test with new interface @test L(v, u, p, t; scale) ≈ expected - + # Test in-place with new interface copy!(w, zeros(N, K)) - L(w, v, u, p, t; scale) + L(w, v, u, p, t; scale) @test w ≈ expected - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) @@ -349,7 +350,7 @@ end β_val = rand() L(w, v, u, p, t, α_val, β_val; scale) @test w ≈ α_val * expected + β_val * orig_w - + # Test that outputs aren't accidentally mutated v1 = rand(N, K) v2 = rand(N, K) @@ -359,17 +360,17 @@ end # Expected results with different vectors result1 = A * v1 result2 = A * v2 - + # Test output consistency w1 = zeros(N, K) w2 = zeros(N, K) - + L(w1, v1, u, p, t; scale) L(w2, v2, u, p, t; scale) - + @test w1 ≈ result1 @test w2 ≈ result2 - + # Test matrix-vector multiplication w1 = L * v1 @test w1 ≈ A * v1 @@ -381,21 +382,21 @@ end # Test in-place matrix-vector multiplication v1 .= 0.0 v2 .= 0.0 - + mul!(w1, L, v1) @test w1 ≈ A * v1 mul!(w2, L, v2) @test w2 ≈ A * v2 @test w1 ≈ A * v1 @test w1 + w2 ≈ A * (v1 + v2) - + # Test scaling v1 = rand(N, K) w1 = copy(v1) v2 = rand(N, K) w2 = copy(v2) a1, a2, b1, b2 = rand(4) - + res = copy(w1) mul!(res, L, v1, a1, b1) @test res ≈ a1 * A * v1 + b1 * w1 @@ -405,4 +406,4 @@ end @test res ≈ a1 * A * v1 + b1 * w1 @test res + res2 ≈ (a1 * A * v1 + b1 * w1) + (a2 * A * v2 + b2 * w2) end -end \ No newline at end of file +end diff --git a/test/matrix.jl b/test/matrix.jl index 64170eb4..4b912d61 100644 --- a/test/matrix.jl +++ b/test/matrix.jl @@ -14,7 +14,7 @@ K = 19 u = rand(N, K) # Both update and action vector v = rand(N, K) # Output/action vector w = zeros(N, K) # Output vector - + p = nothing t = 0 α = rand() @@ -58,7 +58,7 @@ K = 19 # Test with new interface - same vector for update and action @test A * u ≈ AA(u, u, p, t) @test At * u ≈ AAt(u, u, p, t) - + # Test with different vectors for update and action @test A * v ≈ AA(v, u, p, t) @test At * v ≈ AAt(v, u, p, t) @@ -70,7 +70,7 @@ K = 19 copy!(w, zeros(N, K)) AA(w, v, u, p, t) @test w ≈ A * v - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) @@ -83,7 +83,7 @@ end u = rand(N, K) # Update vector v = rand(N, K) # Action vector w = zeros(N, K) # Output vector - + p = nothing t = 0 α = rand() @@ -106,7 +106,7 @@ end copy!(w, zeros(N, K)) L(w, v, u, p, t) @test w ≈ d .* v - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) @@ -117,7 +117,7 @@ end copy!(w, zeros(N, K)) ldiv!(w, L, u) @test w ≈ d .\ u - + # Existing test for in-place ldiv! v = copy(u) ldiv!(L, v) @@ -129,7 +129,7 @@ end u = rand(N, K) # Update vector v = rand(N, K) # Action vector w = zeros(N, K) # Output vector - + p = rand(N) t = rand() α = rand() @@ -143,30 +143,30 @@ end # Expected matrix after update A = p * p' - + # Test with new interface - same vector for update and action @test L(u, u, p, t) ≈ A * u - + # Test with different vectors for update and action @test L(v, u, p, t) ≈ A * v - + # Test in-place operation copy!(w, zeros(N, K)) L(w, v, u, p, t) @test w ≈ A * v - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) L(w, v, u, p, t, α, β) @test w ≈ α * (A * v) + β * orig_w - A = [-2.0 1 0 0 0 - 1 -2 1 0 0 - 0 1 -2 1 0 - 0 0 1 -2 1 - 0 0 0 1 -2] - v = [3.0,2.0,1.0,2.0,3.0] + A = [-2.0 1 0 0 0 + 1 -2 1 0 0 + 0 1 -2 1 0 + 0 0 1 -2 1 + 0 0 0 1 -2] + v = [3.0, 2.0, 1.0, 2.0, 3.0] opA = MatrixOperator(A) function update_function!(B, u, p, t) @@ -174,28 +174,31 @@ end B .= A .* u + dt*I end - u = Array(1:1.0:5); p = 0.1; t = 0.0 + u = Array(1:1.0:5); + p = 0.1; + t = 0.0 opB = MatrixOperator(copy(A); update_func! = update_function!) - function Bfunc!(w,v,u,p,t) + function Bfunc!(w, v, u, p, t) dt = p w[1] = -(2*u[1]-dt)*v[1] + v[2]*u[1] for i in 2:4 - w[i] = v[i-1]*u[i] - (2*u[i]-dt)*v[i] + v[i+1]*u[i] + w[i] = v[i - 1]*u[i] - (2*u[i]-dt)*v[i] + v[i + 1]*u[i] end w[5] = v[4]*u[5] - (2*u[5]-dt)*v[5] nothing end - function Bfunc!(v,u,p,t) + function Bfunc!(v, u, p, t) w = zeros(5) - Bfunc!(w,v,u,p,t) + Bfunc!(w, v, u, p, t) w end - mfopB = FunctionOperator(Bfunc!, zeros(5), zeros(5); u, p, t, isconstant=false) + mfopB = FunctionOperator(Bfunc!, zeros(5), zeros(5); u, p, t, isconstant = false) - @test iszero(opB(v, Array(2:1.0:6), 0.5, nothing) - mfopB(v, Array(2:1.0:6), 0.5, nothing)) + @test iszero(opB(v, Array(2:1.0:6), 0.5, nothing) - + mfopB(v, Array(2:1.0:6), 0.5, nothing)) end @testset "DiagonalOperator update test" begin @@ -203,7 +206,7 @@ end u = rand(N, K) # Update vector v = rand(N, K) # Action vector w = zeros(N, K) # Output vector - + p = rand(N) t = rand() α = rand() @@ -219,15 +222,15 @@ end # Expected result after update expected = (p * t) .* v - + # Test with new interface - different vectors for update and action @test D(v, u, p, t) ≈ expected - + # Test in-place operation copy!(w, zeros(N, K)) D(w, v, u, p, t) @test w ≈ expected - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) @@ -240,7 +243,7 @@ end u = rand(N, K) # Update vector v = rand(N, K) # Action vector w = zeros(N, K) # Output vector - + d = rand(N, K) α = rand() β = rand() @@ -256,12 +259,12 @@ end # Test with new interface @test L(v, u, p, t) ≈ d .* v - + # Test in-place operation copy!(w, zeros(N, K)) L(w, v, u, p, t) @test w ≈ d .* v - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) @@ -270,11 +273,11 @@ end # Test division operations @test L \ u ≈ d .\ u - + copy!(w, zeros(N, K)) ldiv!(w, L, u) @test w ≈ d .\ u - + # Existing test for in-place ldiv! v = copy(u) ldiv!(L, u) @@ -286,7 +289,7 @@ end u = rand(N, K) # Update vector v = rand(N, K) # Action vector w = zeros(N, K) # Output vector - + d = zeros(N, K) p = rand(N, K) t = rand() @@ -301,10 +304,10 @@ end # Expected result after update expected = (p * t) .* v - + # Test with new interface - different vectors for update and action @test D(v, u, p, t) ≈ expected - + # Test in-place operation copy!(w, zeros(N, K)) D(w, v, u, p, t) @@ -316,7 +319,7 @@ end u = rand(N, K) # Update vector v = rand(N, K) # Action vector w = zeros(N, K) # Output vector - + A = rand(N, N) B = rand(N, N) D = Diagonal(A) @@ -334,12 +337,12 @@ end # Test with new interface expected = A * v + B * b @test L(v, u, p, t) ≈ expected - + # Test in-place operation copy!(w, zeros(N, K)) L(w, v, u, p, t) @test w ≈ expected - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) @@ -352,11 +355,11 @@ end # Test division operations @test L \ u ≈ D \ (u - B * b) - + copy!(w, zeros(N, K)) ldiv!(w, L, u) @test w ≈ D \ (u - B * b) - + # Existing test for in-place ldiv! v = copy(u) ldiv!(L, u) @@ -370,23 +373,23 @@ end expected = v + b @test L(v, u, p, t) ≈ expected @test L \ u ≈ u - b - + # Test in-place operation copy!(w, zeros(N, K)) L(w, v, u, p, t) @test w ≈ expected - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) L(w, v, u, p, t, α, β) @test w ≈ α * expected + β * orig_w - + # Test division operations copy!(w, zeros(N, K)) ldiv!(w, L, u) @test w ≈ u - b - + # Existing test for in-place ldiv! v = copy(u) ldiv!(L, u) @@ -400,23 +403,23 @@ end expected = v + B * b @test L(v, u, p, t) ≈ expected @test L \ u ≈ u - B * b - + # Test in-place operation copy!(w, zeros(N, K)) L(w, v, u, p, t) @test w ≈ expected - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) L(w, v, u, p, t, α, β) @test w ≈ α * expected + β * orig_w - + # Test division operations copy!(w, zeros(N, K)) ldiv!(w, L, u) @test w ≈ u - B * b - + # Existing test for in-place ldiv! v = copy(u) ldiv!(L, u) @@ -428,7 +431,7 @@ end u = rand(N, K) # Update vector v = rand(N, K) # Action vector w = zeros(N, K) # Output vector - + A = rand(N, N) B = rand(N, N) p = rand(N, K) @@ -445,15 +448,15 @@ end # Expected updated bias and result b = p * t expected = A * v + B * b - + # Test with new interface - different vectors for update and action @test L(v, u, p, t) ≈ expected - + # Test in-place operation copy!(w, zeros(N, K)) L(w, v, u, p, t) @test w ≈ expected - + # Test in-place with scaling copy!(w, rand(N, K)) orig_w = copy(w) @@ -496,11 +499,11 @@ end # Inputs/Update vectors u2 = rand(n1 * n2, K) u3 = rand(n1 * n2 * n3, K) - + # Action vectors (same as update vectors initially) v2 = copy(u2) v3 = copy(u3) - + # Output vectors w2 = zeros(m1 * m2, K) w3 = zeros(m1 * m2 * m3, K) @@ -545,11 +548,11 @@ end @test AB * v2 ≈ opAB(v2, u2, p, t) @test ABC * v3 ≈ opABC(v3, u3, p, t) - @test AB \ w2 ≈ opAB \ w2 + @test AB \ w2 ≈ opAB \ w2 @test AB \ w2 ≈ opAB_F \ w2 - @test ABC \ w3 ≈ opABC \ w3 + @test ABC \ w3 ≈ opABC \ w3 @test ABC \ w3 ≈ opABC_F \ w3 - + @test !iscached(opAB) @test !iscached(opABC) @@ -578,7 +581,7 @@ end w2 = zeros(M2, K) # Output vector opAB(w2, v2, u2, p, t) @test w2 ≈ AB * v2 - + v3 = rand(n1 * n2 * n3, K) # Action vector w3 = zeros(M3, K) # Output vector opABC(w3, v3, u3, p, t) @@ -590,7 +593,7 @@ end orig_w2 = copy(w2) opAB(w2, v2, u2, p, t, α, β) @test w2 ≈ α * AB * v2 + β * orig_w2 - + v3 = rand(n1 * n2 * n3, K) # Action vector w3 = rand(M3, K) # Output vector orig_w3 = copy(w3) @@ -602,23 +605,23 @@ end v2 = rand(M2, K) # Action vector (size of output space) u2 = rand(N2, K) # Update vector (size of input space) w2 = zeros(N2, K) # Output vector (size of input space) - + # ldiv! with new interface ldiv!(w2, opAB_F, v2) @test w2 ≈ AB \ v2 - + v3 = rand(M3, K) # Action vector u3 = rand(N3, K) # Update vector w3 = zeros(N3, K) # Output vector ldiv!(w3, opABC_F, v3) @test w3 ≈ ABC \ v3 - + # In-place ldiv! (original style) v2 = rand(M2, K) u2 = copy(v2) ldiv!(opAB_F, v2) @test v2 ≈ AB \ u2 - + v3 = rand(M3, K) u3 = copy(v3) ldiv!(opABC_F, v3) @@ -627,16 +630,16 @@ end v2 = rand(M2, K) # Action vector u2 = rand(N2, K) # Update vector w2 = zeros(N2, K) # Output vector - + if VERSION < v"1.9-" @test_broken ldiv!(w2, opAB_F, v2) ≈ AB \ v2 else @test ldiv!(w2, opAB_F, v2) ≈ AB \ v2 end - + v3 = rand(M3, K) # Action vector u3 = rand(N3, K) # Update vector w3 = zeros(N3, K) # Output vector @test_broken ldiv!(w3, opABC_F, v3) ≈ ABC \ v3 # errors end -end \ No newline at end of file +end diff --git a/test/scalar.jl b/test/scalar.jl index 1bd64066..a0a7b06b 100644 --- a/test/scalar.jl +++ b/test/scalar.jl @@ -60,15 +60,15 @@ K = 12 # Tests with the new interface v = copy(u) # Action vector w = zeros(N, K) # Output vector - + # Test with new interface result = α(v, u, nothing, 0.0) @test result ≈ v * x - + # Test in-place operations α(w, v, u, nothing, 0.0) @test w ≈ v * x - + # Test in-place operations with scaling orig_w = rand(N, K) copy!(w, orig_w) @@ -81,7 +81,7 @@ end α = ScalarOperator(x) u = rand(N, K) # Update vector v = rand(N, K) # Action vector - + # Test scalar operator combinations β = α + α @test β isa AddedScalarOperator @@ -103,7 +103,7 @@ end @test β * u ≈ (1 / x) * u # Original style test @inferred convert(Float32, β) @test convert(Number, β) ≈ 1 / x - + β = α * inv(α) @test β isa ComposedScalarOperator @test β * u ≈ u @@ -115,30 +115,30 @@ end @test β * u ≈ u @inferred convert(Float32, β) @test convert(Number, β) ≈ true - + # Test combination with other operators for op in (MatrixOperator(rand(N, N)), SciMLOperators.IdentityOperator(N)) @test α + op isa SciMLOperators.AddedOperator @test (α + op) * u ≈ x * u + op * u - + L = α + op @test L isa SciMLOperators.AddedOperator @test L(v, u, nothing, 0.0) ≈ x * v + op * v - + @test α * op isa SciMLOperators.ScaledOperator @test (α * op) * u ≈ x * (op * u) - + L = α * op @test L isa SciMLOperators.ScaledOperator @test L(v, u, nothing, 0.0) ≈ x * (op * v) - + # Division tests from original @test all(map(T -> (T isa SciMLOperators.ScaledOperator), (α / op, op / α, op \ α, α \ op))) @test (α / op) * u ≈ (op \ α) * u ≈ α * (op \ u) @test (op / α) * u ≈ (α \ op) * u ≈ 1 / α * op * u end - + # Test for ComposedScalarOperator nesting (from original) α_new = ScalarOperator(rand()) L = α_new * (α_new * α_new) * α_new @@ -167,7 +167,7 @@ end w = zeros(N, K) # Output vector p = 2.0 t = 4.0 - + c = rand() d = rand() @@ -197,7 +197,7 @@ end # Tests with new interface @test α(v, u, p, t) ≈ p * v @test β(v, u, p, t) ≈ t * v - + # Test in-place with scaling orig_w = rand(N, K) copy!(w, orig_w) @@ -214,40 +214,40 @@ end # Test operator combinations num = α + 2 / β * 3 - 4 val = p + 2 / t * 3 - 4 - + @test convert(Number, num) ≈ val # Test with keyword arguments γ = ScalarOperator(0.0; update_func = (args...; dtgamma) -> dtgamma, - accepted_kwargs = (:dtgamma,)) - + accepted_kwargs = (:dtgamma,)) + dtgamma = rand() # Original tests @test_throws MethodError γ(u, p, t; dtgamma) ≈ dtgamma * u - + # New interface tests @test γ(v, u, p, t; dtgamma) ≈ dtgamma * v - + # In-place test with keywords w_test = zeros(N, K) γ(w_test, v, u, p, t; dtgamma) @test w_test ≈ dtgamma * v - + γ_added = γ + α # Original tests @test_throws MethodError γ_added(u, p, t; dtgamma) ≈ (dtgamma + p) * u - + # New interface tests @test γ_added(v, u, p, t; dtgamma) ≈ (dtgamma + p) * v - + # In-place test with keywords for combined operator w_test = zeros(N, K) γ_added(w_test, v, u, p, t; dtgamma) @test w_test ≈ (dtgamma + p) * v - + # In-place test with scaling and keywords w_test = rand(N, K) w_orig = copy(w_test) γ_added(w_test, v, u, p, t, c, d; dtgamma) @test w_test ≈ c * (dtgamma + p) * v + d * w_orig -end \ No newline at end of file +end diff --git a/test/total.jl b/test/total.jl index 6033f5f9..ca513458 100644 --- a/test/total.jl +++ b/test/total.jl @@ -126,26 +126,26 @@ end # Update operator @test_nowarn update_coefficients!(op, u, p, t; diag, matrix) - + # Form dense operator manually dense_T1 = kron(A, p * ones(N, N)) dense_T2 = kron(_C, (p * t) .* matrix) dense_DD = Diagonal(vcat(p * ones(N2), p * t * diag)) dense_op = hcat(dense_T1', dense_T2') * dense_DD * vcat(dense_T1, dense_T2) - + # Test correctness of op @test op * v ≈ dense_op * v - + # Test consistency with three-arg mul! w = rand(N2, K) @test mul!(w, op, v) ≈ op * v - + # Test consistency with in-place five-arg mul! w = rand(N2, K) w2 = copy(w) @test mul!(w, op, v, α, β) ≈ α * (op * v) + β * w2 - - # Create a fresh operator for each test + + # Create a fresh operator for each test op_fresh = TT' * DD * TT op_fresh = cache_operator(op_fresh, v) # Use in-place update directly in test @@ -155,7 +155,6 @@ end @test result1 ≈ dense_op * v end - @testset "Resize! test" begin M1 = 4 M2 = 12 diff --git a/test/zygote.jl b/test/zygote.jl index b5d3defa..8121dec6 100644 --- a/test/zygote.jl +++ b/test/zygote.jl @@ -66,7 +66,7 @@ for (LType, L) in ((IdentityOperator, IdentityOperator(N)), (AddedScalarOperator, α + α), (ComposedScalarOperator, α * α)) @assert L isa LType - + # Cache the operator for efficient application L_cached = cache_operator(L, u0) @@ -81,11 +81,11 @@ for (LType, L) in ((IdentityOperator, IdentityOperator(N)), loss_div = function (p) v = Diagonal(p) * u0 - + # Update coefficients first, then apply inverse L_updated = update_coefficients(L_cached, u0, p, t) w = L_updated \ v - + l = sum(w) end @@ -106,4 +106,4 @@ for (LType, L) in ((IdentityOperator, IdentityOperator(N)), @test !isa(g_div, Nothing) end end -end \ No newline at end of file +end