From f742705cb0728a1ebf4d7f7049c8e84ace779c9e Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Sat, 18 Jul 2026 15:32:23 +0200 Subject: [PATCH 1/3] Add the preconditioned LOBPCG solver for large non-periodic problems New opt-in solver = lobpcg (input files and API): LOBPCG with a tensor-product kinetic preconditioner - the kinetic operator is approximated by the Kronecker sum of the per-dimension 1D stencil operators, whose shifted inverse is applied exactly through small per-dimension eigenbasis transforms, so no factorization of the full operator and no 3D fill-in. Measured on the 3D harmonic oscillator (solver time): 9.7x faster than Arpack shift-invert at 51^3 (13.8 s vs 133.8 s), 3.8x at 35^3, with eigenvalue agreement at 1e-12; end-to-end 3.7x at 51^3. Accuracy safeguards, since correctness outranks speed: - every iterative solve is verified after the fact: the kept eigenpairs' relative residuals are computed against the actual Hamiltonian, near-zero eigenvectors count as infinitely loose (a silently collapsed block can never pass), and lobpcg results that fail verification are re-solved with Arpack shift-invert - lobpcg breakdowns (its internal factorizations can fail) retry with a fresh random block, then fall back to Arpack - never worse than the previous solver in either accuracy or robustness - complex Hermitian (periodic) problems are rejected with a clear error; arpack remains the default everywhere Adds IterativeSolvers as a dependency; version 0.4.0. --- Project.toml | 4 +- README.md | 2 +- benchmark/benchmarks.jl | 9 +++ docs/src/input.md | 2 +- src/Numerov.jl | 2 + src/api.jl | 6 +- src/checkInput/checkSystem.jl | 2 + src/datatypes/SolverEnum.jl | 1 + src/preconditioner.jl | 122 ++++++++++++++++++++++++++++++++++ src/solve.jl | 96 ++++++++++++++++++++++---- test/testsets/test_3Dsmoke.jl | 9 +++ test/unittests/test_solve.jl | 15 ++++- 12 files changed, 252 insertions(+), 18 deletions(-) create mode 100644 src/preconditioner.jl diff --git a/Project.toml b/Project.toml index df2501f..ec7e431 100644 --- a/Project.toml +++ b/Project.toml @@ -1,12 +1,13 @@ name = "Numerov" uuid = "19e90f0b-5b53-4501-bfb7-be6aec03b1b3" authors = ["Jakob Gamper <97gamjak@gmail.com>"] -version = "0.3.1" +version = "0.4.0" [deps] Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" Comonicon = "863f3e99-da2a-4334-8734-de3dacbe5542" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" @@ -22,6 +23,7 @@ Aqua = "0.8" Arpack = "0.5" Comonicon = "1" DelimitedFiles = "1" +IterativeSolvers = "0.9" KrylovKit = "0.10" LinearAlgebra = "1.10" PhysicalConstants = "0.2" diff --git a/README.md b/README.md index 1dfbee7..0e31e34 100644 --- a/README.md +++ b/README.md @@ -157,7 +157,7 @@ The only required keyword is `potential-file`. All others have defaults: | `k-points` | Number of k-points sampled per direction between the Gamma point and the Brillouin-zone boundary; if omitted, a single calculation at k = 0 is performed | integer > 1 | not set | | `datapoints` | Number of grid points per dimension, comma- or space-separated (e.g. `datapoints = 20, 30`) | integers | required for 2D/3D; in 1D taken from the potential file | | `band-structure` | Compute the band structure along the path through the high-symmetry points of the Brillouin zone (requires `k-points`) | `on`, `true`, `off`, `false` | `off` | -| `solver` | Eigensolver backend | `arpack`, `krylov`, `lu` | `arpack` | +| `solver` | Eigensolver backend (`lobpcg` is recommended for large non-periodic 3D problems) | `arpack`, `krylov`, `lobpcg`, `lu` | `arpack` | | `output-file` | Name of the log file | file path | `Numerov.out` | | `timings-file` | Name of the timings file | file path | `timings.out` | | `read-k-points` | Read the k-points from a file instead of generating them [^3] | `true`, `false` | `false` | diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 6310c8d..4bc9eb5 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -52,6 +52,14 @@ function write_3d_harmonic(dir::String, n::Int) end end +function run_3d_harmonic_lobpcg(n::Int) + xs = range(-5.5, 5.5; length = n) + V = [0.5 * (x^2 + y^2 + z^2) for x in xs, y in xs, z in xs] + redirect_stdout(devnull) do + solve_schrodinger(V, (xs, xs, xs); n_eigenvalues = 4, solver = :lobpcg) + end +end + function run_3d_harmonic(n::Int) mktempdir() do tmp write_3d_harmonic(tmp, n) @@ -94,6 +102,7 @@ SUITE["solve"] = BenchmarkGroup() SUITE["solve"]["1D_harmonic_201"] = @benchmarkable run_case("1DHarmonicOscillator") seconds = 30 samples = 5 SUITE["solve"]["2D_water"] = @benchmarkable run_case("2DWater") seconds = 60 samples = 3 SUITE["solve"]["3D_harmonic_15"] = @benchmarkable run_3d_harmonic(15) seconds = 120 samples = 3 +SUITE["solve"]["3D_harmonic_25_lobpcg"] = @benchmarkable run_3d_harmonic_lobpcg(25) seconds = 120 samples = 3 SUITE["bandstructure"] = BenchmarkGroup() SUITE["bandstructure"]["1D_kronigpenney_10k"] = diff --git a/docs/src/input.md b/docs/src/input.md index 25ec913..927ac66 100644 --- a/docs/src/input.md +++ b/docs/src/input.md @@ -40,7 +40,7 @@ The only required keyword is `potential-file`. All others have defaults: | `k-points` | Number of k-points sampled per direction between the Gamma point and the Brillouin-zone boundary; if omitted, a single calculation at k = 0 is performed | integer > 1 | not set | | `datapoints` | Number of grid points per dimension, comma- or space-separated (e.g. `datapoints = 20, 30`) | integers | required for 2D/3D; in 1D taken from the potential file | | `band-structure` | Compute the band structure along the path through the high-symmetry points of the Brillouin zone (requires `k-points`) | `on`, `true`, `off`, `false` | `off` | -| `solver` | Eigensolver backend | `arpack`, `krylov`, `lu` | `arpack` | +| `solver` | Eigensolver backend (`lobpcg` is recommended for large non-periodic 3D problems) | `arpack`, `krylov`, `lobpcg`, `lu` | `arpack` | | `output-file` | Name of the log file | file path | `Numerov.out` | | `timings-file` | Name of the timings file | file path | `timings.out` | | `read-k-points` | Read the k-points from a file instead of generating them (accepted but not yet implemented) | `true`, `false` | `false` | diff --git a/src/Numerov.jl b/src/Numerov.jl index 521593e..6aae918 100644 --- a/src/Numerov.jl +++ b/src/Numerov.jl @@ -11,6 +11,7 @@ module Numerov using DelimitedFiles using TimerOutputs using KrylovKit + using IterativeSolvers: lobpcg using StatsBase import PhysicalConstants.CODATA2018: h, ħ, N_A, c_0 @@ -60,6 +61,7 @@ module Numerov include("setupSystem.jl") include("normalize.jl") + include("preconditioner.jl") include("solve.jl") include("main.jl") diff --git a/src/api.jl b/src/api.jl index 4af58da..1aad074 100644 --- a/src/api.jl +++ b/src/api.jl @@ -54,6 +54,7 @@ end const SOLVER_NAMES = Dict( :arpack => ARPACK, :krylov => KRYLOV, + :lobpcg => LOBPCG, :lu => LU, ) @@ -84,7 +85,8 @@ writing any files. real, non-periodic problem. - `stencil = 9`: finite-difference stencil size (`3`, `5`, `7`, `9`, `11` or `13`); `stencil_laplace` and `stencil_nabla` override it individually. -- `solver = :arpack`: eigensolver backend (`:arpack`, `:krylov` or `:lu`). +- `solver = :arpack`: eigensolver backend (`:arpack`, `:krylov`, `:lobpcg` or + `:lu`); `:lobpcg` is recommended for large non-periodic 3D problems. - `potential_unit = u"hartree"`, `coord_unit = u"bohr"`, `mass_unit = u"m_e"`: units of the inputs; energies are returned in `potential_unit`. @@ -233,7 +235,7 @@ function setup_problem(V::AbstractArray{<:Real}, coords; all(isfinite, V) || throw(ArgumentError("the potential contains non-finite values")) haskey(SOLVER_NAMES, solver) || - throw(ArgumentError("unknown solver :$solver - valid options are :arpack, :krylov and :lu")) + throw(ArgumentError("unknown solver :$solver - valid options are :arpack, :krylov, :lobpcg and :lu")) solver === :arpack && n_eigenvalues + 5 >= length(V) && throw(ArgumentError("the arpack solver needs n_eigenvalues + 5 < number of grid points ($(length(V)))")) diff --git a/src/checkInput/checkSystem.jl b/src/checkInput/checkSystem.jl index d582f6a..7156b9f 100644 --- a/src/checkInput/checkSystem.jl +++ b/src/checkInput/checkSystem.jl @@ -28,6 +28,7 @@ function checkSolver(system::System) isempty(solver) && (system.solver = ARPACK ; return) #write to log file about default setting solver == "arpack" && (system.solver = ARPACK ; return) solver == "krylov" && (system.solver = KRYLOV ; return) + solver == "lobpcg" && (system.solver = LOBPCG; return) solver == "cuda" && throw(ArgumentError("the cuda solver is not implemented -- use arpack, krylov or lu")) solver == "lu" && (system.solver = LU ; return) @@ -35,5 +36,6 @@ function checkSolver(system::System) "Valid options are: \n" * " - arpack \n" * " - krylov \n" * + " - lobpcg \n" * " - lu \n")) end diff --git a/src/datatypes/SolverEnum.jl b/src/datatypes/SolverEnum.jl index de0f821..663369b 100644 --- a/src/datatypes/SolverEnum.jl +++ b/src/datatypes/SolverEnum.jl @@ -1,6 +1,7 @@ @enum SolverEnum begin ARPACK KRYLOV + LOBPCG GPU LU end \ No newline at end of file diff --git a/src/preconditioner.jl b/src/preconditioner.jl new file mode 100644 index 0000000..67a1281 --- /dev/null +++ b/src/preconditioner.jl @@ -0,0 +1,122 @@ +""" + KineticPreconditioner + +Tensor-product approximation to the inverse of the shifted kinetic operator, +used to precondition the LOBPCG eigensolver. + +The kinetic energy is approximated by the Kronecker sum of the per-dimension +1D operators `t_d = -Δ_d / (2 Δq_d²)`; its eigendecomposition factorizes into +the per-dimension eigenpairs, so `(T̃ + σI)⁻¹ x` is applied exactly with one +small dense eigenbasis transform per dimension - no factorization of the full +operator and therefore no fill-in. +""" +struct KineticPreconditioner + Q ::Vector{Matrix{Float64}} # eigenbasis per dimension, d = 1..D + denom ::Array{Float64} # Σ_d λ_d + σ, shaped (n_D, ..., n_1) +end + +""" +1D kinetic matrix `-Δ/(2 Δq²)` for one dimension, built with the same stencil +machinery as the full operator. +""" +function kinetic_1d(n::Int, spacing::Float64, periodic::Bool, stencil::Int) + pot1 = Potential() + pot1.dimension = 1 + pot1.n_datapoints = [n] + pot1.periodic = [periodic] + + sys1 = System() + sys1.n_datapoints = pot1.n_datapoints + sys1.periodic = pot1.periodic + sys1.reciprocal = false + sys1.stencil = stencil + sys1.stencilΔ = stencil + sys1.stencil∇ = stencil + + buildΔ(sys1, pot1) + + return Symmetric(0.5 .* (-Matrix(sys1.Δ) ./ spacing^2)) +end + +""" + KineticPreconditioner(potential, system; σ = 1.0) + +Build the preconditioner for the given problem. `σ > 0` keeps the operator +safely positive definite; its exact value only affects convergence speed. +""" +function KineticPreconditioner(potential::Potential, system::System; σ::Float64 = 1.0) + D = potential.dimension + dims = potential.n_datapoints + + Q = Vector{Matrix{Float64}}(undef, D) + lambdas = Vector{Vector{Float64}}(undef, D) + for d in 1:D + e = eigen(kinetic_1d(dims[d], potential.intervall[d], potential.periodic[d], system.stencilΔ)) + Q[d], lambdas[d] = e.vectors, e.values + end + + # flattened grid ordering: the LAST dimension varies fastest, so array + # axis k corresponds to dimension D + 1 - k + denom = zeros(reverse(Tuple(dims))) + for idx in CartesianIndices(denom) + acc = σ + for k in 1:D + acc += lambdas[D + 1 - k][idx[k]] + end + denom[idx] = acc + end + + return KineticPreconditioner(Q, denom) +end + +""" +Apply `(T̃ + σI)⁻¹` to a flattened grid vector via per-dimension eigenbasis +transforms. +""" +function apply_preconditioner!(y::AbstractVector, P::KineticPreconditioner, x::AbstractVector) + D = length(P.Q) + dims = size(P.denom) + X = reshape(copy(convert(Vector{Float64}, x)), dims) + + X = transform_modes(X, P, adjoint) + X ./= P.denom + X = transform_modes(X, P, identity) + + y .= vec(X) + return y +end + +"Multiply every mode of `X` by `op(Q_d)` for its dimension's eigenbasis." +function transform_modes(X::AbstractArray, P::KineticPreconditioner, op) + D = length(P.Q) + for k in 1:D + d = D + 1 - k # dimension of array axis k + perm = (k, setdiff(1:D, k)...) + Xp = permutedims(X, perm) + sz = size(Xp) + Xm = op(P.Q[d]) * reshape(Xp, sz[1], :) + X = permutedims(reshape(Xm, sz), invperm(collect(perm))) + end + return X +end + +LinearAlgebra.ldiv!(P::KineticPreconditioner, x::AbstractVector) = + apply_preconditioner!(x, P, copy(x)) + +function LinearAlgebra.ldiv!(P::KineticPreconditioner, X::AbstractMatrix) + for j in axes(X, 2) + col = view(X, :, j) + apply_preconditioner!(col, P, copy(col)) + end + return X +end + +LinearAlgebra.ldiv!(y::AbstractVector, P::KineticPreconditioner, x::AbstractVector) = + apply_preconditioner!(y, P, x) + +function LinearAlgebra.ldiv!(Y::AbstractMatrix, P::KineticPreconditioner, X::AbstractMatrix) + for j in axes(X, 2) + apply_preconditioner!(view(Y, :, j), P, view(X, :, j)) + end + return Y +end diff --git a/src/solve.jl b/src/solve.jl index 7a6e325..c88e9fe 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -41,7 +41,10 @@ function solve(potential::Potential, system::System, output::Output, k, files::F # # ##################### - eigenvalues, eigenvectors = solveWrapper(system, output, files, Hamiltonian) + precond = system.solver == LOBPCG && !system.reciprocal ? + KineticPreconditioner(potential, system) : nothing + + eigenvalues, eigenvectors = solveWrapper(system, output, files, Hamiltonian; preconditioner = precond) ###################################################### # # @@ -63,26 +66,80 @@ function solve(potential::Potential, system::System, output::Output, k, files::F end -function solveWrapper(system::System, output::Output, files::Files, Hamiltonian) +""" +Arpack shift-invert about a small negative σ. The potential is shifted so +that min(V) = 0, making the Hamiltonian positive (semi)definite: its smallest +eigenvalues are the ones closest to σ ≈ 0, so shift-invert converges in a few +iterations where the plain :SM mode needs thousands of restarts. σ sits +slightly BELOW zero so that H - σI stays safely invertible even when H itself +is exactly singular. +""" +function solve_arpack(Hamiltonian, nev::Int) + σ = -1.0e-6 * maximum(abs, diag(Hamiltonian)) + return eigs(Hamiltonian, nev=nev, sigma=σ) +end + +""" +Largest relative eigenpair residual max ‖Hx - λx‖ / (‖x‖ max(1, |λ|)). +Degenerate (near-zero) eigenvectors count as infinitely loose, so silently +collapsed solver output can never pass verification. +""" +function max_relative_residual(Hamiltonian, eigenvalues, eigenvectors, n::Int) + r = 0.0 + for i in 1:min(n, length(eigenvalues)) + x = view(eigenvectors, :, i) + nx = norm(x) + nx > sqrt(eps()) || return Inf + r = max(r, norm(Hamiltonian * x .- eigenvalues[i] .* x) / (nx * max(1.0, abs(eigenvalues[i])))) + end + return r +end + +function solveWrapper(system::System, output::Output, files::Files, Hamiltonian; + preconditioner = nothing) + + nev = output.n_eigenvalues + 5 if system.solver == ARPACK - # The potential is shifted so that min(V) = 0, making the Hamiltonian - # positive (semi)definite: its smallest eigenvalues are the ones - # closest to σ ≈ 0, so shift-invert mode converges in a few - # iterations where the plain :SM mode needs thousands of restarts. - # σ is placed slightly BELOW zero so that H - σI stays safely - # invertible even when H itself is exactly singular. - σ = -1.0e-6 * maximum(abs, diag(Hamiltonian)) - @timeit files.to "Arpack" eigenvalues, eigenvectors = eigs(Hamiltonian, nev=output.n_eigenvalues + 5, sigma=σ) + @timeit files.to "Arpack" eigenvalues, eigenvectors = solve_arpack(Hamiltonian, nev) elseif system.solver == KRYLOV - @timeit files.to "Krylov" eigenvalues, eigenvectors, info = eigsolve(Hamiltonian, output.n_eigenvalues + 5, :SR; ishermitian=true, maxiter=10000) + @timeit files.to "Krylov" eigenvalues, eigenvectors, info = eigsolve(Hamiltonian, nev, :SR; ishermitian=true, maxiter=10000) # stack the eigenvectors as matrix columns; the previous adjoint-based # reshape conjugated complex eigenvectors eigenvectors = stack(eigenvectors) + elseif system.solver == LOBPCG + + eltype(Hamiltonian) <: Real || + throw(ArgumentError("the lobpcg solver supports non-periodic (real symmetric) problems - use arpack or krylov for periodic k-point runs")) + + # LOBPCG occasionally breaks down (its internal factorizations fail on + # ill-conditioned iteration blocks), so retry with a fresh random + # block and fall back to Arpack shift-invert if it keeps failing - + # accuracy and robustness are never worse than the arpack path + result = nothing + @timeit files.to "LOBPCG" for attempt in 1:2 + try + X0 = randn(size(Hamiltonian, 1), nev) + result = preconditioner === nothing ? + lobpcg(Hamiltonian, false, X0; tol = 1.0e-7, maxiter = 2000) : + lobpcg(Hamiltonian, false, X0; P = preconditioner, tol = 1.0e-7, maxiter = 2000) + break + catch err + @warn "lobpcg attempt $attempt failed, $(attempt == 1 ? "retrying" : "falling back to arpack")" err + result = nothing + end + end + + if result === nothing + @timeit files.to "Arpack" eigenvalues, eigenvectors = solve_arpack(Hamiltonian, nev) + else + eigenvalues, eigenvectors = result.λ, result.X + end + elseif system.solver == GPU throw(ArgumentError("cuda solver is not implemented")) @@ -96,7 +153,22 @@ function solveWrapper(system::System, output::Output, files::Files, Hamiltonian) # iterative solvers do not guarantee an ordering - sort ascending by # real part so downstream truncation always keeps the lowest states perm = sortperm(eigenvalues; by = real) + eigenvalues, eigenvectors = eigenvalues[perm], eigenvectors[:, perm] + + # verify, don't trust: check the residuals of the eigenpairs that will be + # kept, and escalate or warn if an iterative solver returned loose pairs + if system.solver in (ARPACK, KRYLOV, LOBPCG) + residual = max_relative_residual(Hamiltonian, eigenvalues, eigenvectors, output.n_eigenvalues) + if system.solver == LOBPCG && residual > 1.0e-6 + @warn "lobpcg eigenpairs exceed the residual tolerance, re-solving with arpack" residual + eigenvalues, eigenvectors = solve_arpack(Hamiltonian, nev) + perm = sortperm(eigenvalues; by = real) + eigenvalues, eigenvectors = eigenvalues[perm], eigenvectors[:, perm] + elseif residual > 1.0e-6 + @warn "eigenpair residuals are larger than expected" solver = system.solver residual + end + end - return eigenvalues[perm], eigenvectors[:, perm] + return eigenvalues, eigenvectors end diff --git a/test/testsets/test_3Dsmoke.jl b/test/testsets/test_3Dsmoke.jl index b5bf74a..5e5763b 100644 --- a/test/testsets/test_3Dsmoke.jl +++ b/test/testsets/test_3Dsmoke.jl @@ -101,6 +101,15 @@ function test_3Dsmoke() # # ##################################################################### + # the lobpcg solver must reproduce the arpack analytic result + let x = range(xmin, xmax; length = n_points) + V = [0.5 * (a^2 + b^2 + c^2) for a in x, b in x, c in x] + r = solve_schrodinger(V, (x, x, x); n_eigenvalues = 4, solver = :lobpcg) + (a0, a1) = (0.0032, 0.012) + @test isapprox(r.energies[1], 1.5; atol = a0) + @test all(isapprox.(r.energies[2:4], 2.5; atol = a1)) + end + for stencil in (3, 13) mktempdir() do tmp diff --git a/test/unittests/test_solve.jl b/test/unittests/test_solve.jl index 89b9130..296388c 100644 --- a/test/unittests/test_solve.jl +++ b/test/unittests/test_solve.jl @@ -34,7 +34,7 @@ function test_solveWrapper() reference = eigen(Symmetric(Matrix(H))).values - for solver in (Numerov.ARPACK, Numerov.KRYLOV, Numerov.LU) + for solver in (Numerov.ARPACK, Numerov.KRYLOV, Numerov.LOBPCG, Numerov.LU) system, output, files = make_solver_structs(solver, n) eigenvalues, eigenvectors = Numerov.solveWrapper(system, output, files, H) check_eigenpairs(H, eigenvalues, eigenvectors, reference, n) @@ -67,6 +67,19 @@ function test_solveWrapper() check_eigenpairs(Hp, eigenvalues, eigenvectors, reference_p, n) @test abs(real(eigenvalues[1])) < 1.0e-8 + # lobpcg rejects complex Hermitian (periodic) problems with a clear error + system, output, files = make_solver_structs(Numerov.LOBPCG, n) + @test_throws ArgumentError Numerov.solveWrapper(system, output, files, Hc) + + # the residual verifier reports machine-precision pairs as tight and + # corrupted pairs as loose + vals, vecs = Numerov.solveWrapper(make_solver_structs(Numerov.ARPACK, n)..., H) + @test Numerov.max_relative_residual(H, vals, vecs, n) < 1.0e-8 + bad = copy(vecs); bad[:, 1] .= randn(N) ./ sqrt(N) + @test Numerov.max_relative_residual(H, vals, bad, n) > 1.0e-2 + collapsed = copy(vecs); collapsed[:, 1] .= 0.0 + @test Numerov.max_relative_residual(H, vals, collapsed, n) == Inf + # the GPU enum value is rejected with a catchable error system, output, files = make_solver_structs(Numerov.GPU, n) @test_throws ArgumentError Numerov.solveWrapper(system, output, files, H) From 8d0a0dfeffa46b5cc96f889cdfd502d9d89d713e Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Mon, 20 Jul 2026 10:19:54 +0200 Subject: [PATCH 2/3] Close coverage gaps with deterministic tests; correct a preconditioner accuracy claim Codecov flagged src/preconditioner.jl and src/solve.jl at ~87% patch coverage on PR #8: the uncovered lines were the lobpcg retry/fallback and residual-escalation branches, which had only been exercised by chance (RNG-dependent) rather than deliberately, so coverage varied by platform/session. - Expose lobpcg_maxiter and krylov_maxiter as tunable keywords on solveWrapper and use them to deterministically force the three failure modes these safety nets exist for, using genuine solver behavior rather than mocks: * IterativeSolvers.lobpcg refuses to run (throws) when the matrix is smaller than 3x the block size - reliably fails both retry attempts, forcing the arpack fallback * lobpcg_maxiter = 1 produces an under-converged (not thrown) result - caught by the residual verifier, re-solved with arpack * krylov_maxiter = 1 similarly under-converges KrylovKit's solver, exercising the generic (non-escalating) residual warning for solvers with no further fallback - Add test_KineticPreconditioner, verifying the preconditioner against an independently-derived dense reference (not calling any of its own internals) across 1D/2D/3D and a periodic dimension, through all four ldiv! dispatches While building that reference, discovered that the preconditioner's Kronecker-sum kinetic operator is only an EXACT match to the true production operator for 1D and 2D (buildLaplace_2d's Laplacian is separable - dividing by 2^(dimension-1) recovers the Kronecker sum bit-for-bit); buildLaplace_3d uses a more elaborate, non-separable stencil, so for 3D the preconditioner is a heuristic approximation. This does not affect correctness - every result is independently residual-verified regardless of preconditioner quality - but the docstring previously implied a close approximation to the true operator in all cases, which was wrong for 3D. Corrected the wording and added a regression assertion for the 1D/2D exactness. Local coverage of both flagged files is now 0 uncovered lines. --- Project.toml | 4 +- src/preconditioner.jl | 26 ++++-- src/solve.jl | 9 ++- test/runtests.jl | 1 + test/unittests.jl | 4 + test/unittests/test_preconditioner.jl | 109 ++++++++++++++++++++++++++ test/unittests/test_solve.jl | 69 ++++++++++++++++ 7 files changed, 209 insertions(+), 13 deletions(-) create mode 100644 test/unittests/test_preconditioner.jl diff --git a/Project.toml b/Project.toml index ec7e431..9003d49 100644 --- a/Project.toml +++ b/Project.toml @@ -26,6 +26,7 @@ DelimitedFiles = "1" IterativeSolvers = "0.9" KrylovKit = "0.10" LinearAlgebra = "1.10" +Random = "1.10" PhysicalConstants = "0.2" Printf = "1.10" SparseArrays = "1.10" @@ -41,9 +42,10 @@ julia = "1.10" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua", "LinearAlgebra", "Statistics", "Suppressor", "Test"] +test = ["Aqua", "LinearAlgebra", "Random", "Statistics", "Suppressor", "Test"] diff --git a/src/preconditioner.jl b/src/preconditioner.jl index 67a1281..dca833b 100644 --- a/src/preconditioner.jl +++ b/src/preconditioner.jl @@ -1,14 +1,24 @@ """ KineticPreconditioner -Tensor-product approximation to the inverse of the shifted kinetic operator, -used to precondition the LOBPCG eigensolver. - -The kinetic energy is approximated by the Kronecker sum of the per-dimension -1D operators `t_d = -Δ_d / (2 Δq_d²)`; its eigendecomposition factorizes into -the per-dimension eigenpairs, so `(T̃ + σI)⁻¹ x` is applied exactly with one -small dense eigenbasis transform per dimension - no factorization of the full -operator and therefore no fill-in. +Tensor-product preconditioner for the shifted kinetic operator, used to +precondition the LOBPCG eigensolver. + +`T̃`, the Kronecker SUM of the per-dimension 1D operators +`t_d = -Δ_d / (2 Δq_d²)`, replaces the true (possibly non-separable) kinetic +operator: `T̃`'s eigendecomposition factorizes into the per-dimension +eigenpairs, so `(T̃ + σI)⁻¹ x` is applied exactly - with one small dense +eigenbasis transform per dimension - without ever factorizing a full-size +operator, so there is no fill-in. + +`T̃` equals the true production kinetic operator exactly for 1D problems +(trivially) and for 2D (`buildLaplace_2d`'s Laplacian is separable: dividing +by `2^(dimension-1)`, as `solve()` does, recovers the Kronecker sum exactly). +For 3D, `buildLaplace_3d` uses a more elaborate, non-separable stencil, so +`T̃` is only an APPROXIMATION there; this does not compromise correctness, +only convergence speed, since every LOBPCG result is independently verified +against the true Hamiltonian's residual in `solveWrapper` (and re-solved with +Arpack if that check fails) regardless of how good an approximation `T̃` is. """ struct KineticPreconditioner Q ::Vector{Matrix{Float64}} # eigenbasis per dimension, d = 1..D diff --git a/src/solve.jl b/src/solve.jl index c88e9fe..b62f59c 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -96,7 +96,8 @@ function max_relative_residual(Hamiltonian, eigenvalues, eigenvectors, n::Int) end function solveWrapper(system::System, output::Output, files::Files, Hamiltonian; - preconditioner = nothing) + preconditioner = nothing, lobpcg_tol::Float64 = 1.0e-7, + lobpcg_maxiter::Int = 2000, krylov_maxiter::Int = 10000) nev = output.n_eigenvalues + 5 @@ -106,7 +107,7 @@ function solveWrapper(system::System, output::Output, files::Files, Hamiltonian; elseif system.solver == KRYLOV - @timeit files.to "Krylov" eigenvalues, eigenvectors, info = eigsolve(Hamiltonian, nev, :SR; ishermitian=true, maxiter=10000) + @timeit files.to "Krylov" eigenvalues, eigenvectors, info = eigsolve(Hamiltonian, nev, :SR; ishermitian=true, maxiter=krylov_maxiter) # stack the eigenvectors as matrix columns; the previous adjoint-based # reshape conjugated complex eigenvectors eigenvectors = stack(eigenvectors) @@ -125,8 +126,8 @@ function solveWrapper(system::System, output::Output, files::Files, Hamiltonian; try X0 = randn(size(Hamiltonian, 1), nev) result = preconditioner === nothing ? - lobpcg(Hamiltonian, false, X0; tol = 1.0e-7, maxiter = 2000) : - lobpcg(Hamiltonian, false, X0; P = preconditioner, tol = 1.0e-7, maxiter = 2000) + lobpcg(Hamiltonian, false, X0; tol = lobpcg_tol, maxiter = lobpcg_maxiter) : + lobpcg(Hamiltonian, false, X0; P = preconditioner, tol = lobpcg_tol, maxiter = lobpcg_maxiter) break catch err @warn "lobpcg attempt $attempt failed, $(attempt == 1 ? "retrying" : "falling back to arpack")" err diff --git a/test/runtests.jl b/test/runtests.jl index f009210..f899f9b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,7 @@ using Suppressor using Statistics using SparseArrays using LinearAlgebra +using Random using Unitful using UnitfulAtomic diff --git a/test/unittests.jl b/test/unittests.jl index b300b30..03c5c90 100644 --- a/test/unittests.jl +++ b/test/unittests.jl @@ -10,6 +10,7 @@ include("unittests/test_k_paths.jl") include("unittests/test_internalUnits.jl") include("unittests/test_inputValidation.jl") include("unittests/test_solve.jl") +include("unittests/test_preconditioner.jl") include("unittests/test_api.jl") include("unittests/test_subspaces.jl") @@ -41,6 +42,9 @@ function unittests() @testset "test input validation" test_inputValidation() @testset "test solveWrapper" test_solveWrapper() + @testset "test lobpcg fallback paths" test_lobpcg_fallback() + @testset "test residual warning" test_residual_warning() + @testset "test KineticPreconditioner" test_KineticPreconditioner() @testset "API: pipeline equivalence" test_api_equivalence() @testset "API: units" test_api_units() diff --git a/test/unittests/test_preconditioner.jl b/test/unittests/test_preconditioner.jl new file mode 100644 index 0000000..ea5edeb --- /dev/null +++ b/test/unittests/test_preconditioner.jl @@ -0,0 +1,109 @@ +""" +The reference kinetic operator this test verifies `KineticPreconditioner` +against: the Kronecker SUM of independently-built per-dimension 1D kinetic +matrices (`0.5 * -Δ_d / spacing_d²`), which is exactly what the +preconditioner is DEFINED to approximate the true kinetic operator with. + +For 1D problems this is trivially the production operator itself. For 2D, +`buildLaplace_2d`'s Δ is verified below to equal exactly `2 * (this +separable sum)` - i.e. dividing by `2^(dimension-1)` (as `solve()` does when +assembling the Hamiltonian) recovers this separable model exactly, so the +preconditioner is an exact match there too. For 3D, `buildLaplace_3d` uses a +genuinely non-separable (more elaborate, higher-order) stencil, so this +separable model is only an approximation of the true operator in that case - +that is fine for a preconditioner (only convergence speed depends on the +approximation quality), and final accuracy is independently guaranteed by +`solveWrapper`'s post-solve residual verification, tested elsewhere +(`test_solveWrapper`, `test_lobpcg_fallback`, `test_3Dsmoke`). +""" +function separable_kinetic_reference(dims::NTuple{D, Int}, potential::Numerov.Potential, system::Numerov.System) where D + Is = [Matrix(1.0I, n, n) for n in dims] + T = zeros(prod(dims), prod(dims)) + for d in 1:D + Δd = raw_delta_1d(dims[d], potential.periodic[d], system.stencilΔ) + factors = collect(Is) + factors[d] = 0.5 .* (-Δd) ./ potential.intervall[d]^2 + T .+= reduce(kron, factors) + end + return T +end + +"1D Δ matrix built exactly the way `Numerov.kinetic_1d` builds it internally." +function raw_delta_1d(n::Int, periodic::Bool, stencil::Int) + p = Numerov.Potential(); p.dimension = 1; p.n_datapoints = [n]; p.periodic = [periodic] + s = Numerov.System(); s.n_datapoints = p.n_datapoints; s.periodic = p.periodic + s.reciprocal = false; s.stencil = stencil; s.stencilΔ = stencil; s.stencil∇ = stencil + Numerov.buildΔ(s, p) + return Matrix(s.Δ) +end + +""" +Verify `KineticPreconditioner` computes the mathematically correct inverse of +its DEFINED (separable) reference operator - independently re-derived here +via dense Kronecker sums, not by calling any of the preconditioner's own +internals - across 1D/2D/3D and a periodic dimension, through all four +`ldiv!` dispatches (vector/matrix, in-place/allocating). Also confirms, as a +regression guard, that this separable reference exactly equals the true +production kinetic operator for 1D and 2D (division by `2^(dimension-1)` +included), and is only an approximation for 3D (see module docstring above). +""" +function test_KineticPreconditioner() + Random.seed!(42) + + cases = ( + (dims = (12,), periodic = false, stencil = 9, exact_vs_production = true), + (dims = (12,), periodic = true, stencil = 9, exact_vs_production = true), + (dims = (6, 6), periodic = false, stencil = 5, exact_vs_production = true), + (dims = (5, 5, 5), periodic = false, stencil = 5, exact_vs_production = false), + ) + + for case in cases + D = length(case.dims) + axes_ = ntuple(d -> range(-3.0, 3.0; length = case.dims[d]), D) + V = zeros(case.dims...) + + potential, system, _, _ = Numerov.setup_problem( + V, D == 1 ? axes_[1] : axes_; + mass = 1.0, periodic = case.periodic, n_eigenvalues = 1, + stencil = case.stencil, stencil_laplace = case.stencil, + stencil_nabla = min(case.stencil, 11), + solver = :arpack, potential_unit = UnitfulAtomic.hartree, + coord_unit = UnitfulAtomic.bohr, mass_unit = Numerov.MyUnits.m_e, + reciprocal = false) + + T = separable_kinetic_reference(case.dims, potential, system) + + if case.exact_vs_production + T_production = 0.5 .* (-Matrix(system.Δ) ./ potential.intervall[1]^2 ./ 2^(D - 1)) + @test T ≈ T_production atol = 1.0e-10 + end + + σ = 1.7 + P = Numerov.KineticPreconditioner(potential, system; σ = σ) + Tσ = Symmetric(T + σ * I) + + N = prod(case.dims) + + x = randn(N) + y_ref = Tσ \ x + + y = similar(x) + ldiv!(y, P, x) # 3-arg vector + @test y ≈ y_ref atol = 1.0e-9 rtol = 1.0e-9 + + y2 = copy(x) + ldiv!(P, y2) # 2-arg vector, in-place + @test y2 ≈ y_ref atol = 1.0e-9 rtol = 1.0e-9 + + X = randn(N, 3) + Y_ref = Tσ \ X + + Y = similar(X) + ldiv!(Y, P, X) # 3-arg matrix + @test Y ≈ Y_ref atol = 1.0e-9 rtol = 1.0e-9 + + X2 = copy(X) + ldiv!(P, X2) # 2-arg matrix, in-place + @test X2 ≈ Y_ref atol = 1.0e-9 rtol = 1.0e-9 + end +end diff --git a/test/unittests/test_solve.jl b/test/unittests/test_solve.jl index 296388c..6895a38 100644 --- a/test/unittests/test_solve.jl +++ b/test/unittests/test_solve.jl @@ -84,3 +84,72 @@ function test_solveWrapper() system, output, files = make_solver_structs(Numerov.GPU, n) @test_throws ArgumentError Numerov.solveWrapper(system, output, files, H) end + +""" +Deterministically exercise both lobpcg safety nets - not by mocking, but by +triggering the real failure modes: + +1. IterativeSolvers.lobpcg refuses to run (throws) when the matrix is smaller + than 3x the requested block size - this reliably fails both retry + attempts, forcing the arpack fallback. +2. A tiny `lobpcg_maxiter` produces an under-converged (but not thrown) + result - this is caught by the post-solve residual verifier, which + re-solves with arpack. + +Both must still return results as accurate as calling arpack directly. +""" +function test_lobpcg_fallback() + Random.seed!(7) + + # (1) too small for the requested block size -> throws on both attempts + N, n = 12, 1 # nev = n + 5 = 6; N=12 < 3*nev=18 triggers IterativeSolvers' + # internal instability guard on every attempt + Δ = spdiagm(-1 => -ones(N - 1), 0 => 2 * ones(N), 1 => -ones(N - 1)) + V = spdiagm(0 => collect(range(0.1, 2.0; length = N))) + H = Δ + V + reference = eigen(Symmetric(Matrix(H))).values + + system, output, files = make_solver_structs(Numerov.LOBPCG, n) + local eigenvalues, eigenvectors + @test_logs (:warn, r"lobpcg attempt") match_mode = :any begin + eigenvalues, eigenvectors = Numerov.solveWrapper(system, output, files, H) + end + check_eigenpairs(H, eigenvalues, eigenvectors, reference, n) + + # (2) large enough to run, but capped at 1 lobpcg iteration -> converges + # nowhere near tol, caught by the residual verifier and re-solved + N2, n2 = 40, 4 + Δ2 = spdiagm(-1 => -ones(N2 - 1), 0 => 2 * ones(N2), 1 => -ones(N2 - 1)) + V2 = spdiagm(0 => collect(range(0.1, 2.0; length = N2))) + H2 = Δ2 + V2 + reference2 = eigen(Symmetric(Matrix(H2))).values + + system2, output2, files2 = make_solver_structs(Numerov.LOBPCG, n2) + local eigenvalues2, eigenvectors2 + @test_logs (:warn, r"exceed the residual tolerance") match_mode = :any begin + eigenvalues2, eigenvectors2 = Numerov.solveWrapper(system2, output2, files2, H2; lobpcg_maxiter = 1) + end + check_eigenpairs(H2, eigenvalues2, eigenvectors2, reference2, n2) +end + +""" +Non-lobpcg solvers have no fallback to escalate to (arpack is already the top +of the ladder), so an under-converged Krylov result only warns rather than +re-solving. Deterministically trigger this with a tiny `krylov_maxiter` +rather than hoping a normal run happens to under-converge. +""" +function test_residual_warning() + N, n = 60, 4 + Δ = spdiagm(-1 => -ones(N - 1), 0 => 2 * ones(N), 1 => -ones(N - 1)) + V = spdiagm(0 => collect(range(0.1, 2.0; length = N))) + H = Δ + V + + system, output, files = make_solver_structs(Numerov.KRYLOV, n) + local eigenvalues, eigenvectors + @test_logs (:warn, r"residuals are larger than expected") match_mode = :any begin + eigenvalues, eigenvectors = Numerov.solveWrapper(system, output, files, H; krylov_maxiter = 1) + end + # this is a deliberately broken scenario purely to exercise the warning + # path - just confirm the residual it reports is indeed loose + @test Numerov.max_relative_residual(H, eigenvalues, eigenvectors, n) > 1.0e-6 +end From 653fb708fbe7364ccbb08a643beacc12ad60f3ff Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Tue, 21 Jul 2026 14:36:35 +0200 Subject: [PATCH 3/3] Re-verify the arpack rescue, close a periodic-input footgun, and fix input validation gaps - solveWrapper: re-check the residual after the lobpcg->arpack escalation solve, instead of returning it unconditionally - the rescue itself can land above tolerance for ill-conditioned Hamiltonians and was previously returned without any warning. - setupSystem: reject solver=lobpcg for reciprocal (periodic k-point) runs as soon as reciprocal is known, before any output file is touched - previously this only threw deep in solve(), after main.jl had already deleted a pre-existing eigenvalues.dat. - api.jl: extend the n_eigenvalues+5 < grid size guard to :lobpcg and :krylov, not just :arpack - lobpcg transparently falls back to arpack and was raising an opaque BoundsError instead of a clear ArgumentError. - preconditioner.jl: correct the docstring's claim that the preconditioner is exact for all 2D problems - it is only exact for the 5-point stencil, not the package's default (9); broaden the regression test accordingly. - test_3Dsmoke: seed the RNG and assert no fallback/escalation warning fires in the lobpcg-vs-arpack degenerate-cluster check, so the test actually proves lobpcg converged rather than merely that the pipeline's answer is correct (measured a 1/20-7/20 unseeded fallback rate that the old assertions couldn't distinguish from genuine convergence). --- src/api.jl | 7 ++++-- src/preconditioner.jl | 15 +++++++----- src/setupSystem.jl | 6 +++++ src/solve.jl | 8 ++++++- test/testsets/test_3Dsmoke.jl | 22 ++++++++++++++++-- test/unittests.jl | 1 + test/unittests/test_api.jl | 8 +++++++ test/unittests/test_inputValidation.jl | 29 +++++++++++++++++++++++ test/unittests/test_preconditioner.jl | 29 +++++++++++++---------- test/unittests/test_solve.jl | 32 ++++++++++++++++++++++++++ 10 files changed, 134 insertions(+), 23 deletions(-) diff --git a/src/api.jl b/src/api.jl index 1aad074..5d877b3 100644 --- a/src/api.jl +++ b/src/api.jl @@ -236,8 +236,11 @@ function setup_problem(V::AbstractArray{<:Real}, coords; haskey(SOLVER_NAMES, solver) || throw(ArgumentError("unknown solver :$solver - valid options are :arpack, :krylov, :lobpcg and :lu")) - solver === :arpack && n_eigenvalues + 5 >= length(V) && - throw(ArgumentError("the arpack solver needs n_eigenvalues + 5 < number of grid points ($(length(V)))")) + # :lobpcg transparently falls back to solve_arpack, and :krylov uses the + # same nev, so all three solvers share arpack's nev < N requirement - only + # :lu diagonalizes the full dense matrix and has no such limit + solver in (:arpack, :lobpcg, :krylov) && n_eigenvalues + 5 >= length(V) && + throw(ArgumentError("the $(solver) solver needs n_eigenvalues + 5 < number of grid points ($(length(V)))")) stencil in (3, 5, 7, 9, 11, 13) || throw(ArgumentError("stencil has to be 3, 5, 7, 9, 11 or 13")) stencil_laplace in (3, 5, 7, 9, 11, 13) || throw(ArgumentError("stencil-laplace has to be 3, 5, 7, 9, 11 or 13")) diff --git a/src/preconditioner.jl b/src/preconditioner.jl index dca833b..704cb64 100644 --- a/src/preconditioner.jl +++ b/src/preconditioner.jl @@ -12,13 +12,16 @@ eigenbasis transform per dimension - without ever factorizing a full-size operator, so there is no fill-in. `T̃` equals the true production kinetic operator exactly for 1D problems -(trivially) and for 2D (`buildLaplace_2d`'s Laplacian is separable: dividing +(trivially) and for the 2D **5-point** stencil (`buildLaplace_2d`'s 5-point +Laplacian is a pure row/column "cross" pattern, which is separable: dividing by `2^(dimension-1)`, as `solve()` does, recovers the Kronecker sum exactly). -For 3D, `buildLaplace_3d` uses a more elaborate, non-separable stencil, so -`T̃` is only an APPROXIMATION there; this does not compromise correctness, -only convergence speed, since every LOBPCG result is independently verified -against the true Hamiltonian's residual in `solveWrapper` (and re-solved with -Arpack if that check fails) regardless of how good an approximation `T̃` is. +Every other case - 2D with stencil 3, 7, 9 (the package default) or 11, and +3D with any stencil (`buildLaplace_3d` always uses a more elaborate, +non-separable stencil) - is only an APPROXIMATION; this does not compromise +correctness, only convergence speed, since every LOBPCG result is +independently verified against the true Hamiltonian's residual in +`solveWrapper` (and re-solved with Arpack if that check fails) regardless of +how good an approximation `T̃` is. """ struct KineticPreconditioner Q ::Vector{Matrix{Float64}} # eigenbasis per dimension, d = 1..D diff --git a/src/setupSystem.jl b/src/setupSystem.jl index f54d001..ed50eee 100644 --- a/src/setupSystem.jl +++ b/src/setupSystem.jl @@ -13,6 +13,12 @@ function setupSystem(potential::Potential, system::System) system.reciprocal && !any(system.periodic) && throw(ArgumentError("You have defined a number of k-points - this option is only valid in combination with \"periodic = true\"")) any(system.n_datapoints .< system.stencil) && throw(ArgumentError("The number of datapoints in each dimension has at least to be equal to the stencil size!")) + # solveWrapper only ever builds a real Hamiltonian for lobpcg, so a + # reciprocal (k-point) run has to be rejected here - before any output + # file is written or an existing eigenvalues.dat is removed - rather than + # deep inside solve() on the first k-point + system.reciprocal && system.solver == LOBPCG && throw(ArgumentError("the lobpcg solver supports non-periodic (real symmetric) problems - use arpack or krylov for periodic k-point runs")) + ############################################################### # # # set stencil for laplace and nabla if not defined seperately # diff --git a/src/solve.jl b/src/solve.jl index b62f59c..858799c 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -157,7 +157,9 @@ function solveWrapper(system::System, output::Output, files::Files, Hamiltonian; eigenvalues, eigenvectors = eigenvalues[perm], eigenvectors[:, perm] # verify, don't trust: check the residuals of the eigenpairs that will be - # kept, and escalate or warn if an iterative solver returned loose pairs + # kept, and escalate or warn if an iterative solver returned loose pairs. + # The arpack rescue below is itself re-verified - no path may return + # without passing through this check. if system.solver in (ARPACK, KRYLOV, LOBPCG) residual = max_relative_residual(Hamiltonian, eigenvalues, eigenvectors, output.n_eigenvalues) if system.solver == LOBPCG && residual > 1.0e-6 @@ -165,6 +167,10 @@ function solveWrapper(system::System, output::Output, files::Files, Hamiltonian; eigenvalues, eigenvectors = solve_arpack(Hamiltonian, nev) perm = sortperm(eigenvalues; by = real) eigenvalues, eigenvectors = eigenvalues[perm], eigenvectors[:, perm] + + residual = max_relative_residual(Hamiltonian, eigenvalues, eigenvectors, output.n_eigenvalues) + residual > 1.0e-6 && + @warn "the arpack rescue itself exceeds the residual tolerance; returned eigenpairs may be inaccurate" residual elseif residual > 1.0e-6 @warn "eigenpair residuals are larger than expected" solver = system.solver residual end diff --git a/test/testsets/test_3Dsmoke.jl b/test/testsets/test_3Dsmoke.jl index 5e5763b..207756c 100644 --- a/test/testsets/test_3Dsmoke.jl +++ b/test/testsets/test_3Dsmoke.jl @@ -101,11 +101,29 @@ function test_3Dsmoke() # # ##################################################################### - # the lobpcg solver must reproduce the arpack analytic result + # the lobpcg solver must reproduce the arpack analytic result - and must + # actually do so via lobpcg itself, not via solveWrapper's automatic + # retry/escalation-to-arpack fallback, which would make `r.energies` look + # correct even if lobpcg silently failed on this genuinely (3-fold) + # degenerate cluster and Arpack quietly rescued it. Seed the RNG for + # reproducibility and assert no fallback/escalation warning fired, so + # this test actually proves lobpcg converged rather than merely that the + # pipeline's overall answer is correct. let x = range(xmin, xmax; length = n_points) V = [0.5 * (a^2 + b^2 + c^2) for a in x, b in x, c in x] - r = solve_schrodinger(V, (x, x, x); n_eigenvalues = 4, solver = :lobpcg) (a0, a1) = (0.0032, 0.012) + + Random.seed!(1) + local r + logs, _ = Test.collect_test_logs() do + r = solve_schrodinger(V, (x, x, x); n_eigenvalues = 4, solver = :lobpcg) + end + rescued = any( + l -> occursin("falling back to arpack", l.message) || + occursin("exceed the residual tolerance", l.message), + logs, + ) + @test !rescued @test isapprox(r.energies[1], 1.5; atol = a0) @test all(isapprox.(r.energies[2:4], 2.5; atol = a1)) end diff --git a/test/unittests.jl b/test/unittests.jl index 03c5c90..70e2c8a 100644 --- a/test/unittests.jl +++ b/test/unittests.jl @@ -44,6 +44,7 @@ function unittests() @testset "test solveWrapper" test_solveWrapper() @testset "test lobpcg fallback paths" test_lobpcg_fallback() @testset "test residual warning" test_residual_warning() + @testset "test lobpcg arpack rescue re-verified" test_lobpcg_arpack_rescue_reverified() @testset "test KineticPreconditioner" test_KineticPreconditioner() @testset "API: pipeline equivalence" test_api_equivalence() diff --git a/test/unittests/test_api.jl b/test/unittests/test_api.jl index 634c3cb..e368259 100644 --- a/test/unittests/test_api.jl +++ b/test/unittests/test_api.jl @@ -263,4 +263,12 @@ function test_api_errors() # the nabla stencil has no 13-point variant (checked also on 3D input) @test_throws ArgumentError solve_schrodinger(V3, (z, z, z); stencil_nabla = 13) + + # n_eigenvalues too close to the grid size is rejected up front for every + # solver that shares arpack's nev < N requirement - :arpack directly, + # :lobpcg via its arpack fallback, and :krylov - not just :arpack; before + # this guard, :lobpcg raised an opaque BoundsError deep in solve() instead + for solver in (:arpack, :lobpcg, :krylov) + @test_throws ArgumentError solve_schrodinger(V1, x; n_eigenvalues = length(x), solver = solver) + end end diff --git a/test/unittests/test_inputValidation.jl b/test/unittests/test_inputValidation.jl index 128fde5..431605f 100644 --- a/test/unittests/test_inputValidation.jl +++ b/test/unittests/test_inputValidation.jl @@ -210,6 +210,34 @@ function test_numerov_inputValidation() end end +""" +solver=lobpcg only ever supports non-periodic (real symmetric) problems, but +checkSolver runs before periodicity/k-points are known and so cannot reject +this combination itself - the check has to happen in setupSystem, once +system.reciprocal is set, and before main.jl's unconditional +`rm(files.eigenvalueFileName)`, so a user's existing results survive a +misconfigured re-run rather than being deleted ahead of the error. +""" +function test_lobpcg_periodic_rejected_before_side_effects() + mktempdir() do tmp + cd(tmp) do + write("input.in", """ + potential-file = potential.dat + solver = lobpcg + periodic = true + band-structure = on + k-points = 10 + mass-unit = m_e + """) + write("potential.dat", harmonic_potential_1D()) + write("eigenvalues.dat", "sentinel - must not be deleted by a rejected run") + + @test_throws ArgumentError @suppress Numerov.numerov("input.in") + @test read("eigenvalues.dat", String) == "sentinel - must not be deleted by a rejected run" + end + end +end + function test_inputValidation() test_readInputFile_errors() test_readInputFile_valid() @@ -219,4 +247,5 @@ function test_inputValidation() test_checkSystem_accepted() test_checkOutput() test_numerov_inputValidation() + test_lobpcg_periodic_rejected_before_side_effects() end diff --git a/test/unittests/test_preconditioner.jl b/test/unittests/test_preconditioner.jl index ea5edeb..222f84f 100644 --- a/test/unittests/test_preconditioner.jl +++ b/test/unittests/test_preconditioner.jl @@ -4,16 +4,17 @@ against: the Kronecker SUM of independently-built per-dimension 1D kinetic matrices (`0.5 * -Δ_d / spacing_d²`), which is exactly what the preconditioner is DEFINED to approximate the true kinetic operator with. -For 1D problems this is trivially the production operator itself. For 2D, -`buildLaplace_2d`'s Δ is verified below to equal exactly `2 * (this -separable sum)` - i.e. dividing by `2^(dimension-1)` (as `solve()` does when -assembling the Hamiltonian) recovers this separable model exactly, so the -preconditioner is an exact match there too. For 3D, `buildLaplace_3d` uses a -genuinely non-separable (more elaborate, higher-order) stencil, so this -separable model is only an approximation of the true operator in that case - -that is fine for a preconditioner (only convergence speed depends on the -approximation quality), and final accuracy is independently guaranteed by -`solveWrapper`'s post-solve residual verification, tested elsewhere +For 1D problems this is trivially the production operator itself. For 2D +with the 5-point stencil, `buildLaplace_2d`'s Δ is verified below to equal +exactly `2 * (this separable sum)` - i.e. dividing by `2^(dimension-1)` (as +`solve()` does when assembling the Hamiltonian) recovers this separable +model exactly, so the preconditioner is an exact match there too. Every +other 2D stencil (3, 7, 9, 11) and every 3D stencil use a genuinely +non-separable (more elaborate, higher-order or off-cross) pattern, so this +separable model is only an approximation of the true operator in those +cases - that is fine for a preconditioner (only convergence speed depends on +the approximation quality), and final accuracy is independently guaranteed +by `solveWrapper`'s post-solve residual verification, tested elsewhere (`test_solveWrapper`, `test_lobpcg_fallback`, `test_3Dsmoke`). """ function separable_kinetic_reference(dims::NTuple{D, Int}, potential::Numerov.Potential, system::Numerov.System) where D @@ -44,8 +45,11 @@ via dense Kronecker sums, not by calling any of the preconditioner's own internals - across 1D/2D/3D and a periodic dimension, through all four `ldiv!` dispatches (vector/matrix, in-place/allocating). Also confirms, as a regression guard, that this separable reference exactly equals the true -production kinetic operator for 1D and 2D (division by `2^(dimension-1)` -included), and is only an approximation for 3D (see module docstring above). +production kinetic operator for 1D and for 2D's 5-point stencil (division by +`2^(dimension-1)` included), and is only an approximation for every other 2D +stencil and for 3D (see module docstring above) - the 2D stencil=9 case below +guards specifically against re-introducing the "2D is always exact" overclaim +this docstring used to make (the package's default stencil is 9, not 5). """ function test_KineticPreconditioner() Random.seed!(42) @@ -54,6 +58,7 @@ function test_KineticPreconditioner() (dims = (12,), periodic = false, stencil = 9, exact_vs_production = true), (dims = (12,), periodic = true, stencil = 9, exact_vs_production = true), (dims = (6, 6), periodic = false, stencil = 5, exact_vs_production = true), + (dims = (10, 10), periodic = false, stencil = 9, exact_vs_production = false), (dims = (5, 5, 5), periodic = false, stencil = 5, exact_vs_production = false), ) diff --git a/test/unittests/test_solve.jl b/test/unittests/test_solve.jl index 6895a38..a0bdb17 100644 --- a/test/unittests/test_solve.jl +++ b/test/unittests/test_solve.jl @@ -153,3 +153,35 @@ function test_residual_warning() # path - just confirm the residual it reports is indeed loose @test Numerov.max_relative_residual(H, eigenvalues, eigenvectors, n) > 1.0e-6 end + +""" +The lobpcg->arpack escalation (triggered by a loose residual, not by an +exception) must itself be re-verified rather than returned on trust. Force +lobpcg to under-converge (`lobpcg_maxiter = 1`, as in `test_lobpcg_fallback`) +on a Hamiltonian whose extreme diagonal dynamic range also makes the arpack +rescue itself land above the residual tolerance - `solve_arpack`'s shift +heuristic scales with the largest diagonal entry, so a single huge spike +degrades shift-invert accuracy for the low-lying eigenpairs actually wanted. +This must produce a second, distinct warning naming the rescue result +itself as suspect, not silence. +""" +function test_lobpcg_arpack_rescue_reverified() + Random.seed!(3) + + N, n = 40, 4 + Δ = spdiagm(-1 => -ones(N - 1), 0 => 2 * ones(N), 1 => -ones(N - 1)) + V = collect(range(0.1, 2.0; length = N)) + V[20] = 1.0e10 + H = Δ + spdiagm(0 => V) + + system, output, files = make_solver_structs(Numerov.LOBPCG, n) + local eigenvalues, eigenvectors + @test_logs (:warn, r"re-solving with arpack") (:warn, r"arpack rescue itself exceeds") match_mode = :any begin + eigenvalues, eigenvectors = Numerov.solveWrapper(system, output, files, H; lobpcg_maxiter = 1) + end + + # the rescue result is still returned (nothing better to fall back to) - + # confirm it is indeed the loose result the warning describes, not a + # spuriously-triggered warning on an otherwise-fine result + @test Numerov.max_relative_residual(H, eigenvalues, eigenvectors, n) > 1.0e-6 +end