Merge pull request #139 from vpuri3/islinear

remove abstractscimllinearoperator in favor or trait islinear(L)

Author: xtalax

src/SciMLOperators.jl

@@ -26,12 +26,7 @@ abstract type AbstractSciMLOperator{T} end
 """
 $(TYPEDEF)
 """
-abstract type AbstractSciMLLinearOperator{T} <: AbstractSciMLOperator{T} end
-
-"""
-$(TYPEDEF)
-"""
-abstract type AbstractSciMLScalarOperator{T} <: AbstractSciMLLinearOperator{T} end
+abstract type AbstractSciMLScalarOperator{T} <: AbstractSciMLOperator{T} end
 
 include("utils.jl")
 include("interface.jl")
@@ -58,6 +53,9 @@ export update_coefficients!,
 
        cache_operator,
 
+       issquare,
+       islinear,
+
        has_adjoint,
        has_expmv,
        has_expmv!,

src/basic.jl

@@ -1,7 +1,7 @@
 """
 $(TYPEDEF)
 """
-struct IdentityOperator{N} <: AbstractSciMLLinearOperator{Bool} end
+struct IdentityOperator{N} <: AbstractSciMLOperator{Bool} end
 
 # constructors
 IdentityOperator(u::AbstractArray) = IdentityOperator{size(u,1)}()
@@ -28,7 +28,7 @@ end
 
 getops(::IdentityOperator) = ()
 isconstant(::IdentityOperator) = true
-islinear(L::IdentityOperator) = true
+islinear(::IdentityOperator) = true
 has_adjoint(::IdentityOperator) = true
 has_mul!(::IdentityOperator) = true
 has_ldiv(::IdentityOperator) = true
@@ -92,7 +92,7 @@ end
 """
 $(TYPEDEF)
 """
-struct NullOperator{N} <: AbstractSciMLLinearOperator{Bool} end
+struct NullOperator{N} <: AbstractSciMLOperator{Bool} end
 
 # constructors
 NullOperator(u::AbstractArray) = NullOperator{size(u,1)}()
@@ -120,7 +120,7 @@ LinearAlgebra.isposdef(::NullOperator) = false
 
 getops(::NullOperator) = ()
 isconstant(::NullOperator) = true
-islinear(L::NullOperator) = true
+islinear(::NullOperator) = true
 Base.iszero(::NullOperator) = true
 has_adjoint(::NullOperator) = true
 has_mul!(::NullOperator) = true
@@ -229,8 +229,8 @@ Base.conj(L::ScaledOperator) = conj(L.λ) * conj(L.L)
 LinearAlgebra.opnorm(L::ScaledOperator, p::Real=2) = abs(L.λ) * opnorm(L.L, p)
 
 getops(L::ScaledOperator) = (L.λ, L.L,)
-islinear(L::ScaledOperator) = all(islinear, L.ops)
 isconstant(L::ScaledOperator) = isconstant(L.L) & isconstant(L.λ)
+islinear(L::ScaledOperator) = islinear(L.L)
 Base.iszero(L::ScaledOperator) = iszero(L.L) | iszero(L.λ)
 has_adjoint(L::ScaledOperator) = has_adjoint(L.L)
 has_mul(L::ScaledOperator) = has_mul(L.L)
@@ -375,6 +375,7 @@ end
 Base.conj(L::AddedOperator) = AddedOperator(conj.(L.ops))
 
 getops(L::AddedOperator) = L.ops
+islinear(L::AddedOperator) = all(islinear, getops(L))
 Base.iszero(L::AddedOperator) = all(iszero, getops(L))
 has_adjoint(L::AddedOperator) = all(has_adjoint, L.ops)
 
