From 2782977ada14a6e2ab9443ef03fab34b928c8132 Mon Sep 17 00:00:00 2001
From: Gaurav Arya <gauravarya272@gmail.com>
Date: Fri, 21 Jul 2023 00:22:45 -0400
Subject: [PATCH 01/11] Progress towards ConcretizedOperator

---
 src/basic.jl | 97 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 97 insertions(+)

diff --git a/src/basic.jl b/src/basic.jl
index 3e5dd74f..1dd5351d 100644
--- a/src/basic.jl
+++ b/src/basic.jl
@@ -819,4 +819,101 @@ function LinearAlgebra.ldiv!(L::InvertedOperator, u::AbstractVecOrMat)
     copy!(L.cache, u)
     mul!(u, L.L, L.cache)
 end
+
+struct ConcretizedOperator{T, LType, AType} <: AbstractSciMLOperator{T}
+    L::LType
+    A::AType
+    function ConcretizedOperator(L::AbstractSciMLOperator{T}, A::AbstractMatrix) where {T}
+        new{T,typeof(L),typeof(A)}(L, A)
+    end
+end
+
+
+"""
+    ConcretizedOperator(L::AbstractSciMLOperator)
+
+Concretization of a SciMLOperator `L`, with a concrete backing `A = concretize(L)` that is used
+for all linear algebra operations. Unlike `A` itself, a concretized operator correctly supports
+updates to the operator state, by first updating `L` and then updating `A` accordingly.
+"""
+function ConcretizedOperator(L::AbstractSciMLOperator{T}) where {T}
+    ConcretizedOperator(L, convert(AbstractMatrix, L))
+end
+
+Base.convert(::Type{AbstractMatrix}, L::ConcretizedOperator) = L.A
+
+function Base.show(io::IO, L::ConcretizedOperator)
+    print(io, "ConcretizedOperator(")
+    show(io, L.L)
+    print(io, ")")
+end
+Base.size(L::ConcretizedOperator) = size(L.A)
+function Base.resize!(L::ConcretizedOperator, n::Integer)
+
+    resize!(L.L, n)
+    resize!(L.cache, n)
+    # TODO: these next two lines seem dangerous... in which cases do they make sense?
+    resize!(L.A, n)
+    concretize!(L.A, L.L)
+end
+
+function update_coefficients(L::ConcretizedOperator, u, p, t)
+    @set! L.L = update_coefficients(L.L, u, p, t)
+    @set! L.A = concretize(L.L)
+end
+
+function update_coefficients!(L::ConcretizedOperator, u, p, t)
+    for op in getops(L)
+        update_coefficients!(op, u, p, t; kwargs...)
+    end
+    concretize!(L.A, L.L) # TODO: this needs to be supported
+end
+
+getops(L::ConcretizedOperator) = (L.L,)
+islinear(L::ConcretizedOperator) = islinear(L.L)
+isconvertible(::ConcretizedOperator) = true 
+
+for op in (
+           :adjoint,
+           :transpose,
+           :conj,
+          )
+    @eval Base.$op(L::ConcretizedOperator) = ConcretizedOperator($op(L.L), $op(L.A))
+end
+
+@forward ConcretizedOperator.L (
+                             # LinearAlgebra
+                             LinearAlgebra.issymmetric,
+                             LinearAlgebra.ishermitian,
+                             LinearAlgebra.isposdef,
+                             LinearAlgebra.opnorm,
+
+                             # SciML
+                             isconstant,
+                             has_adjoint,
+                             has_mul,
+                             has_mul!,
+                             has_ldiv,
+                             has_ldiv!
+                            )
+
+Base.:*(L::ConcretizedOperator, u::AbstractVecOrMat) = L.A * u
+Base.:\(L::ConcretizedOperator, u::AbstractVecOrMat) = L.A \ u
+
+function LinearAlgebra.mul!(v::AbstractVecOrMat, L::ConcretizedOperator, u::AbstractVecOrMat)
+    mul!(v, L.A, u)
+end
+
+function LinearAlgebra.mul!(v::AbstractVecOrMat, L::ConcretizedOperator, u::AbstractVecOrMat, α, β)
+    mul!(v, L.A, u, α, β)
+end
+
+function LinearAlgebra.ldiv!(v::AbstractVecOrMat, L::ConcretizedOperator, u::AbstractVecOrMat)
+    ldiv!(v, L.A, u)
+end
+
+function LinearAlgebra.ldiv!(L::ConcretizedOperator, u::AbstractVecOrMat)
+    ldiv!(L.A, u)
+end
+
 #

