fixes

Author: oscardssmith

src/woperator.jl

@@ -1,4 +1,4 @@
-abstract type AbstractWOperator{T} <: AbstractSciMLOperator{T} 
+abstract type AbstractWOperator{T} <: AbstractSciMLOperator{T} end
 
 struct StaticWOperator{isinv, T, F} <: AbstractWOperator{T}
     W::T
@@ -8,9 +8,7 @@ struct StaticWOperator{isinv, T, F} <: AbstractWOperator{T}
         isinv = n <= 7 # cutoff where LU has noticable overhead
 
         F = if isinv && callinv
-            # this should be in ArrayInterface but can't be for silly reasons
-            # doing to how StaticArrays and StaticArraysCore are split up
-            StaticArrays.LU(LowerTriangular(W), UpperTriangular(W), SVector{n}(1:n))
+            ArrayInterface.lu_instance(W)
         else
             lu(W, check = false)
         end
@@ -36,18 +34,10 @@ W = \\frac{1}{\\gamma}MM - J
 ```
 
 where `MM` is the mass matrix (a regular `AbstractMatrix` or a `UniformScaling`),
-`γ` is a real number proportional to the time step, and `J` is the Jacobian
-operator (must be a `AbstractSciMLOperator`). A `WOperator` can also be
-constructed using a `*DEFunction` directly as
-
-    WOperator(f,gamma)
-
-`f` needs to have a jacobian and `jac_prototype`, but the prototype does not need
-to be a diffeq operator --- it will automatically be converted to one.
+`γ` is a real number, and `J` is the Jacobian operator (must be a `AbstractSciMLOperator`).
 
 `WOperator` supports lazy `*` and `mul!` operations, the latter utilizing an
 internal cache (can be specified in the constructor; default to regular `Vector`).
-It supports all of `AbstractSciMLOperator`'s interface.
 """
 mutable struct WOperator{IIP, T,
     MType,
@@ -64,50 +54,11 @@ mutable struct WOperator{IIP, T,
     jacvec::JV
 
     function WOperator{IIP}(mass_matrix, gamma, J, u, jacvec = nothing) where {IIP}
-        # TODO: there is definitely a missing interface.
-        # Tentative interface: `has_concrete` and `concertize(A)`
-        if J isa Union{Number, ScalarOperator}
-            _concrete_form = -mass_matrix / gamma + convert(Number, J)
-            _func_cache = nothing
-        else
-            AJ = J isa MatrixOperator ? convert(AbstractMatrix, J) : J
-            if AJ isa AbstractMatrix
-                mm = mass_matrix isa MatrixOperator ?
-                     convert(AbstractMatrix, mass_matrix) : mass_matrix
-                if is_sparse(AJ)
-
-                    # If gamma is zero, then it's just an initialization and we want to make sure
-                    # we get the right sparsity pattern. If gamma is not zero, then it's a case where
-                    # a new W is created (as part of an out-of-place solve) and thus the actual
-                    # values actually matter!
-                    #
-                    # Constant operators never refactorize so always use the correct values there
-                    # as well
-                    if gamma == 0 && !(J isa MatrixOperator && isconstant(J))
-                        # Workaround https://github.com/JuliaSparse/SparseArrays.jl/issues/190
-                        # Hopefully `rand()` does not match any value in the array (prob ~ 0, with a check)
-                        # Then `one` is required since gamma is zero
-                        # Otherwise this will not pick up the union sparsity pattern
-                        # But instead drop the runtime zeros (i.e. all values) of the AJ pattern!
-                        AJn = nonzeros(AJ)
-                        x = rand()
-                        @assert all(!isequal(x), AJn)
-
-                        fill!(AJn, rand())
-                        _concrete_form = -mm / one(gamma) + AJ
-                        fill!(_concrete_form, false) # safety measure, throw singular error if not filled
-                    else
-                        _concrete_form = -mm / gamma + AJ
-                    end
-                else
-                    _concrete_form = -mm / gamma + AJ
-                end
-
-            else
-                _concrete_form = nothing
-            end
-            _func_cache = zero(u)
-        end
+        AJ = J isa MatrixOperator ? convert(AbstractMatrix, J) : J
+        mm = mass_matrix isa MatrixOperator ?
+             convert(AbstractMatrix, mass_matrix) : mass_matrix
+        _concrete_form = -mm / gamma + AJ
+        _func_cache = zero(u)
         T = eltype(_concrete_form)
         MType = typeof(mass_matrix)
         GType = typeof(gamma)
@@ -116,52 +67,42 @@ mutable struct WOperator{IIP, T,
         C = typeof(_concrete_form)
         JV = typeof(jacvec)
         return new{IIP, T, MType, GType, JType, F, C, JV}(mass_matrix, gamma, J,
-            _func_cache, _concrete_form,
-            jacvec)
+            _func_cache, _concrete_form, jacvec)
     end
 
     function Base.copy(W::WOperator{IIP, T, MType, GType, JType, F, C, JV}) where {IIP, T, MType, GType, JType, F, C, JV}
         return new{IIP, T, MType, GType, JType, F, C, JV}(
             W.mass_matrix, 
             W.gamma, 
             W.J, 
-            W._func_cache === nothing ? nothing : copy(W._func_cache),
-            W._concrete_form === nothing ? nothing : copy(W._concrete_form),
+            copy(W._func_cache),
+            copy(W._concrete_form),
             W.jacvec
         )
     end
 end
 Base.eltype(W::WOperator) = eltype(W.J)
 
-# In WOperator update_coefficients!, accept both missing u/p/t and missing dtgamma and don't update them in that case.
+# In WOperator update_coefficients!, accept both missing u/p/t and missing gamma and don't update them in that case.
 # This helps support partial updating logic used with Newton solvers.
 function update_coefficients!(W::WOperator,
         u = nothing,
         p = nothing,
         t = nothing;
-        dtgamma = nothing)
+        gamma = nothing)
     if (u !== nothing) && (p !== nothing) && (t !== nothing)
         update_coefficients!(W.J, u, p, t)
         update_coefficients!(W.mass_matrix, u, p, t)
         !isnothing(W.jacvec) && update_coefficients!(W.jacvec, u, p, t)
     end
-    dtgamma !== nothing && (W.gamma = dtgamma)
+    gamma !== nothing && (W.gamma = gamma)
     W
 end
 
-function update_coefficients!(J::UJacobianWrapper, u, p, t)
-    J.p = p
-    J.t = t
-end
-
 function Base.convert(::Type{AbstractMatrix}, W::WOperator{IIP}) where {IIP}
     if !IIP
         # Allocating
         W._concrete_form = -W.mass_matrix / W.gamma + convert(AbstractMatrix, W.J)
-    else
-        # Non-allocating already updated
-        #_W = W._concrete_form
-        #jacobian2W!(_W, W.mass_matrix, W.gamma, W.J)
     end
     return W._concrete_form
 end
@@ -202,19 +143,18 @@ function LinearAlgebra.mul!(Y::AbstractVecOrMat, W::WOperator, B::AbstractVecOrM
     # Compute mass_matrix * B
     if isa(W.mass_matrix, UniformScaling)
         a = -W.mass_matrix.λ / W.gamma
-        @.. broadcast=false Y=a*B
+        @. Y=a*B
     else
         mul!(_vec(Y), W.mass_matrix, _vec(B))
         lmul!(-inv(W.gamma), Y)
     end
     # Compute J * B and add
-    if W.jacvec !== nothing
-        mul!(_vec(W._func_cache), W.jacvec, _vec(B))
-    else
+    if isnothing(W.jacvec)
         mul!(_vec(W._func_cache), W.J, _vec(B))
+    else
+        mul!(_vec(W._func_cache), W.jacvec, _vec(B))
     end
     _vec(Y) .+= _vec(W._func_cache)
 end
 
-
 has_concretization(::AbstractWOperator) = true

test/WOperator.jl

@@ -0,0 +1,14 @@
+using SciMLOperators, LinearAlgebra
+using Random
+using Test
+
+Random.seed!(0)
+@testset "WOperator" begin
+    J = rand(12, 12)
+    u = rand(12)
+    M = I(12)
+    gamma = 1/123
+    W = WOperator{true}(M, gamma, J, u)
+    
+    @test convert(AbstractMatrix, W) ≈ J - M/gamma
+end