From 89ff014cf9429f0e50839cf58629bf35dc4119bb Mon Sep 17 00:00:00 2001
From: Valentin Churavy <v.churavy@gmail.com>
Date: Tue, 10 Dec 2024 06:08:50 +0100
Subject: [PATCH 1/2] Extract eigenvalues

---
 examples/Project.toml   |  1 +
 examples/eigenvalues.jl | 55 +++++++++++++++++++++++++++++++++++++++++
 src/NewtonKrylov.jl     | 10 +++++++-
 3 files changed, 65 insertions(+), 1 deletion(-)
 create mode 100644 examples/eigenvalues.jl

diff --git a/examples/Project.toml b/examples/Project.toml
index 3671d5c..64aea47 100644
--- a/examples/Project.toml
+++ b/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"
diff --git a/examples/eigenvalues.jl b/examples/eigenvalues.jl
new file mode 100644
index 0000000..363e2cb
--- /dev/null
+++ b/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)
\ No newline at end of file
diff --git a/src/NewtonKrylov.jl b/src/NewtonKrylov.jl
index db993aa..023a173 100644
--- a/src/NewtonKrylov.jl
+++ b/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)
 

From 7499157f120b65bf2bb5fbacf85e0307c5c674ef Mon Sep 17 00:00:00 2001
From: Valentin Churavy <v.churavy@gmail.com>
Date: Tue, 17 Dec 2024 12:26:49 +0100
Subject: [PATCH 2/2] make symmetric configurable for JacobianOperator

---
 examples/eigenvalues.jl | 8 +++-----
 src/NewtonKrylov.jl     | 7 +++++--
 2 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/examples/eigenvalues.jl b/examples/eigenvalues.jl
index 363e2cb..65b46c0 100644
--- a/examples/eigenvalues.jl
+++ b/examples/eigenvalues.jl
@@ -39,13 +39,11 @@ u₀ = sin.(x .* π)
 JOp = NewtonKrylov.JacobianOperator(
     (res, u) -> bratu!(res, u, dx, λ),
     similar(u₀),
-    copy(u₀)
+    copy(u₀);
+    symmetric = true
 )
 
-# Q: is this generally true
-LinearAlgebra.issymmetric(::NewtonKrylov.JacobianOperator) = true
-
-# issymmetric(collect(JOp))
+@assert issymmetric(collect(JOp))
 
 l, ϕ = eigs(JOp; nev = 10)
 l2, ϕ2 = eigs(collect(JOp); nev = 10)
diff --git a/src/NewtonKrylov.jl b/src/NewtonKrylov.jl
index 023a173..68cfd5b 100644
--- a/src/NewtonKrylov.jl
+++ b/src/NewtonKrylov.jl
@@ -24,15 +24,18 @@ struct JacobianOperator{F, A}
     f_cache::F
     res::A
     u::A
-    function JacobianOperator(f::F, res, u) where {F}
+    symmetric::Union{Bool,Missing}
+    function JacobianOperator(f::F, res, u; symmetric=missing) where {F}
         f_cache = Enzyme.make_zero(f)
-        return new{F, typeof(u)}(f, f_cache, res, u)
+        return new{F, typeof(u)}(f, f_cache, res, u, symmetric, hermetian)
     end
 end
 
 Base.size(J::JacobianOperator) = (length(J.res), length(J.u))
 Base.eltype(J::JacobianOperator) = eltype(J.u)
 
+LinearAlgebra.issymmetric(J::JacobianOperator) = J.symmetric
+
 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