Merge pull request #542 from vpuri3/smo

Fix method ambiguities arising due to SciMLOperators integration

Author: ChrisRackauckas

.gitignore

@@ -5,4 +5,6 @@ deps/deps.jl
 Manifest.toml
 *.png
 *.gif
-docs/build
+docs/build
+.*.swo
+.*.swp

Project.toml

@@ -34,7 +34,7 @@ LoopVectorization = "0.12"
 NNlib = "0.6, 0.7, 0.8"
 Requires = "1"
 RuntimeGeneratedFunctions = "0.4, 0.5"
-SciMLBase = "1.11"
+SciMLBase = "1.83"
 SparseDiffTools = "1.17"
 StaticArrays = "0.10, 0.11, 0.12, 1.0, 1.10, 1.3.2"
 julia = "1.6"

src/composite_operators.jl

@@ -61,10 +61,18 @@ getindex(L::DiffEqOperatorCombination, i::Int) = sum(op -> op[i], L.ops)
 function getindex(L::DiffEqOperatorCombination, I::Vararg{Int, N}) where {N}
     sum(op -> op[I...], L.ops)
 end
+Base.:*(L::DiffEqOperatorCombination, x::AbstractVecOrMat) = sum(op -> op * x, L.ops)
 Base.:*(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op * x, L.ops)
+
+Base.:*(x::AbstractVecOrMat, L::DiffEqOperatorCombination) = sum(op -> x * op, L.ops)
 Base.:*(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x * op, L.ops)
+
+/(L::DiffEqOperatorCombination, x::AbstractVecOrMat) = sum(op -> op / x, L.ops)
 /(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op / x, L.ops)
+
+\(x::AbstractVecOrMat, L::DiffEqOperatorCombination) = sum(op -> x \ op, L.ops)
 \(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x \ op, L.ops)
+
 function mul!(y::AbstractVector, L::DiffEqOperatorCombination, b::AbstractVector)
     mul!(y, L.ops[1], b)
     for op in L.ops[2:end]
