Merge pull request #363 from SebastianM-C/move_scale

Move out `DiffEqScaledOperator`

Author: ChrisRackauckas

Project.toml

@@ -14,6 +14,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

src/DiffEqOperators.jl

@@ -3,15 +3,16 @@ module DiffEqOperators
 import Base: +, -, *, /, \, size, getindex, setindex!, Matrix, convert, ==
 using DiffEqBase, StaticArrays, LinearAlgebra
 import LinearAlgebra: mul!, ldiv!, lmul!, rmul!, axpy!, opnorm, factorize, I
-import DiffEqBase: AbstractDiffEqLinearOperator, update_coefficients!, isconstant
+import DiffEqBase: update_coefficients!, isconstant
+using SciMLBase: AbstractDiffEqLinearOperator, AbstractDiffEqCompositeOperator, DiffEqScaledOperator
+import SciMLBase: getops
 using SparseArrays, ForwardDiff, BandedMatrices, NNlib, LazyArrays, BlockBandedMatrices
 using LazyBandedMatrices, ModelingToolkit
 using RuntimeGeneratedFunctions
 RuntimeGeneratedFunctions.init(@__MODULE__)
 
 abstract type AbstractDiffEqAffineOperator{T} end
 abstract type AbstractDerivativeOperator{T} <: AbstractDiffEqLinearOperator{T} end
-abstract type AbstractDiffEqCompositeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 abstract type AbstractMatrixFreeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 
 ### Matrix-free Operators

src/composite_operators.jl

@@ -2,52 +2,6 @@
 # derivative) using arithmetic or other operator compositions. The composite
 # operator types are lazy and maintain the structure used to build them.
 
-# Recursive routines that use `getops`
-function update_coefficients!(L::AbstractDiffEqCompositeOperator,u,p,t)
-  for op in getops(L)
-    update_coefficients!(op,u,p,t)
-  end
-  L
-end
-isconstant(L::AbstractDiffEqCompositeOperator) = all(isconstant, getops(L))
-
-# Scaled operator (α * A)
-struct DiffEqScaledOperator{T,F,OpType<:AbstractDiffEqLinearOperator{T}} <: AbstractDiffEqCompositeOperator{T}
-  coeff::DiffEqScalar{T,F}
-  op::OpType
-end
-*(α::DiffEqScalar{T,F}, L::AbstractDiffEqLinearOperator{T}) where {T,F} = DiffEqScaledOperator(α, L)
-*(α::Number, L::AbstractDiffEqLinearOperator{T}) where T = DiffEqScaledOperator(DiffEqScalar(convert(T,α)), L)
--(L::AbstractDiffEqLinearOperator{T}) where {T} = DiffEqScalar(-one(T)) * L
-getops(L::DiffEqScaledOperator) = (L.coeff, L.op)
-Matrix(L::DiffEqScaledOperator) = L.coeff * Matrix(L.op)
-convert(::Type{AbstractMatrix}, L::DiffEqScaledOperator) = L.coeff * convert(AbstractMatrix, L.op)
-
-size(L::DiffEqScaledOperator, args...) = size(L.op, args...)
-opnorm(L::DiffEqScaledOperator, p::Real=2) = abs(L.coeff) * opnorm(L.op, p)
-getindex(L::DiffEqScaledOperator, i::Int) = L.coeff * L.op[i]
-getindex(L::DiffEqScaledOperator, I::Vararg{Int, N}) where {N} =
-  L.coeff * L.op[I...]
-*(L::DiffEqScaledOperator, x::AbstractArray) = L.coeff * (L.op * x)
-*(x::AbstractArray, L::DiffEqScaledOperator) = (L.op * x) * L.coeff
-/(L::DiffEqScaledOperator, x::AbstractArray) = L.coeff * (L.op / x)
-/(x::AbstractArray, L::DiffEqScaledOperator) = 1/L.coeff * (x / L.op)
-\(L::DiffEqScaledOperator, x::AbstractArray) = 1/L.coeff * (L.op \ x)
-\(x::AbstractArray, L::DiffEqScaledOperator) = L.coeff * (x \ L)
-for N in (2,3)
-  @eval begin
-    mul!(Y::AbstractArray{T,$N}, L::DiffEqScaledOperator{T}, B::AbstractArray{T,$N}) where {T} =
-        lmul!(Y, L.coeff, mul!(Y, L.op, B))
-  end
-end
-ldiv!(Y::AbstractArray, L::DiffEqScaledOperator, B::AbstractArray) =
-  lmul!(1/L.coeff, ldiv!(Y, L.op, B))
-factorize(L::DiffEqScaledOperator) = L.coeff * factorize(L.op)
-for fact in (:lu, :lu!, :qr, :qr!, :cholesky, :cholesky!, :ldlt, :ldlt!,
-  :bunchkaufman, :bunchkaufman!, :lq, :lq!, :svd, :svd!)
-  @eval LinearAlgebra.$fact(L::DiffEqScaledOperator, args...) =
-    L.coeff * fact(L.op, args...)
-end
 
 # Linear Combination
 struct DiffEqOperatorCombination{T,O<:Tuple{Vararg{AbstractDiffEqLinearOperator{T}}},

test/DerivativeOperators/2D_3D_fast_multiplication.jl

@@ -1352,7 +1352,7 @@ end
     mul!(M_temp, A, M)
 
     # Numerical errors accumulating for this test case, more so than the other tests
-    @test_broken M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])
+    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:]) rtol=1e-5
 
     # Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M, no coefficient
     A = Lx2 + Lz2 + Lx3 + Lz3 + Lx4 + Lz4