@@ -667,6 +668,7 @@ Base.adjoint(L::InvertedOperator) = InvertedOperator(adjoint(L.L); cache = L.iss
 Base.conj(L::InvertedOperator) = InvertedOperator(conj(L.L); cache=L.cache)
 
 getops(L::InvertedOperator) = (L.L,)
+islinear(L::InvertedOperator) = islinear(L.L)
 
 has_mul(L::InvertedOperator) = has_ldiv(L.L)
 has_mul!(L::InvertedOperator) = has_ldiv!(L.L)

src/batch.jl

@@ -8,7 +8,7 @@ by `update_coefficients!` and is assumed to have the following signature:
 
     update_func(diag::AbstractVector,u,p,t) -> [modifies diag]
 """
-struct BatchedDiagonalOperator{T,D,F} <: AbstractSciMLLinearOperator{T}
+struct BatchedDiagonalOperator{T,D,F} <: AbstractSciMLOperator{T}
     diag::D
     update_func::F
 
@@ -58,7 +58,7 @@ end
 LinearAlgebra.isposdef(L::BatchedDiagonalOperator) = isposdef(Diagonal(vec(L.diag)))
 
 isconstant(L::BatchedDiagonalOperator) = L.update_func == DEFAULT_UPDATE_FUNC
-issquare(L::BatchedDiagonalOperator) = true
+islinear(::BatchedDiagonalOperator) = true
 has_adjoint(L::BatchedDiagonalOperator) = true
 has_ldiv(L::BatchedDiagonalOperator) = all(x -> !iszero(x), L.diag)
 has_ldiv!(L::BatchedDiagonalOperator) = has_ldiv(L)

src/func.jl

@@ -104,10 +104,12 @@ function FunctionOperator(op,
                           t::Union{Number,Nothing}=nothing,
 
                           # traits
-                          opnorm=nothing,
-                          issymmetric::Bool=false,
-                          ishermitian::Bool=false,
-                          isposdef::Bool=false,
+                          islinear::Bool = false,
+
+                          opnorm = nothing,
+                          issymmetric::Bool = false,
+                          ishermitian::Bool = false,
+                          isposdef::Bool = false,
                          )
 
     sz = (size(output, 1), size(input, 1))
@@ -148,6 +150,8 @@ function FunctionOperator(op,
     end
 
     traits = (;
+              islinear = islinear,
+
               opnorm = opnorm,
               issymmetric = issymmetric,
               ishermitian = ishermitian,
@@ -282,6 +286,7 @@ LinearAlgebra.ishermitian(L::FunctionOperator) = L.traits.ishermitian
 LinearAlgebra.isposdef(L::FunctionOperator) = L.traits.isposdef
 
 getops(::FunctionOperator) = ()
+islinear(L::FunctionOperator) = L.traits.islinear
 has_adjoint(L::FunctionOperator) = !(L.op_adjoint isa Nothing)
 has_mul(L::FunctionOperator{iip}) where{iip} = true
 has_mul!(L::FunctionOperator{iip}) where{iip} = iip

src/interface.jl

@@ -75,17 +75,12 @@ Base.eltype(::AbstractSciMLOperator{T}) where T = T
 Base.oneunit(L::AbstractSciMLOperator) = one(L)
 Base.oneunit(LType::Type{<:AbstractSciMLOperator}) = one(LType)
 
-issquare(L::AbstractSciMLOperator) = isequal(size(L)...)
-
-Base.iszero(::AbstractSciMLOperator) = false
-
-islinear(::AbstractSciMLOperator) = false
-islinear(::AbstractSciMLLinearOperator) = true
+Base.iszero(::AbstractSciMLOperator) = false # TODO
 
 has_adjoint(L::AbstractSciMLOperator) = false # L', adjoint(L)
 has_expmv!(L::AbstractSciMLOperator) = false # expmv!(v, L, t, u)
 has_expmv(L::AbstractSciMLOperator) = false # v = exp(L, t, u)
-has_exp(L::AbstractSciMLOperator) = false # v = exp(L, t)*u
+has_exp(L::AbstractSciMLOperator) = islinear(L)
 has_mul(L::AbstractSciMLOperator) = true # du = L*u
 has_mul!(L::AbstractSciMLOperator) = false # mul!(du, L, u)
 has_ldiv(L::AbstractSciMLOperator) = false # du = L\u
@@ -97,7 +92,9 @@ isconstant(L) = true
 isconstant(L::AbstractSciMLOperator) = all(isconstant, getops(L))
 #isconstant(L::AbstractSciMLOperator) = L.update_func = DEFAULT_UPDATE_FUNC
 
-islinear(L) = false
+#islinear(L) = false
+islinear(::AbstractSciMLOperator) = false
+
 islinear(::Union{
                  # LinearAlgebra
                  AbstractMatrix,
@@ -106,23 +103,20 @@ islinear(::Union{
 
                  # Base
                  Number,
-
-                 # SciMLOperators
-                 AbstractSciMLLinearOperator,
                 }
         ) = true
 
 has_mul(L) = false
 has_mul(::Union{
-                 # LinearAlgebra
-                 AbstractVecOrMat,
-                 AbstractMatrix,
-                 UniformScaling,
+                # LinearAlgebra
+                AbstractVecOrMat,
+                AbstractMatrix,
+                UniformScaling,
 
-                 # Base
-                 Number,
-                }
-        ) = true
+                # Base
+                Number,
+               }
+       ) = true
 
 has_mul!(L) = false
 has_mul!(::Union{
@@ -152,7 +146,7 @@ has_ldiv!(::Union{
                  }
          ) = true
 
-has_adjoint(L) = false
+has_adjoint(L) = islinear(L)
 has_adjoint(::Union{
                     # LinearAlgebra
                     AbstractMatrix,
@@ -161,13 +155,10 @@ has_adjoint(::Union{
 
                     # Base
                     Number,
-
-                    # SciMLOperators
-                    AbstractSciMLLinearOperator,
                    }
            ) = true
 
-issquare(A) = size(A,1) === size(A,2)
+issquare(L) = isequal(size(L)...)
 issquare(::Union{
                  # LinearAlgebra
                  UniformScaling,
@@ -182,7 +173,6 @@ Base.isreal(L::AbstractSciMLOperator{T}) where{T} = T <: Real
 Base.Matrix(L::AbstractSciMLOperator) = Matrix(convert(AbstractMatrix, L))
 
 LinearAlgebra.exp(L::AbstractSciMLOperator,t) = exp(t*L)
-has_exp(L::AbstractSciMLLinearOperator) = true
 expmv(L::AbstractSciMLOperator,u,p,t) = exp(L,t)*u
 expmv!(v,L::AbstractSciMLOperator,u,p,t) = mul!(v,exp(L,t),u)
 

src/left.jl

@@ -60,6 +60,10 @@ AbstractAdjointVecOrMat    = Adjoint{  T,<:AbstractVecOrMat} where{T}
 AbstractTransposedVecOrMat = Transpose{T,<:AbstractVecOrMat} where{T}
 
 has_adjoint(::AdjointOperator) = true
+#has_adjoint(::TransposedOperator) = ??
+
+islinear(L::AdjointOperator) = islinear(L.L)
+islinear(L::TransposedOperator) = islinear(L.L)
 
 Base.transpose(L::AdjointOperator) = conj(L.L)
 Base.adjoint(L::TransposedOperator) = conj(L.L)

src/matrix.jl

@@ -8,7 +8,7 @@ the following signature:
 
     update_func(A::AbstractMatrix,u,p,t) -> [modifies A]
 """
-struct MatrixOperator{T,AType<:AbstractMatrix{T},F} <: AbstractSciMLLinearOperator{T}
+struct MatrixOperator{T,AType<:AbstractMatrix{T},F} <: AbstractSciMLOperator{T}
     A::AType
     update_func::F
     MatrixOperator(A::AType; update_func=DEFAULT_UPDATE_FUNC) where{AType} =
@@ -28,6 +28,8 @@ Base.similar(L::MatrixOperator, ::Type{T}, dims::Dims) where{T} = MatrixOperator
                            has_ldiv,
                            has_ldiv!,
                           )