@@ -138,24 +146,49 @@ SparseArrays.sparse(L::DiffEqOperatorComposition) = prod(sparse1, reverse(L.ops)
 size(L::DiffEqOperatorComposition) = (size(L.ops[end], 1), size(L.ops[1], 2))
 size(L::DiffEqOperatorComposition, m::Integer) = size(L)[m]
 opnorm(L::DiffEqOperatorComposition) = prod(opnorm, L.ops)
+
+function Base.:*(L::DiffEqOperatorComposition, x::AbstractVecOrMat)
+    foldl((acc, op) -> op * acc, L.ops; init = x)
+end
 function Base.:*(L::DiffEqOperatorComposition, x::AbstractArray)
     foldl((acc, op) -> op * acc, L.ops; init = x)
 end
+
+function Base.:*(x::AbstractVecOrMat, L::DiffEqOperatorComposition)
+    foldl((acc, op) -> acc * op, reverse(L.ops); init = x)
+end
 function Base.:*(x::AbstractArray, L::DiffEqOperatorComposition)
     foldl((acc, op) -> acc * op, reverse(L.ops); init = x)
 end
+
+function /(L::DiffEqOperatorComposition, x::AbstractVecOrMat)
+    foldl((acc, op) -> op / acc, L.ops; init = x)
+end
 function /(L::DiffEqOperatorComposition, x::AbstractArray)
     foldl((acc, op) -> op / acc, L.ops; init = x)
 end
+
+function /(x::AbstractVecOrMat, L::DiffEqOperatorComposition)
+    foldl((acc, op) -> acc / op, L.ops; init = x)
+end
 function /(x::AbstractArray, L::DiffEqOperatorComposition)
     foldl((acc, op) -> acc / op, L.ops; init = x)
 end
+
+function \(L::DiffEqOperatorComposition, x::AbstractVecOrMat)
+    foldl((acc, op) -> op \ acc, reverse(L.ops); init = x)
+end
 function \(L::DiffEqOperatorComposition, x::AbstractArray)
     foldl((acc, op) -> op \ acc, reverse(L.ops); init = x)
 end
+
+function \(x::AbstractVecOrMat, L::DiffEqOperatorComposition)
+    foldl((acc, op) -> acc \ op, reverse(L.ops); init = x)
+end
 function \(x::AbstractArray, L::DiffEqOperatorComposition)
     foldl((acc, op) -> acc \ op, reverse(L.ops); init = x)
 end
+
 function mul!(y::AbstractVector, L::DiffEqOperatorComposition, b::AbstractVector)
     mul!(L.caches[1], L.ops[1], b)
     for i in 2:(length(L.ops) - 1)

src/derivative_operators/derivative_operator_functions.jl

@@ -15,8 +15,17 @@
 #
 
 # Fallback mul! implementation for a single DerivativeOperator operating on an AbstractArray
+function LinearAlgebra.mul!(x_temp::AbstractVecOrMat{T}, A::DerivativeOperator{T, N},
+                            M::AbstractVecOrMat{T}; overwrite = true) where {T, N}
+    _mul!(x_temp, A, M; overwrite = overwrite)
+end
 function LinearAlgebra.mul!(x_temp::AbstractArray{T}, A::DerivativeOperator{T, N},
                             M::AbstractArray{T}; overwrite = true) where {T, N}
+    _mul!(x_temp, A, M; overwrite = overwrite)
+end
+
+function _mul!(x_temp::AbstractArray{T}, A::DerivativeOperator{T, N},
+                            M::AbstractArray{T}; overwrite = true) where {T, N}
 
     # Check that x_temp has valid dimensions, allowing unnecessary padding in M
     v = zeros(ndims(x_temp))
@@ -138,7 +147,14 @@ end
 
 ###########################################
 
-function *(A::DerivativeOperator{T, N}, M::AbstractArray{T}) where {T <: Real, N}
+function Base.:*(A::DerivativeOperator{T, N}, M::AbstractVecOrMat{T}) where {T <: Real, N}
+    _mul(A, M)
+end
+function Base.:*(A::DerivativeOperator{T, N}, M::AbstractArray{T}) where {T <: Real, N}
+    _mul(A, M)
+end
+
+function _mul(A::DerivativeOperator{T, N}, M::AbstractArray{T}) where {T <: Real, N}
     size_x_temp = [size(M)...]
     size_x_temp[N] -= 2
     x_temp = zeros(promote_type(eltype(A), eltype(M)), size_x_temp...)

src/derivative_operators/ghost_derivative_operator.jl

@@ -12,32 +12,49 @@ function *(L::AbstractDiffEqCompositeOperator{T}, Q::AbstractBC{T}) where {T}
     return sum(map(op -> op * Q, L.ops))
 end
 
+function LinearAlgebra.mul!(out::AbstractVecOrMat, A::GhostDerivativeOperator,
+                            u::AbstractVecOrMat)
+    padded = A.Q * u  # assume: creates boundary padded array w/o realloc
+    LinearAlgebra.mul!(out, A.L, padded)
+end
 function LinearAlgebra.mul!(out::AbstractArray, A::GhostDerivativeOperator,
                             u::AbstractArray)
     padded = A.Q * u  # assume: creates boundary padded array w/o realloc
     LinearAlgebra.mul!(out, A.L, padded)
 end
 
+function *(A::GhostDerivativeOperator{T1}, u::AbstractVecOrMat{T2}) where {T1, T2}
+    #TODO Implement a function domaincheck(L::AbstractDiffEqLinearOperator, u) to see if components of L along each dimension match the size of u
+    x = zeros(promote_type(T1, T2), unpadded_size(u))
+    LinearAlgebra.mul!(x, A, u)
+    return x
+end
 function *(A::GhostDerivativeOperator{T1}, u::AbstractArray{T2}) where {T1, T2}
     #TODO Implement a function domaincheck(L::AbstractDiffEqLinearOperator, u) to see if components of L along each dimension match the size of u
     x = zeros(promote_type(T1, T2), unpadded_size(u))
     LinearAlgebra.mul!(x, A, u)
     return x
 end
 
-function \(A::GhostDerivativeOperator, u::AbstractArray) # FIXME should have T1,T2 and promote result
+function \(A::GhostDerivativeOperator, u::AbstractVector) # FIXME as above, should promote_type(T1,T2)
+    # TODO: is this specialization to u::AbstractVector really any faster?
+    @assert length(u) == size(A.L, 1)
+    (A_l, A_b) = sparse(A, length(u))
+    A_l \ Vector(u .- A_b)
+end
+function \(A::GhostDerivativeOperator, u::AbstractMatrix) # FIXME should have T1,T2 and promote result
     #TODO implement check that A has compatible size with u
     s = size(u)
     (A_l, A_b) = sparse(A, s)
     x = A_l \ Vector(reshape(u, length(u)) .- A_b) #Has to be converted to vector to work, A_b being sparse was causing a conversion to sparse.
     return reshape(x, s)
 end
-
-function \(A::GhostDerivativeOperator, u::AbstractVector) # FIXME as above, should promote_type(T1,T2)
-    # TODO: is this specialization to u::AbstractVector really any faster?
-    @assert length(u) == size(A.L, 1)
-    (A_l, A_b) = sparse(A, length(u))
-    A_l \ Vector(u .- A_b)
+function \(A::GhostDerivativeOperator, u::AbstractArray) # FIXME should have T1,T2 and promote result
+    #TODO implement check that A has compatible size with u
+    s = size(u)
+    (A_l, A_b) = sparse(A, s)
+    x = A_l \ Vector(reshape(u, length(u)) .- A_b) #Has to be converted to vector to work, A_b being sparse was causing a conversion to sparse.
+    return reshape(x, s)
 end
 
 # update coefficients

src/jacvec_operators.jl

@@ -1,6 +1,6 @@
 ### Operator Implementation
 mutable struct JacVecOperator{T, F, T1, T2, uType, P, tType, O} <:
-               DiffEqBase.AbstractDiffEqLinearOperator{T}
+               AbstractDiffEqLinearOperator{T}
     f::F
     cache1::T1
     cache2::T2
@@ -151,7 +151,7 @@ end
 ### AnalyticalOperator Implementation
 
 mutable struct AnalyticalJacVecOperator{T, F, uType, P, tType, O} <:
-               DiffEqBase.AbstractDiffEqLinearOperator{T}
+               AbstractDiffEqLinearOperator{T}
     f::F
     u::uType
     p::P

src/matrixfree_operators.jl

@@ -20,7 +20,7 @@ function Base.size(M::MatrixFreeOperator)
         error("M.size is nothing, please define size as a tuple of integers")
     return M.size
 end
-@inline function Base.size(M::MatrixFreeOperator, n)
+@inline function Base.size(M::MatrixFreeOperator, n::Integer)
     M.size === nothing &&
         error("M.size is nothing, please define size as a tuple of integers")
     n <= 0 && error("dimension out of range")

src/vecjac_operators.jl

@@ -1,5 +1,5 @@
 mutable struct VecJacOperator{T, F, T1, T2, uType, P, tType, O} <:
-               DiffEqBase.AbstractDiffEqLinearOperator{T}
+               AbstractDiffEqLinearOperator{T}
     f::F
     cache1::T1
     cache2::T2