Merge pull request #44 from SebastianM-C/move_scale

Move `DiffEqScaledOperator`

Author: ChrisRackauckas

src/SciMLBase.jl

@@ -442,6 +442,11 @@ $(TYPEDEF)
 """
 abstract type AbstractDiffEqLinearOperator{T} <: AbstractDiffEqOperator{T} end
 
+"""
+$(TYPEDEF)
+"""
+abstract type AbstractDiffEqCompositeOperator{T} <: AbstractDiffEqLinearOperator{T} end
+
 """
 $(TYPEDEF)
 
@@ -488,9 +493,9 @@ include("scimlfunctions.jl")
 include("alg_traits.jl")
 
 include("operators/operators.jl")
+include("operators/basic_operators.jl")
 include("operators/diffeq_operator.jl")
 include("operators/common_defaults.jl")
-include("operators/basic_operators.jl")
 
 include("problems/problem_utils.jl")
 include("problems/discrete_problems.jl")
@@ -573,7 +578,7 @@ export EnsembleAnalysis, EnsembleSummary
 
 export tuples, intervals, TimeChoiceIterator
 
-export AffineDiffEqOperator
+export AffineDiffEqOperator, DiffEqScaledOperator
 
 export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity
 

src/operators/common_defaults.jl

@@ -10,16 +10,16 @@ LinearAlgebra.opnorm(L::AbstractDiffEqLinearOperator, p::Real=2) = opnorm(conver
 Base.@propagate_inbounds Base.getindex(L::AbstractDiffEqLinearOperator, I::Vararg{Any,N}) where {N} = convert(AbstractMatrix,L)[I...]
 Base.getindex(L::AbstractDiffEqLinearOperator, I::Vararg{Int, N}) where {N} =
   convert(AbstractMatrix,L)[I...]
-for op in (:*, :/, :\), T in (:AbstractArray, :Number)
-  @eval Base.$op(L::AbstractDiffEqLinearOperator, x::$T) = $op(convert(AbstractMatrix,L), x)
-  @eval Base.$op(x::$T, L::AbstractDiffEqLinearOperator) = $op(x, convert(AbstractMatrix,L))
+for op in (:*, :/, :\)
+  @eval Base.$op(L::AbstractDiffEqLinearOperator, x::AbstractArray) = $op(convert(AbstractMatrix,L), x)
+  @eval Base.$op(x::AbstractArray, L::AbstractDiffEqLinearOperator) = $op(x, convert(AbstractMatrix,L))
+  @eval Base.$op(L::DiffEqArrayOperator, x::Number) = $op(convert(AbstractMatrix,L), x)
+  @eval Base.$op(x::Number, L::DiffEqArrayOperator) = $op(x, convert(AbstractMatrix,L))
 end
 LinearAlgebra.mul!(Y::AbstractArray, L::AbstractDiffEqLinearOperator, B::AbstractArray) =
   mul!(Y, convert(AbstractMatrix,L), B)
 LinearAlgebra.mul!(Y::AbstractArray, L::AbstractDiffEqLinearOperator, B::AbstractArray, α::Number, β::Number) =
   mul!(Y, convert(AbstractMatrix,L), B, α, β)
-LinearAlgebra.ldiv!(Y::AbstractVecOrMat, L::AbstractDiffEqLinearOperator, B::AbstractVecOrMat) =
-  ldiv!(Y, convert(AbstractMatrix,L), B)
 for pred in (:isreal, :issymmetric, :ishermitian, :isposdef)
   @eval LinearAlgebra.$pred(L::AbstractDiffEqLinearOperator) = $pred(convert(AbstractArray, L))
 end

src/operators/diffeq_operator.jl

@@ -74,3 +74,52 @@ function update_coefficients!(L::AffineDiffEqOperator,u,p,t)
 end
 
 @deprecate is_constant(L::AbstractDiffEqOperator) isconstant(L)
+
+# Scaled operator (α * A)
+struct DiffEqScaledOperator{T,F,OpType<:AbstractDiffEqLinearOperator{T}} <: AbstractDiffEqCompositeOperator{T}
+  coeff::DiffEqScalar{T,F}
+  op::OpType
+end
+
+# 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
+
+Base.:*(α::DiffEqScalar{T,F}, L::AbstractDiffEqLinearOperator{T}) where {T,F} = DiffEqScaledOperator(α, L)
+Base.:*(α::Number, L::AbstractDiffEqLinearOperator{T}) where T = DiffEqScaledOperator(DiffEqScalar(convert(T,α)), L)
+Base.:-(L::AbstractDiffEqLinearOperator{T}) where {T} = DiffEqScalar(-one(T)) * L
+getops(L::DiffEqScaledOperator) = (L.coeff, L.op)
+Base.Matrix(L::DiffEqScaledOperator) = L.coeff * Matrix(L.op)
+Base.convert(::Type{AbstractMatrix}, L::DiffEqScaledOperator) = L.coeff * convert(AbstractMatrix, L.op)
+
+Base.size(L::DiffEqScaledOperator, args...) = size(L.op, args...)
+LinearAlgebra.opnorm(L::DiffEqScaledOperator, p::Real=2) = abs(L.coeff) * opnorm(L.op, p)
+Base.getindex(L::DiffEqScaledOperator, i::Int) = L.coeff * L.op[i]
+Base.getindex(L::DiffEqScaledOperator, I::Vararg{Int, N}) where {N} =
+  L.coeff * L.op[I...]
+Base.:*(L::DiffEqScaledOperator, x::AbstractArray) = L.coeff * (L.op * x)
+Base.:*(x::AbstractArray, L::DiffEqScaledOperator) = (L.op * x) * L.coeff
+Base.:/(L::DiffEqScaledOperator, x::AbstractArray) = L.coeff * (L.op / x)
+Base.:/(x::AbstractArray, L::DiffEqScaledOperator) = 1/L.coeff * (x / L.op)
+Base.:\(L::DiffEqScaledOperator, x::AbstractArray) = 1/L.coeff * (L.op \ x)
+Base.:\(x::AbstractArray, L::DiffEqScaledOperator) = L.coeff * (x \ L)
+for N in (2,3)
+  @eval begin
+    LinearAlgebra.mul!(Y::AbstractArray{T,$N}, L::DiffEqScaledOperator{T}, B::AbstractArray{T,$N}) where {T} =
+        LinearAlgebra.lmul!(Y, L.coeff, mul!(Y, L.op, B))
+  end
+end
+LinearAlgebra.ldiv!(Y::AbstractArray, L::DiffEqScaledOperator, B::AbstractArray) =
+  lmul!(1/L.coeff, ldiv!(Y, L.op, B))
+LinearAlgebra.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
+
+isconstant(L::AbstractDiffEqCompositeOperator) = all(isconstant, getops(L))