+islinear(::MatrixOperator) = true
+
 Base.size(L::MatrixOperator) = size(L.A)
 for op in (
            :adjoint,
@@ -121,7 +123,7 @@ Like MatrixOperator, but stores a Factorization instead.
 
 Supports left division and `ldiv!` when applied to an array.
 """
-struct InvertibleOperator{T,FType} <: AbstractSciMLLinearOperator{T}
+struct InvertibleOperator{T,FType} <: AbstractSciMLOperator{T}
     F::FType
 
     function InvertibleOperator(F)
@@ -131,7 +133,7 @@ struct InvertibleOperator{T,FType} <: AbstractSciMLLinearOperator{T}
 end
 
 # constructor
-function LinearAlgebra.factorize(L::AbstractSciMLLinearOperator)
+function LinearAlgebra.factorize(L::AbstractSciMLOperator)
     fact = factorize(convert(AbstractMatrix, L))
     InvertibleOperator(fact)
 end
@@ -172,6 +174,7 @@ LinearAlgebra.opnorm(L::InvertibleOperator{T}, p=2) where{T} = one(T) / opnorm(L
 LinearAlgebra.issuccess(L::InvertibleOperator) = issuccess(L.F)
 
 getops(L::InvertibleOperator) = (L.F,)
+islinear(L::InvertibleOperator) = islinear(L.F)
 
 @forward InvertibleOperator.F (
                                # LinearAlgebra

src/scalar.jl

@@ -26,7 +26,8 @@ Base.adjoint(α::AbstractSciMLScalarOperator) = conj(α)
 Base.transpose(α::AbstractSciMLScalarOperator) = α
 
 has_mul!(::AbstractSciMLScalarOperator) = true
-issquare(L::AbstractSciMLScalarOperator) = true
+issquare(::AbstractSciMLScalarOperator) = true
+islinear(::AbstractSciMLScalarOperator) = true
 has_adjoint(::AbstractSciMLScalarOperator) = true
 
 Base.:*(α::AbstractSciMLScalarOperator, u::AbstractVecOrMat) = convert(Number, α) * u

test/basic.jl

@@ -27,6 +27,8 @@ K = 12
     β = rand()
     Id = IdentityOperator{N}()
 
+    @test issquare(Id)
+    @test islinear(Id)
     @test IdentityOperator(u) isa IdentityOperator{N}
     @test one(A) isa IdentityOperator{N}
     @test convert(AbstractMatrix, Id) == Matrix(I, N, N)
@@ -61,6 +63,8 @@ end
     β = rand()
     Z = NullOperator{N}()
 
+    @test issquare(Z)
+    @test islinear(Z)
     @test NullOperator(u) isa NullOperator{N}
     @test zero(A) isa NullOperator{N}
     @test convert(AbstractMatrix, Z) == zeros(size(Z))
@@ -99,6 +103,8 @@ end
     op = ScaledOperator(α, MatrixOperator(A))
 
     @test op isa ScaledOperator
+    @test issquare(op)
+    @test islinear(op)
 
     @test α * A * u ≈ op * u
     @test (β * op) * u ≈ β * α * A * u
@@ -164,10 +170,14 @@ end
 
     @test op isa ComposedOperator
     @test *(op.ops...) isa ComposedOperator
+    @test issquare(op)
+    @test islinear(op)
 
     opF = factorize(op)
 
     @test opF isa ComposedOperator
+    @test issquare(opF)
+    @test islinear(opF)
 
     @test ABCmulu ≈ op * u
     @test ABCdivu ≈ op \ u ≈ opF \ u
@@ -202,10 +212,14 @@ end
 
 @testset "Adjoint, Transpose" begin
 
-    for (op, LType, VType) in (
-                               (adjoint,   AdjointOperator,    AbstractAdjointVecOrMat   ),
-                               (transpose, TransposedOperator, AbstractTransposedVecOrMat),
-                              )
+    for (op,
+         LType,
+         VType
+        ) in (
+              (adjoint,   AdjointOperator,    AbstractAdjointVecOrMat   ),
+              (transpose, TransposedOperator, AbstractTransposedVecOrMat),
+             )
+
         A = rand(N,N)
         D = Bidiagonal(rand(N,N), :L)
         u = rand(N,K)
@@ -215,7 +229,13 @@ end
         b = rand()
 
         At = op(A)
-        Dt = op(A)
+        Dt = op(D)
+
+        @test issquare(At)
+        @test issquare(Dt)
+
+        @test islinear(At)
+        @test islinear(Dt)
 
         AA = MatrixOperator(A)
         DD = MatrixOperator(D)
@@ -247,6 +267,9 @@ end
 
     Di = cache_operator(Di, u)
 
+    @test issquare(Di)
+    @test islinear(Di)
+
     @test Di * u ≈ u ./ s
     v=rand(N); @test mul!(v, Di, u) ≈ u ./ s
     v=rand(N); w=copy(v); @test mul!(v, Di, u, α, β) ≈ α *(u ./ s) + β*w