Allow matrix-backed operator Jacobians in LHL defaults

docs/src/tutorials/shifted_systems.md

@@ -189,9 +189,9 @@ uc ≈ (J - LA.I / γc) \ bc
 ## When it is chosen automatically
 
 Wrapping a matrix in a `WOperator` is itself the statement that the shift will move while
-`J` stays put, so `defaultalg` treats it as one: for a `WOperator` with a dense Jacobian, a
-scalar multiple of `I` as the mass matrix, and `n ≥ LinearSolve.LHL_DEFAULT_MIN_SIZE`, no
-algorithm argument is needed.
+`J` stays put, so `defaultalg` treats it as one: for a `WOperator` with a dense matrix or
+matrix-backed operator as its Jacobian, a scalar multiple of `I` as the mass matrix, and
+`n ≥ LinearSolve.LHL_DEFAULT_MIN_SIZE`, no algorithm argument is needed.
 
 ```@example shifted
 LS.defaultalg(W, b, LS.OperatorAssumptions(true))

src/default.jl

@@ -1,4 +1,6 @@
 needs_concrete_A(alg::DefaultLinearSolver) = true
+# Jacobian staleness belongs to the wrapper, so copying it would disconnect the signal.
+default_alias_A(::DefaultLinearSolver, ::WOperator, ::Any) = true
 
 # Every algorithm the default can dispatch to either ignores the tolerances
 # (the factorizations) or reads `cache.abstol`/`cache.reltol` at solve time (the
@@ -233,11 +235,11 @@ const LHL_DEFAULT_MIN_SIZE = 32
 # but left uninitialized the solve fails, and if it is initialized but never selected the
 # buffers are wasted. Both ask here.
 function _lhl_defaultable(A::WOperator, assump::OperatorAssumptions)
-    # A `MatrixOperator` is updated in place while its identity and `jac_stale` remain
-    # unchanged, so the reduction cannot detect that its contents moved.
     (assump.issq && _lhl_scalar_massmatrix(A.mass_matrix)) || return false
-    A.J isa DenseMatrix && return size(A, 1) >= LHL_DEFAULT_MIN_SIZE
-    return issparsematrixcsc(A.J)
+    (A.J isa AbstractMatrix || A.J isa MatrixOperator) || return false
+    J = _lhl_jacobian(A)
+    J isa DenseMatrix && return size(A, 1) >= LHL_DEFAULT_MIN_SIZE
+    return issparsematrixcsc(J)
 end
 _lhl_defaultable(A, assump) = false
 

src/lhl.jl

@@ -19,9 +19,9 @@ fixed for many steps.
 
 Give it the system matrix unassembled, as a `SciMLOperators.WOperator` — the split
 `J - M/γ` an implicit solver already builds — and move the shift with
-[`update_gamma!`](@ref).  A `WOperator` whose Jacobian is a dense matrix is also what
-`defaultalg` selects this algorithm for, at sizes where the reduction pays.  The mass
-matrix must be a multiple of `I`.
+[`update_gamma!`](@ref).  A `WOperator` whose Jacobian is a dense matrix or matrix-backed
+operator is also what `defaultalg` selects this algorithm for, at sizes where the
+reduction pays.  The mass matrix must be a multiple of `I`.
 
 Handed an ordinary matrix `A`, it solves `Ax = b` as `Z H⁻¹ Z⁻¹ b`; that works, but it is
 strictly worse than an LU and the algorithm has no reason to be chosen.

test/Core/lhl.jl

@@ -393,11 +393,7 @@ end
     end
 end
 
-@testset "an operator Jacobian is never claimed by the reduction" begin
-    # `update_coefficients!` moves a `MatrixOperator`'s numbers in place, leaving both the
-    # object identity and `jac_stale` untouched. Nothing tells the reduction it went
-    # stale, so the split form must not be claimed at all for an operator `J` — it would
-    # answer with the previous Jacobian and raise nothing.
+@testset "a dense operator Jacobian uses LHL by default" begin
     n = 40
     γ = 0.1
     A = randn(MersenneTwister(21), n, n)
@@ -406,16 +402,18 @@ end
     W = WOperator{true}(I, γ, Aop, zeros(n))
     assump = LinearSolve.OperatorAssumptions(true)
 
-    @test !LinearSolve._lhl_defaultable(W, assump)
-    @test LinearSolve.defaultalg(W, b, assump) !=
+    @test LinearSolve._lhl_defaultable(W, assump)
+    @test LinearSolve.defaultalg(W, b, assump) ==
         LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LHLFactorization)
 
     cache = init(LinearProblem(W, b))
+    @test cache.A === W
+    @test cache.cacheval.LHLFactorization !== nothing
     @test solve!(cache).u ≈ (A - I / γ) \ b rtol = 1.0e-6
 
     Anew = randn(MersenneTwister(23), n, n)
     copyto!(Aop.A, Anew)
-    cache.isfresh = true
+    mark_jacobian_updated!(W)
     @test solve!(cache).u ≈ (Anew - I / γ) \ b rtol = 1.0e-6
 end