Extract eigenvalues

examples/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 KrylovPreconditioners = "45d422c2-293f-44ce-8315-2cb988662dec"

examples/eigenvalues.jl

@@ -0,0 +1,55 @@
+using NewtonKrylov
+using Arpack
+
+using SparseArrays, LinearAlgebra
+using CairoMakie
+
+# ## 1D Bratu equations
+# $y′′ + λ * exp(y) = 0$
+
+# F(u) = 0
+
+function bratu!(res, y, Δx, λ)
+    N = length(y)
+    for i in 1:N
+        y_l = i == 1 ? zero(eltype(y)) : y[i - 1]
+        y_r = i == N ? zero(eltype(y)) : y[i + 1]
+        y′′ = (y_r - 2y[i] + y_l) / Δx^2
+
+        res[i] = y′′ + λ * exp(y[i]) # = 0
+    end
+    return nothing
+end
+
+function bratu(y, dx, λ)
+    res = similar(y)
+    bratu!(res, y, dx, λ)
+    return res
+end
+
+# ## Choice of parameters
+const N = 500
+const λ = 3.51382
+const dx = 1 / (N + 1) # Grid-spacing
+
+# ### Domain and Inital condition
+x = LinRange(0.0 + dx, 1.0 - dx, N)
+u₀ = sin.(x .* π)
+
+JOp = NewtonKrylov.JacobianOperator(
+    (res, u) -> bratu!(res, u, dx, λ),
+    similar(u₀),
+    copy(u₀)
+)
+
+# Q: is this generally true
+LinearAlgebra.issymmetric(::NewtonKrylov.JacobianOperator) = true
+
+# issymmetric(collect(JOp))
+
+l, ϕ = eigs(JOp; nev = 10)
+l2, ϕ2 = eigs(collect(JOp); nev = 10)
+
+
+@time  eigs(collect(JOp); nev = 300)
+@time  eigs(JOp; nev = 300)

src/NewtonKrylov.jl

@@ -33,7 +33,7 @@ end
 Base.size(J::JacobianOperator) = (length(J.res), length(J.u))
 Base.eltype(J::JacobianOperator) = eltype(J.u)
 
-function mul!(out, J::JacobianOperator, v)
+function mul!(out::A, J::JacobianOperator{F,A}, v::A) where {F, A}
     # Enzyme.make_zero!(J.f_cache)
     f_cache = Enzyme.make_zero(J.f) # Stop gap until we can zero out mutable values
     autodiff(
@@ -45,6 +45,14 @@ function mul!(out, J::JacobianOperator, v)
     return nothing
 end
 
+function mul!(out, J::JacobianOperator, v)
+    _v = similar(J.u)
+    _v .= v
+    _out = similar(J.res)
+    mul!(_out, J, _v)
+    out .= _out
+end
+
 LinearAlgebra.adjoint(J::JacobianOperator) = Adjoint(J)
 LinearAlgebra.transpose(J::JacobianOperator) = Transpose(J)