From 0afc6efd5c76903499f84797aea686d905744e8d Mon Sep 17 00:00:00 2001
From: Gaurav Arya <gauravarya272@gmail.com>
Date: Fri, 21 Jul 2023 00:26:54 -0400
Subject: [PATCH 02/11] Add concretize! fallback

---
 src/basic.jl     |  2 +-
 src/interface.jl | 13 +++++++++++++
 2 files changed, 14 insertions(+), 1 deletion(-)

diff --git a/src/basic.jl b/src/basic.jl
index 1dd5351d..d5711ea0 100644
--- a/src/basic.jl
+++ b/src/basic.jl
@@ -866,7 +866,7 @@ function update_coefficients!(L::ConcretizedOperator, u, p, t)
     for op in getops(L)
         update_coefficients!(op, u, p, t; kwargs...)
     end
-    concretize!(L.A, L.L) # TODO: this needs to be supported
+    concretize!(L.A, L.L) # TODO: this needs to be supported. Also, problematic if L.A is scalar, should we only support matrix?
 end
 
 getops(L::ConcretizedOperator) = (L.L,)
diff --git a/src/interface.jl b/src/interface.jl
index 8437f0f3..364c0cdb 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -311,6 +311,19 @@ concretize(L::Union{
                    }
           ) = convert(Number, L)
 
+function concretize!(A, L::AbstractMatrix)
+    A .= L 
+end
+
+function concretize!(A, L::Union{Factorization, AbstractSciMLOperator}) 
+    @warn """using concretize-based fallback for concretize!"""
+    A .= concretize(L)
+end
+
+"""
+
+"""
+
 """
 $SIGNATURES
 

From b503519d9a256a76479f49879ebe5e1eaf37c366 Mon Sep 17 00:00:00 2001
From: Gaurav Arya <gauravarya272@gmail.com>
Date: Fri, 21 Jul 2023 00:28:37 -0400
Subject: [PATCH 03/11] Add concretization for MatrixOperator

---
 src/matrix.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/matrix.jl b/src/matrix.jl
index 30687f2c..d5156274 100644
--- a/src/matrix.jl
+++ b/src/matrix.jl
@@ -159,6 +159,10 @@ function update_coefficients!(L::MatrixOperator, u, p, t; kwargs...)
     L.update_func!(L.A, u, p, t; kwargs...)
 end
 
+function concretize!(A, L:::MatrixOperator)
+    return concretize!(A, L.A)
+end
+
 SparseArrays.sparse(L::MatrixOperator) = sparse(L.A)
 SparseArrays.issparse(L::MatrixOperator) = issparse(L.A)
 

From 3157835702ce66c476b6d81ea7abb56b72308957 Mon Sep 17 00:00:00 2001
From: Gaurav Arya <gauravarya272@gmail.com>
Date: Fri, 21 Jul 2023 00:31:22 -0400
Subject: [PATCH 04/11] Fix typos

---
 src/interface.jl | 4 ----
 src/matrix.jl    | 2 +-
 2 files changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/interface.jl b/src/interface.jl
index 364c0cdb..82913df5 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -320,10 +320,6 @@ function concretize!(A, L::Union{Factorization, AbstractSciMLOperator})
     A .= concretize(L)
 end
 
-"""
-
-"""
-
 """
 $SIGNATURES
 
diff --git a/src/matrix.jl b/src/matrix.jl
index d5156274..5988297d 100644
--- a/src/matrix.jl
+++ b/src/matrix.jl
@@ -159,7 +159,7 @@ function update_coefficients!(L::MatrixOperator, u, p, t; kwargs...)
     L.update_func!(L.A, u, p, t; kwargs...)
 end
 
-function concretize!(A, L:::MatrixOperator)
+function concretize!(A, L::MatrixOperator)
     return concretize!(A, L.A)
 end
 

From e8c120d43075c30475501cbd6e01e8cc70fb18df Mon Sep 17 00:00:00 2001
From: Gaurav Arya <gauravarya272@gmail.com>
Date: Fri, 21 Jul 2023 00:36:59 -0400
Subject: [PATCH 05/11] tweak

