Merge pull request #179 from vpuri3/fwdbwd

Combine a forward operator with a backward operator for fast evaluaiton

Author: xtalax

src/func.jl

@@ -225,7 +225,8 @@ function iscached(L::FunctionOperator)
 end
 
 function cache_self(L::FunctionOperator, u::AbstractVecOrMat, v::AbstractVecOrMat)
-    !L.traits.ifcache && @warn "you are allocating cache for a FunctionOperator for which ifcache = false."
+    !L.traits.ifcache && @warn """Cache is being allocated for a FunctionOperator
+        created with kwarg ifcache = false."""
     @set! L.cache = zero.((u, v))
     L
 end

src/interface.jl

@@ -126,7 +126,7 @@ has_expmv!(L::AbstractSciMLOperator) = false # expmv!(v, L, t, u)
 has_expmv(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_mul!(L::AbstractSciMLOperator) = true # mul!(du, L, u)
 has_ldiv(L::AbstractSciMLOperator) = false # du = L\u
 has_ldiv!(L::AbstractSciMLOperator) = false # ldiv!(du, L, u)
 

src/matrix.jl

@@ -158,25 +158,28 @@ const AdjointFact = isdefined(LinearAlgebra, :AdjointFactorization) ? LinearAlge
 const TransposeFact = isdefined(LinearAlgebra, :TransposeFactorization) ? LinearAlgebra.TransposeFactorization : Transpose
 
 """
-    InvertibleOperator(F)
+    InvertibleOperator(L, F)
 
-Like MatrixOperator, but stores a Factorization instead.
-
-Supports left division and `ldiv!` when applied to an array.
+Stores an operator and its factorization (or inverse operator).
+Supports left division and `ldiv!` via `F`, and operator application
+via `L`.
 """
-struct InvertibleOperator{T,FType} <: AbstractSciMLOperator{T}
-    F::FType
+struct InvertibleOperator{T,LT,FT} <: AbstractSciMLOperator{T}
+    L::LT
+    F::FT
 
-    function InvertibleOperator(F)
+    function InvertibleOperator(L, F)
         @assert has_ldiv(F) | has_ldiv!(F) "$F is not invertible"
-        new{eltype(F),typeof(F)}(F)
+        T = promote_type(eltype(L), eltype(F))
+
+        new{T,typeof(L),typeof(F)}(L, F)
     end
 end
 
 # constructor
 function LinearAlgebra.factorize(L::AbstractSciMLOperator)
     fact = factorize(convert(AbstractMatrix, L))
-    InvertibleOperator(fact)
+    InvertibleOperator(L, fact)
 end
 
 for fact in (
@@ -190,40 +193,36 @@ for fact in (
             )
 
     @eval LinearAlgebra.$fact(L::AbstractSciMLOperator, args...) =
-        InvertibleOperator($fact(convert(AbstractMatrix, L), args...))
+        InvertibleOperator(L, $fact(convert(AbstractMatrix, L), args...))
     @eval LinearAlgebra.$fact(L::AbstractSciMLOperator; kwargs...) =
-        InvertibleOperator($fact(convert(AbstractMatrix, L); kwargs...))
+        InvertibleOperator(L, $fact(convert(AbstractMatrix, L); kwargs...))
 end
 
 function Base.convert(::Type{<:Factorization}, L::InvertibleOperator{T,<:Factorization}) where{T}
     L.F
 end
 
-Base.convert(::Type{AbstractMatrix}, L::InvertibleOperator) =
-    convert(AbstractMatrix, L.F)
-Base.convert(::Type{AbstractMatrix}, L::InvertibleOperator{<:Any,<:Union{Adjoint,AdjointFact}}) =
-    adjoint(convert(AbstractMatrix, adjoint(L.F)))
-Base.convert(::Type{AbstractMatrix}, L::InvertibleOperator{<:Any,<:Union{Transpose,TransposeFact}}) =
-    transpose(convert(AbstractMatrix, transpose(L.F)))
+Base.convert(::Type{AbstractMatrix}, L::InvertibleOperator) = convert(AbstractMatrix, L.L)
 
 # traits
-Base.size(L::InvertibleOperator) = size(L.F)
-Base.transpose(L::InvertibleOperator) = InvertibleOperator(transpose(L.F))
-Base.adjoint(L::InvertibleOperator) = InvertibleOperator(L.F')
-Base.conj(L::InvertibleOperator) = InvertibleOperator(conj(L.F))
-Base.resize!(L::InvertibleOperator, n::Integer) = (resize!(L.F, n); L)
+Base.size(L::InvertibleOperator) = size(L.L)
+Base.transpose(L::InvertibleOperator) = InvertibleOperator(transpose(L.L), transpose(L.F))
+Base.adjoint(L::InvertibleOperator) = InvertibleOperator(L.L', L.F')
+Base.conj(L::InvertibleOperator) = InvertibleOperator(conj(L.L), conj(L.F))
+Base.resize!(L::InvertibleOperator, n::Integer) = (resize!(L.L, n); resize!(L.F, n); L)
 LinearAlgebra.opnorm(L::InvertibleOperator{T}, p=2) where{T} = one(T) / opnorm(L.F)
 LinearAlgebra.issuccess(L::InvertibleOperator) = issuccess(L.F)
 
 function update_coefficients(L::InvertibleOperator, u, p, t)
+    @set! L.L = update_coefficients(L.L, u, p, t)
     @set! L.F = update_coefficients(L.F, u, p, t)
     L
 end
 
-getops(L::InvertibleOperator) = (L.F,)
-islinear(L::InvertibleOperator) = islinear(L.F)
+getops(L::InvertibleOperator) = (L.L, L.F,)
+islinear(L::InvertibleOperator) = islinear(L.L)
 
-@forward InvertibleOperator.F (
+@forward InvertibleOperator.L (
                                # LinearAlgebra
                                LinearAlgebra.issymmetric,
                                LinearAlgebra.ishermitian,
@@ -234,15 +233,16 @@ islinear(L::InvertibleOperator) = islinear(L.F)
                                has_adjoint,
                                has_mul,
                                has_mul!,
-                               has_ldiv,
-                               has_ldiv!,
                               )
 
+has_ldiv(L::InvertibleOperator) = has_mul(L.F)
+has_ldiv!(L::InvertibleOperator) = has_ldiv!(L.F)
+
 # operator application
-Base.:*(L::InvertibleOperator, x::AbstractVecOrMat) = L.F * x
+Base.:*(L::InvertibleOperator, x::AbstractVecOrMat) = L.L * x
 Base.:\(L::InvertibleOperator, x::AbstractVecOrMat) = L.F \ x
-LinearAlgebra.mul!(v::AbstractVecOrMat, L::InvertibleOperator, u::AbstractVecOrMat) = mul!(v, L.F, u)
-LinearAlgebra.mul!(v::AbstractVecOrMat, L::InvertibleOperator, u::AbstractVecOrMat,α, β) = mul!(v, L.F, u, α, β)
+LinearAlgebra.mul!(v::AbstractVecOrMat, L::InvertibleOperator, u::AbstractVecOrMat) = mul!(v, L.L, u)
+LinearAlgebra.mul!(v::AbstractVecOrMat, L::InvertibleOperator, u::AbstractVecOrMat,α, β) = mul!(v, L.L, u, α, β)
 LinearAlgebra.ldiv!(v::AbstractVecOrMat, L::InvertibleOperator, u::AbstractVecOrMat) = ldiv!(v, L.F, u)
 LinearAlgebra.ldiv!(L::InvertibleOperator, u::AbstractVecOrMat) = ldiv!(L.F, u)
 

test/basic.jl

@@ -222,14 +222,8 @@ end
     v=copy(u);   @test ldiv!(op, u)    ≈ (A * B * C) \ v
 
     # Test caching of composed operator when inner ops do not support Base.:*
-    # See issue #129
+    # ComposedOperator caching was modified in PR # 174
     inner_op = qr(MatrixOperator(rand(N, N)))
-    # We use the QR factorization of a non-square matrix, which does
-    # not support * as verified below.
-    @test !has_mul(inner_op)
-    @test has_ldiv(inner_op)
-    @test_throws MethodError inner_op * u
-    # We can now test that caching does not rely on matmul
     op = inner_op * factorize(MatrixOperator(rand(N, N)))
     @test !iscached(op)
     @test_nowarn op = cache_operator(op, rand(N))

test/func.jl

@@ -2,7 +2,7 @@
 using SciMLOperators, LinearAlgebra
 using Random
 
-using SciMLOperators: InvertibleOperator, ⊗
+using SciMLOperators: ⊗
 
 Random.seed!(0)
 N = 8

test/total.jl

@@ -137,7 +137,7 @@ end
                     (SciMLOperators.AdjointOperator(F), SciMLOperators.AdjointOperator),
                     (SciMLOperators.TransposedOperator(F), SciMLOperators.TransposedOperator),
                     (SciMLOperators.InvertedOperator(F), SciMLOperators.InvertedOperator),
-                    (SciMLOperators.InvertibleOperator(F), SciMLOperators.InvertibleOperator),
+                    (SciMLOperators.InvertibleOperator(F, F), SciMLOperators.InvertibleOperator),
                    )
 
         L = deepcopy(L)

test/zygote.jl

@@ -29,7 +29,7 @@ for (op_type, A) in
      (AffineOperator, AffineOperator(rand(N,N), rand(N,N), rand(N,K))),
      (ScaledOperator, rand() * MatrixOperator(rand(N,N))),
      (InvertedOperator, InvertedOperator(rand(N,N) |> MatrixOperator)),
-     (InvertibleOperator, InvertibleOperator(rand(N,N) |> MatrixOperator)),
+     (InvertibleOperator, InvertibleOperator(MatrixOperator(M), MatrixOperator(inv(M)))),
      (BatchedDiagonalOperator, DiagonalOperator(rand(N,K))),
      (AddedOperator, MatrixOperator(rand(N,N)) + MatrixOperator(rand(N,N))),
      (ComposedOperator, MatrixOperator(rand(N,N)) * MatrixOperator(rand(N,N))),