From 0ebe1c613d5f367ec5cf9c23ce965edea51079f1 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 23 Sep 2022 11:40:53 -0400 Subject: [PATCH 1/7] some missing methods --- src/SciMLOperators.jl | 2 +- src/basic.jl | 11 ++++++++++- src/interface.jl | 9 --------- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/SciMLOperators.jl b/src/SciMLOperators.jl index d6cbc437..e1041d58 100644 --- a/src/SciMLOperators.jl +++ b/src/SciMLOperators.jl @@ -40,8 +40,8 @@ include("left.jl") include("multidim.jl") include("scalar.jl") -include("basic.jl") include("matrix.jl") +include("basic.jl") include("batch.jl") include("func.jl") include("tensor.jl") diff --git a/src/basic.jl b/src/basic.jl index 64d1a8c8..69cdd55d 100644 --- a/src/basic.jl +++ b/src/basic.jl @@ -311,7 +311,14 @@ end AddedOperator(L::AbstractSciMLOperator) = L # constructors -Base.:+(ops::AbstractSciMLOperator...) = AddedOperator(ops...) +function Base.:+(ops::Union{AbstractSciMLOperator, AbstractMatrix}...) + ops_ = () + for op in ops + op = op isa AbstractMatrix ? MatrixOperator(op) : op + ops_ = (ops_..., op) + end + AddedOperator(ops_...) +end Base.:+(A::AbstractSciMLOperator, B::AddedOperator) = AddedOperator(A, B.ops...) Base.:+(A::AddedOperator, B::AbstractSciMLOperator) = AddedOperator(A.ops..., B) Base.:+(A::AddedOperator, B::AddedOperator) = AddedOperator(A.ops..., B.ops...) @@ -327,6 +334,8 @@ function Base.:+(Z::NullOperator, A::AddedOperator) end Base.:-(A::AbstractSciMLOperator, B::AbstractSciMLOperator) = AddedOperator(A, -B) +Base.:-(A::AbstractSciMLOperator, B::AbstractMatrix) = A - MatrixOperator(B) +Base.:-(A::AbstractMatrix, B::AbstractSciMLOperator) = MatrixOperator(A) - B for op in ( :+, :-, diff --git a/src/interface.jl b/src/interface.jl index eac4f888..123bdc29 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -214,15 +214,6 @@ for op in ( end end -for op in ( - :+, :-, - ) - - @eval function Base.$op(L::AbstractSciMLLinearOperator, u::AbstractVecOrMat) - $op(convert(AbstractMatrix,L), u) - end -end - function LinearAlgebra.mul!(v::AbstractVecOrMat, L::AbstractSciMLLinearOperator, u::AbstractVecOrMat) mul!(v, convert(AbstractMatrix,L), u) end From b22dbedf5bfa339ed7072c52b68610e296a65db4 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 23 Sep 2022 12:19:20 -0400 Subject: [PATCH 2/7] more fixes; https://github.com/SciML/SciMLOperators.jl/issues/106 --- src/basic.jl | 2 +- src/interface.jl | 24 ++++++++++++------------ src/matrix.jl | 4 ++-- src/multidim.jl | 8 ++++---- 4 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/basic.jl b/src/basic.jl index 69cdd55d..960ce5f2 100644 --- a/src/basic.jl +++ b/src/basic.jl @@ -695,7 +695,7 @@ function LinearAlgebra.mul!(v::AbstractVecOrMat, L::InvertedOperator, u::Abstrac axpy!(β, L.cache, v) end -function LinearAlgebra.ldiv!(v::AbstractVecOrMat, L::InvertedOperator, u) +function LinearAlgebra.ldiv!(v::AbstractVecOrMat, L::InvertedOperator, u::AbstractVecOrMat) mul!(v, L.L, u) end diff --git a/src/interface.jl b/src/interface.jl index 123bdc29..39e596d7 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -167,12 +167,12 @@ issquare(::Union{ issquare(A...) = @. (&)(issquare(A)...) Base.isreal(L::AbstractSciMLOperator{T}) where{T} = T <: Real -Base.Matrix(L::AbstractSciMLLinearOperator) = Matrix(convert(AbstractMatrix, L)) +Base.Matrix(L::AbstractSciMLOperator) = Matrix(convert(AbstractMatrix, L)) -LinearAlgebra.exp(L::AbstractSciMLLinearOperator,t) = exp(t*L) +LinearAlgebra.exp(L::AbstractSciMLOperator,t) = exp(t*L) has_exp(L::AbstractSciMLLinearOperator) = true -expmv(L::AbstractSciMLLinearOperator,u,p,t) = exp(L,t)*u -expmv!(v,L::AbstractSciMLLinearOperator,u,p,t) = mul!(v,exp(L,t),u) +expmv(L::AbstractSciMLOperator,u,p,t) = exp(L,t)*u +expmv!(v,L::AbstractSciMLOperator,u,p,t) = mul!(v,exp(L,t),u) ### # fallback implementations @@ -188,37 +188,37 @@ function Base.:(==)(L1::AbstractSciMLOperator, L2::AbstractSciMLOperator) convert(AbstractMatrix, L1) == convert(AbstractMatrix, L1) end -Base.@propagate_inbounds function Base.getindex(L::AbstractSciMLLinearOperator, I::Vararg{Any,N}) where {N} +Base.@propagate_inbounds function Base.getindex(L::AbstractSciMLOperator, I::Vararg{Any,N}) where {N} convert(AbstractMatrix, L)[I...] end -function Base.getindex(L::AbstractSciMLLinearOperator, I::Vararg{Int, N}) where {N} +function Base.getindex(L::AbstractSciMLOperator, I::Vararg{Int, N}) where {N} convert(AbstractMatrix,L)[I...] end -LinearAlgebra.exp(L::AbstractSciMLLinearOperator) = exp(Matrix(L)) -LinearAlgebra.opnorm(L::AbstractSciMLLinearOperator, p::Real=2) = opnorm(convert(AbstractMatrix,L), p) +LinearAlgebra.exp(L::AbstractSciMLOperator) = exp(Matrix(L)) +LinearAlgebra.opnorm(L::AbstractSciMLOperator, p::Real=2) = opnorm(convert(AbstractMatrix,L), p) for pred in ( :issymmetric, :ishermitian, :isposdef, ) - @eval function LinearAlgebra.$pred(L::AbstractSciMLLinearOperator) + @eval function LinearAlgebra.$pred(L::AbstractSciMLOperator) $pred(convert(AbstractMatrix, L)) end end for op in ( :sum,:prod ) - @eval function LinearAlgebra.$op(L::AbstractSciMLLinearOperator; kwargs...) + @eval function LinearAlgebra.$op(L::AbstractSciMLOperator; kwargs...) $op(convert(AbstractMatrix, L); kwargs...) end end -function LinearAlgebra.mul!(v::AbstractVecOrMat, L::AbstractSciMLLinearOperator, u::AbstractVecOrMat) +function LinearAlgebra.mul!(v::AbstractVecOrMat, L::AbstractSciMLOperator, u::AbstractVecOrMat) mul!(v, convert(AbstractMatrix,L), u) end -function LinearAlgebra.mul!(v::AbstractVecOrMat, L::AbstractSciMLLinearOperator, u::AbstractVecOrMat, α, β) +function LinearAlgebra.mul!(v::AbstractVecOrMat, L::AbstractSciMLOperator, u::AbstractVecOrMat, α, β) mul!(v, convert(AbstractMatrix,L), u, α, β) end # diff --git a/src/matrix.jl b/src/matrix.jl index a27827f9..98839d2a 100644 --- a/src/matrix.jl +++ b/src/matrix.jl @@ -142,9 +142,9 @@ for fact in ( :svd, :svd!, ) - @eval LinearAlgebra.$fact(L::AbstractSciMLLinearOperator, args...) = + @eval LinearAlgebra.$fact(L::AbstractSciMLOperator, args...) = InvertibleOperator($fact(convert(AbstractMatrix, L), args...)) - @eval LinearAlgebra.$fact(L::AbstractSciMLLinearOperator; kwargs...) = + @eval LinearAlgebra.$fact(L::AbstractSciMLOperator; kwargs...) = InvertibleOperator($fact(convert(AbstractMatrix, L); kwargs...)) end diff --git a/src/multidim.jl b/src/multidim.jl index cc4d6124..00e7d084 100644 --- a/src/multidim.jl +++ b/src/multidim.jl @@ -21,7 +21,7 @@ for op in ( end end -function LinearAlgebra.mul!(v::AbstractArray, L::AbstractSciMLLinearOperator, u::AbstractArray) +function LinearAlgebra.mul!(v::AbstractArray, L::AbstractSciMLOperator, u::AbstractArray) u isa AbstractVecOrMat && @error "LinearAlgebra.mul! not defined for $(typeof(L)), $(typeof(u))." sizes = _mat_sizes(L, u) @@ -34,7 +34,7 @@ function LinearAlgebra.mul!(v::AbstractArray, L::AbstractSciMLLinearOperator, u: v end -function LinearAlgebra.mul!(v::AbstractArray, L::AbstractSciMLLinearOperator, u::AbstractArray, α, β) +function LinearAlgebra.mul!(v::AbstractArray, L::AbstractSciMLOperator, u::AbstractArray, α, β) u isa AbstractVecOrMat && @error "LinearAlgebra.mul! not defined for $(typeof(L)), $(typeof(u))." sizes = _mat_sizes(L, u) @@ -47,7 +47,7 @@ function LinearAlgebra.mul!(v::AbstractArray, L::AbstractSciMLLinearOperator, u: v end -function LinearAlgebra.ldiv!(v::AbstractArray, L::AbstractSciMLLinearOperator, u::AbstractArray) +function LinearAlgebra.ldiv!(v::AbstractArray, L::AbstractSciMLOperator, u::AbstractArray) u isa AbstractVecOrMat && @error "LinearAlgebra.ldiv! not defined for $(typeof(L)), $(typeof(u))." sizes = _mat_sizes(L, u) @@ -60,7 +60,7 @@ function LinearAlgebra.ldiv!(v::AbstractArray, L::AbstractSciMLLinearOperator, u v end -function LinearAlgebra.ldiv!(L::AbstractSciMLLinearOperator, u::AbstractArray) +function LinearAlgebra.ldiv!(L::AbstractSciMLOperator, u::AbstractArray) u isa AbstractVecOrMat && @error "LinearAlgebra.ldiv! not defined for $(typeof(L)), $(typeof(u))." sizes = _mat_sizes(L, u) From e29e1a16cfad1d4c640689c9e7d4357e1eb123e8 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 23 Sep 2022 12:33:27 -0400 Subject: [PATCH 3/7] bug fix --- src/basic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/basic.jl b/src/basic.jl index 960ce5f2..d278fcf1 100644 --- a/src/basic.jl +++ b/src/basic.jl @@ -212,7 +212,7 @@ end Base.:-(L::AbstractSciMLOperator) = ScaledOperator(-true, L) Base.:+(L::AbstractSciMLOperator) = L -Base.convert(::Type{AbstractMatrix}, L::ScaledOperator) = L.λ.val * convert(AbstractMatrix, L.L) +Base.convert(::Type{AbstractMatrix}, L::ScaledOperator) = convert(Number,L.λ) * convert(AbstractMatrix, L.L) SparseArrays.sparse(L::ScaledOperator) = L.λ * sparse(L.L) # traits From 1e69def250509be8a6001fc23c27a322d99d129d Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 23 Sep 2022 12:36:13 -0400 Subject: [PATCH 4/7] bump --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8b7961d1..a51757d0 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SciMLOperators" uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" authors = ["xtalax "] -version = "0.1.13" +version = "0.1.14" [deps] ArrayInterfaceCore = "30b0a656-2188-435a-8636-2ec0e6a096e2" From 54a98b07eda7350ba7bce6a954446cc851caeb52 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 23 Sep 2022 12:42:59 -0400 Subject: [PATCH 5/7] convert method for factorization --- src/matrix.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/matrix.jl b/src/matrix.jl index 98839d2a..f915c750 100644 --- a/src/matrix.jl +++ b/src/matrix.jl @@ -148,6 +148,10 @@ for fact in ( InvertibleOperator($fact(convert(AbstractMatrix, L); kwargs...)) end +function Base.convert(::Type{<:Factorization}, L::InvertibleOperator{T,<:Factorization}) where{T} + L.F +end + function Base.convert(::Type{AbstractMatrix}, L::InvertibleOperator) if L.F isa Adjoint convert(AbstractMatrix,L.F')' From 1453743741da2c08042fa8c4df10e75ae28b312f Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 23 Sep 2022 13:37:01 -0400 Subject: [PATCH 6/7] rm _sparse --- src/basic.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/basic.jl b/src/basic.jl index d278fcf1..e0301d48 100644 --- a/src/basic.jl +++ b/src/basic.jl @@ -359,7 +359,7 @@ for op in ( end Base.convert(::Type{AbstractMatrix}, L::AddedOperator) = sum(op -> convert(AbstractMatrix, op), L.ops) -SparseArrays.sparse(L::AddedOperator) = sum(_sparse, L.ops) +SparseArrays.sparse(L::AddedOperator) = sum(sparse, L.ops) # traits Base.size(L::AddedOperator) = size(first(L.ops)) @@ -491,7 +491,7 @@ for op in ( end Base.convert(::Type{AbstractMatrix}, L::ComposedOperator) = prod(op -> convert(AbstractMatrix, op), L.ops) -SparseArrays.sparse(L::ComposedOperator) = prod(_sparse, L.ops) +SparseArrays.sparse(L::ComposedOperator) = prod(sparse, L.ops) # traits Base.size(L::ComposedOperator) = (size(first(L.ops), 1), size(last(L.ops),2)) From 3f36ae2bd07dc3f00ba41703b881cc16ae330f1c Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 23 Sep 2022 15:07:17 -0400 Subject: [PATCH 7/7] https://github.com/SciML/SciMLOperators.jl/issues/107 --- src/SciMLOperators.jl | 2 +- src/matrix.jl | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/SciMLOperators.jl b/src/SciMLOperators.jl index e1041d58..27b15da4 100644 --- a/src/SciMLOperators.jl +++ b/src/SciMLOperators.jl @@ -17,7 +17,7 @@ import Base: +, -, *, /, \, ∘, ==, conj, exp, kron import Base: iszero, inv, adjoint, transpose, size, convert import LinearAlgebra: mul!, ldiv!, lmul!, rmul!, factorize import LinearAlgebra: Matrix, Diagonal -import SparseArrays: sparse +import SparseArrays: sparse, issparse """ $(TYPEDEF) diff --git a/src/matrix.jl b/src/matrix.jl index f915c750..15c848a1 100644 --- a/src/matrix.jl +++ b/src/matrix.jl @@ -52,6 +52,7 @@ isconstant(L::MatrixOperator) = L.update_func == DEFAULT_UPDATE_FUNC Base.iszero(L::MatrixOperator) = iszero(L.A) SparseArrays.sparse(L::MatrixOperator) = sparse(L.A) +SparseArrays.issparse(L::MatrixOperator) = issparse(L.A) # TODO - add tests for MatrixOperator indexing # propagate_inbounds here for the getindex fallback