---
 src/basic.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/basic.jl b/src/basic.jl
index d5711ea0..062e3386 100644
--- a/src/basic.jl
+++ b/src/basic.jl
@@ -849,7 +849,6 @@ function Base.show(io::IO, L::ConcretizedOperator)
 end
 Base.size(L::ConcretizedOperator) = size(L.A)
 function Base.resize!(L::ConcretizedOperator, n::Integer)
-
     resize!(L.L, n)
     resize!(L.cache, n)
     # TODO: these next two lines seem dangerous... in which cases do they make sense?

From 7bba2449a63ac51909bc894278e4a9a2c611be70 Mon Sep 17 00:00:00 2001
From: Gaurav Arya <gauravarya272@gmail.com>
Date: Fri, 21 Jul 2023 00:44:41 -0400
Subject: [PATCH 06/11] Allow scaling in concretize! interface

---
 src/interface.jl | 14 ++++++++++----
 1 file changed, 10 insertions(+), 4 deletions(-)

diff --git a/src/interface.jl b/src/interface.jl
index 82913df5..2c4d1e9b 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -311,13 +311,19 @@ concretize(L::Union{
                    }
           ) = convert(Number, L)
 
-function concretize!(A, L::AbstractMatrix)
-    A .= L 
+function concretize!(A, L)
+    T = eltype(L)
+    concretize!(A, L, one(T), zero(T))
 end
 
-function concretize!(A, L::Union{Factorization, AbstractSciMLOperator}) 
+function concretize!(A, L::AbstractMatrix, α, β)
+    A .= α .* L .+ β .* A
+end
+
+function concretize!(A, L::Union{Factorization, AbstractSciMLOperator}, α, β) 
     @warn """using concretize-based fallback for concretize!"""
-    A .= concretize(L)
+    # TODO: could also use a mul! based fallback on the unit vectors
+    concretize!(A, concertize(L), α, β)
 end
 
 """

From 173de90f728133b5c249d5ad59915d7a5df874e1 Mon Sep 17 00:00:00 2001
From: Gaurav Arya <gauravarya272@gmail.com>
Date: Fri, 21 Jul 2023 00:50:11 -0400
Subject: [PATCH 07/11] Add concretize! for scaled and added operator

---
 src/basic.jl | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/src/basic.jl b/src/basic.jl
index 062e3386..e577acd5 100644
--- a/src/basic.jl
+++ b/src/basic.jl
@@ -259,6 +259,10 @@ has_mul!(L::ScaledOperator) = has_mul!(L.L)
 has_ldiv(L::ScaledOperator) = has_ldiv(L.L) & !iszero(L.λ)
 has_ldiv!(L::ScaledOperator) = has_ldiv!(L.L) & !iszero(L.λ)
 
+function concretize!(A, L::ScaledOperator{T}, α, β) where {T}
+    concretize!(A, L.L, L.λ * α, β)
+end
+
 function cache_internals(L::ScaledOperator, u::AbstractVecOrMat)
     @set! L.L = cache_operator(L.L, u)
     @set! L.λ = cache_operator(L.λ, u)
@@ -420,6 +424,13 @@ function update_coefficients(L::AddedOperator, u, p, t)
     @set! L.ops = ops
 end
 
+function concretize!(A, L::AddedOperator{T}, α, β) where {T}
+    A .*= β
+    for op in L.ops
+        concretize!(A, op, α, one(T))
+    end
+end
+
 getops(L::AddedOperator) = L.ops
 islinear(L::AddedOperator) = all(islinear, getops(L))
 Base.iszero(L::AddedOperator) = all(iszero, getops(L))

From 1fb11e2391833487d106eae3369c2b233bacad21 Mon Sep 17 00:00:00 2001
From: Gaurav Arya <gauravarya272@gmail.com>
Date: Fri, 21 Jul 2023 00:59:40 -0400
Subject: [PATCH 08/11] Add basic test

---
 src/SciMLOperators.jl |  3 ++-
 src/basic.jl          |  6 +++---
 src/interface.jl      |  2 +-
 test/basic.jl         | 10 ++++++++++
 4 files changed, 16 insertions(+), 5 deletions(-)

diff --git a/src/SciMLOperators.jl b/src/SciMLOperators.jl
index 4409e920..5c7aa847 100644
--- a/src/SciMLOperators.jl
+++ b/src/SciMLOperators.jl
@@ -87,7 +87,8 @@ export
        AffineOperator,
        AddVector,
        FunctionOperator,
-       TensorProductOperator
+       TensorProductOperator,
+       ConcretizedOperator
 
 export update_coefficients!,
        update_coefficients,
diff --git a/src/basic.jl b/src/basic.jl
index e577acd5..639877e1 100644
--- a/src/basic.jl
+++ b/src/basic.jl
@@ -867,12 +867,12 @@ function Base.resize!(L::ConcretizedOperator, n::Integer)
     concretize!(L.A, L.L)
 end
 
-function update_coefficients(L::ConcretizedOperator, u, p, t)
-    @set! L.L = update_coefficients(L.L, u, p, t)
+function update_coefficients(L::ConcretizedOperator, u, p, t; kwargs...)
+    @set! L.L = update_coefficients(L.L, u, p, t; kwargs...)
     @set! L.A = concretize(L.L)
 end
 
-function update_coefficients!(L::ConcretizedOperator, u, p, t)
+function update_coefficients!(L::ConcretizedOperator, u, p, t; kwargs...)
     for op in getops(L)
         update_coefficients!(op, u, p, t; kwargs...)
     end
diff --git a/src/interface.jl b/src/interface.jl
index 2c4d1e9b..75f50a91 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -323,7 +323,7 @@ end
 function concretize!(A, L::Union{Factorization, AbstractSciMLOperator}, α, β) 
     @warn """using concretize-based fallback for concretize!"""
     # TODO: could also use a mul! based fallback on the unit vectors
-    concretize!(A, concertize(L), α, β)
+    concretize!(A, concretize(L), α, β)
 end
 
 """
diff --git a/test/basic.jl b/test/basic.jl
index 89e5d919..fddca55c 100644
--- a/test/basic.jl
+++ b/test/basic.jl
@@ -326,4 +326,14 @@ end
     v=rand(N); @test ldiv!(v, Di, u) ≈ u .* s
     v=copy(u); @test ldiv!(Di, u) ≈ v .* s
 end
+
+@testset "ConcretizedOperator" begin
+    A = rand(2, 2); B = rand(2, 2);
+    L = MatrixOperator(A; update_func=(u,p,t)->t * A) + MatrixOperator(B; update_func=(u,p,t)->t * B)
+    C = ConcretizedOperator(L)
+    v = rand(2)
+    C * v ≈ L * v
+    update_coefficients!(C, nothing, nothing, 2.0)
+    C * v ≈ L * v
+end
 #

From 958e00afc03c9a51904a672e6fb423a006654db4 Mon Sep 17 00:00:00 2001
From: Gaurav Arya <gauravarya272@gmail.com>
Date: Fri, 21 Jul 2023 01:03:46 -0400
Subject: [PATCH 09/11] Fix MatrixOperator concretize!

---
 src/interface.jl | 2 +-
 src/matrix.jl    | 4 ++--
 test/basic.jl    | 2 ++
 3 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/interface.jl b/src/interface.jl
index 75f50a91..778e25f1 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -321,7 +321,7 @@ function concretize!(A, L::AbstractMatrix, α, β)
 end
 
 function concretize!(A, L::Union{Factorization, AbstractSciMLOperator}, α, β) 
-    @warn """using concretize-based fallback for concretize!"""
+    @warn """using concretize-based fallback for concretize! for $(typeof(L))"""
     # TODO: could also use a mul! based fallback on the unit vectors
     concretize!(A, concretize(L), α, β)
 end
diff --git a/src/matrix.jl b/src/matrix.jl
index 5988297d..2a1daf3b 100644
--- a/src/matrix.jl
+++ b/src/matrix.jl
@@ -159,8 +159,8 @@ function update_coefficients!(L::MatrixOperator, u, p, t; kwargs...)
     L.update_func!(L.A, u, p, t; kwargs...)
 end
 
-function concretize!(A, L::MatrixOperator)
-    return concretize!(A, L.A)
+function concretize!(A, L::MatrixOperator, α, β)
+    return concretize!(A, L.A, α, β)
 end
 
 SparseArrays.sparse(L::MatrixOperator) = sparse(L.A)
diff --git a/test/basic.jl b/test/basic.jl
index fddca55c..0ed82c12 100644
--- a/test/basic.jl
+++ b/test/basic.jl
@@ -333,7 +333,9 @@ end
     C = ConcretizedOperator(L)
     v = rand(2)
     C * v ≈ L * v
+    @test C.A ≈ convert(AbstractMatrix, L)
     update_coefficients!(C, nothing, nothing, 2.0)
     C * v ≈ L * v
+    @test C.A ≈ convert(AbstractMatrix, L)
 end
 #

From 36d2ef266f36eb0dc40b8c25e90618f9d5f6314c Mon Sep 17 00:00:00 2001
From: Gaurav Arya <gauravarya272@gmail.com>
Date: Fri, 21 Jul 2023 13:09:35 -0400
Subject: [PATCH 10/11] Fix ScaledOperator concretize!

---
 src/basic.jl  | 2 +-
 test/basic.jl | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/basic.jl b/src/basic.jl
index 639877e1..2d1575e5 100644
--- a/src/basic.jl
+++ b/src/basic.jl
@@ -260,7 +260,7 @@ has_ldiv(L::ScaledOperator) = has_ldiv(L.L) & !iszero(L.λ)
 has_ldiv!(L::ScaledOperator) = has_ldiv!(L.L) & !iszero(L.λ)
 
 function concretize!(A, L::ScaledOperator{T}, α, β) where {T}
-    concretize!(A, L.L, L.λ * α, β)
+    concretize!(A, L.L, convert(Number, L.λ) * α, β)
 end
 
 function cache_internals(L::ScaledOperator, u::AbstractVecOrMat)
diff --git a/test/basic.jl b/test/basic.jl
index 0ed82c12..c6110b43 100644
--- a/test/basic.jl
+++ b/test/basic.jl
@@ -329,7 +329,7 @@ end
 
 @testset "ConcretizedOperator" begin
     A = rand(2, 2); B = rand(2, 2);
-    L = MatrixOperator(A; update_func=(u,p,t)->t * A) + MatrixOperator(B; update_func=(u,p,t)->t * B)
+    L = MatrixOperator(A; update_func=(u,p,t)->t * A) + 2 * MatrixOperator(B; update_func=(u,p,t)->t * B)
     C = ConcretizedOperator(L)
     v = rand(2)
     C * v ≈ L * v

From 508da8c19e6292123b83f8e47bef3ee2e14cf9bf Mon Sep 17 00:00:00 2001
From: Gaurav Arya <gauravarya272@gmail.com>
Date: Fri, 21 Jul 2023 15:07:11 -0400
Subject: [PATCH 11/11] Add concretize! for IdentityOperator

---
 src/basic.jl  | 8 ++++++++
 test/basic.jl | 2 +-
 2 files changed, 9 insertions(+), 1 deletion(-)

diff --git a/src/basic.jl b/src/basic.jl
index 2d1575e5..8bd5067a 100644
--- a/src/basic.jl
+++ b/src/basic.jl
@@ -38,6 +38,14 @@ has_mul!(::IdentityOperator) = true
 has_ldiv(::IdentityOperator) = true
 has_ldiv!(::IdentityOperator) = true
 
+function concretize!(A, L::IdentityOperator, α, β) 
+    @assert size(A) == (L.len, L.len)
+    A .*= β
+    for i in 1:L.len
+        A[i, i] += α
+    end
+end
+
 # opeator application
 for op in (
            :*, :\,
diff --git a/test/basic.jl b/test/basic.jl
index c6110b43..38a952c6 100644
--- a/test/basic.jl
+++ b/test/basic.jl
@@ -329,7 +329,7 @@ end
 
 @testset "ConcretizedOperator" begin
     A = rand(2, 2); B = rand(2, 2);
-    L = MatrixOperator(A; update_func=(u,p,t)->t * A) + 2 * MatrixOperator(B; update_func=(u,p,t)->t * B)
+    L = MatrixOperator(A; update_func=(u,p,t)->t * A) + 2 * MatrixOperator(B; update_func=(u,p,t)->t * B) + I
     C = ConcretizedOperator(L)
     v = rand(2)
     C * v ≈ L * v