From c5606b5e1c6b2455541b68e26db755d15337ef0e Mon Sep 17 00:00:00 2001 From: Jash Date: Fri, 3 Jul 2026 21:09:11 +0530 Subject: [PATCH 1/7] Add EigenvalueProblem support with dense, Arpack, ArnoldiMethod, KrylovKit, and JacobiDavidson backends --- Project.toml | 14 ++- ext/LinearSolveArnoldiMethodExt.jl | 49 ++++++++++ ext/LinearSolveArpackExt.jl | 30 ++++++ ext/LinearSolveJacobiDavidsonExt.jl | 65 +++++++++++++ ext/LinearSolveKrylovKitExt.jl | 45 +++++++++ src/LinearSolve.jl | 5 + src/eigenvalue.jl | 136 ++++++++++++++++++++++++++++ test/Core/eigenvalue.jl | 82 +++++++++++++++++ test/runtests.jl | 1 + 9 files changed, 426 insertions(+), 1 deletion(-) create mode 100644 ext/LinearSolveArnoldiMethodExt.jl create mode 100644 ext/LinearSolveArpackExt.jl create mode 100644 ext/LinearSolveJacobiDavidsonExt.jl create mode 100644 src/eigenvalue.jl create mode 100644 test/Core/eigenvalue.jl diff --git a/Project.toml b/Project.toml index b1af7bf35..7ea45d821 100644 --- a/Project.toml +++ b/Project.toml @@ -33,6 +33,8 @@ StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" [weakdeps] AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" +Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" +ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d" BandedMatrices = "aae01518-5342-5314-be14-df237901396f" BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0" cuSOLVER = "887afef0-6a32-4de5-add4-7827692ba8fc" @@ -49,6 +51,7 @@ Ginkgo = "4c8bd3c9-ead9-4b5e-a625-08f1338ba0ec" HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" HSL = "34c5aeac-e683-54a6-a0e9-6e0fdc586c50" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" +JacobiDavidson = "11c68b98-9c9b-11e8-267b-bbb95576cead" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LAPACK_jll = "51474c39-65e3-53ba-86ba-03b1b862ec14" @@ -72,6 +75,8 @@ blis_jll = "6136c539-28a5-5bf0-87cc-b183200dce32" [extensions] LinearSolveAMDGPUExt = "AMDGPU" LinearSolveAlgebraicMultigridExt = "AlgebraicMultigrid" +LinearSolveArpackExt = "Arpack" +LinearSolveArnoldiMethodExt = "ArnoldiMethod" LinearSolveBLISExt = ["blis_jll", "LAPACK_jll"] LinearSolveBandedMatricesExt = "BandedMatrices" LinearSolveBlockDiagonalsExt = "BlockDiagonals" @@ -89,6 +94,7 @@ LinearSolveGinkgoExt = ["Ginkgo", "SparseArrays"] LinearSolveHYPREExt = "HYPRE" LinearSolveHSLExt = ["HSL", "SparseArrays"] LinearSolveIterativeSolversExt = "IterativeSolvers" +LinearSolveJacobiDavidsonExt = "JacobiDavidson" LinearSolveKernelAbstractionsExt = "KernelAbstractions" LinearSolveKrylovKitExt = "KrylovKit" LinearSolveMetalExt = "Metal" @@ -111,6 +117,8 @@ LinearSolveSpecializingFactorizationsExt = "SpecializingFactorizations" AMD = "0.5" AMDGPU = "1.2, 2" AlgebraicMultigrid = "1, 2" +Arpack = "0.5" +ArnoldiMethod = "0.4" ArrayInterface = "7.19" BandedMatrices = "1.8" BlockDiagonals = "0.2" @@ -136,6 +144,7 @@ HYPRE = "1.7" HSL = "0.5" InteractiveUtils = "1.10" IterativeSolvers = "0.9.4" +JacobiDavidson = "0.1" KernelAbstractions = "0.9.30" Krylov = "0.10" KrylovKit = "0.10" @@ -186,6 +195,8 @@ julia = "1.10" [extras] AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" +Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" +ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d" BandedMatrices = "aae01518-5342-5314-be14-df237901396f" BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0" CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" @@ -197,6 +208,7 @@ FixedSizeArrays = "3821ddf9-e5b5-40d5-8e25-6813ab96b5e2" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" +JacobiDavidson = "11c68b98-9c9b-11e8-267b-bbb95576cead" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" KrylovPreconditioners = "45d422c2-293f-44ce-8315-2cb988662dec" MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5" @@ -215,4 +227,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["AlgebraicMultigrid", "BandedMatrices", "BlockDiagonals", "CliqueTrees", "ComponentArrays", "FastAlmostBandedMatrices", "FastLapackInterface", "FiniteDiff", "FixedSizeArrays", "ForwardDiff", "InteractiveUtils", "IterativeSolvers", "KrylovKit", "KrylovPreconditioners", "MultiFloats", "Pkg", "PureUMFPACK", "Random", "RecursiveFactorization", "STRUMPACK_jll", "SafeTestsets", "SciMLTesting", "SparseArrays", "Sparspak", "SpecializingFactorizations", "StaticArrays", "Test", "Zygote"] +test = ["AlgebraicMultigrid", "Arpack", "ArnoldiMethod", "BandedMatrices", "BlockDiagonals", "CliqueTrees", "ComponentArrays", "FastAlmostBandedMatrices", "FastLapackInterface", "FiniteDiff", "FixedSizeArrays", "ForwardDiff", "InteractiveUtils", "IterativeSolvers", "JacobiDavidson", "KrylovKit", "KrylovPreconditioners", "MultiFloats", "Pkg", "PureUMFPACK", "Random", "RecursiveFactorization", "STRUMPACK_jll", "SafeTestsets", "SciMLTesting", "SparseArrays", "Sparspak", "SpecializingFactorizations", "StaticArrays", "Test", "Zygote"] diff --git a/ext/LinearSolveArnoldiMethodExt.jl b/ext/LinearSolveArnoldiMethodExt.jl new file mode 100644 index 000000000..1bc5b3171 --- /dev/null +++ b/ext/LinearSolveArnoldiMethodExt.jl @@ -0,0 +1,49 @@ +module LinearSolveArnoldiMethodExt + +using LinearAlgebra +using LinearSolve +import ArnoldiMethod: partialschur, partialeigen +using SciMLBase: SciMLBase, ReturnCode + +function SciMLBase.solve( + prob::LinearSolve.EigenvalueProblem, + alg::LinearSolve.ArnoldiMethodJL, + args...; kwargs...) + prob.B === nothing || + error("ArnoldiMethod backend currently supports standard eigenvalue problems only.") + nev = LinearSolve.default_nev(prob) + which = prob.sigma === nothing ? prob.which : :LM + A = prob.sigma === nothing ? prob.A : _shift_invert_operator(prob.A, prob.sigma) + kw = (; nev, which, prob.kwargs..., alg.kwargs..., kwargs...) + decomp, history = partialschur(A, alg.args...; kw...) + values, vectors = partialeigen(decomp) + if prob.sigma !== nothing + values = prob.sigma .+ inv.(values) + end + values, vectors = LinearSolve._select_eigenpairs(values, vectors, nev, prob.which, prob.sigma) + retcode = history.converged ? ReturnCode.Success : ReturnCode.ConvergenceFailure + return LinearSolve.build_eigenvalue_solution( + prob, alg, values, vectors; retcode, stats = history) +end + +function _shift_invert_operator(A, sigma) + F = factorize(A - sigma * I) + T = promote_type(eltype(A), typeof(sigma)) + return ShiftInvertMap{typeof(F), T}(F, size(A, 1)) +end + +struct ShiftInvertMap{F, T} + F::F + n::Int +end + +Base.size(A::ShiftInvertMap) = (A.n, A.n) +Base.size(A::ShiftInvertMap, dim::Integer) = dim <= 2 ? A.n : 1 +Base.eltype(::Type{<:ShiftInvertMap{F, T}}) where {F, T} = T + +function LinearAlgebra.mul!(y, A::ShiftInvertMap, x) + copyto!(y, A.F \ x) + return y +end + +end diff --git a/ext/LinearSolveArpackExt.jl b/ext/LinearSolveArpackExt.jl new file mode 100644 index 000000000..31fa5d981 --- /dev/null +++ b/ext/LinearSolveArpackExt.jl @@ -0,0 +1,30 @@ +module LinearSolveArpackExt + +using LinearSolve +using Arpack +using SciMLBase: SciMLBase, ReturnCode + +function SciMLBase.solve( + prob::LinearSolve.EigenvalueProblem, + alg::LinearSolve.ArpackJL, + args...; kwargs...) + nev = LinearSolve.default_nev(prob) + base = (; nev, which = prob.which) + if prob.sigma !== nothing + base = (; base..., sigma = prob.sigma) + end + kw = (; base..., prob.kwargs..., alg.kwargs..., kwargs...) + # `Arpack.eigs` takes the generalized-problem matrix `B` positionally, not + # as a keyword argument. + values, vectors, nconv, niter, nmult, resid = if prob.B === nothing + Arpack.eigs(prob.A, alg.args...; kw...) + else + Arpack.eigs(prob.A, prob.B, alg.args...; kw...) + end + retcode = nconv >= length(values) ? ReturnCode.Success : ReturnCode.ConvergenceFailure + stats = (; nconv, niter, nmult) + return LinearSolve.build_eigenvalue_solution( + prob, alg, values, vectors; retcode, resid, stats) +end + +end diff --git a/ext/LinearSolveJacobiDavidsonExt.jl b/ext/LinearSolveJacobiDavidsonExt.jl new file mode 100644 index 000000000..bbc3a202f --- /dev/null +++ b/ext/LinearSolveJacobiDavidsonExt.jl @@ -0,0 +1,65 @@ +module LinearSolveJacobiDavidsonExt + +using LinearSolve +using LinearAlgebra +using JacobiDavidson +using SciMLBase: SciMLBase, ReturnCode + +function SciMLBase.solve( + prob::LinearSolve.EigenvalueProblem, + alg::LinearSolve.JacobiDavidsonJL, + args...; kwargs...) + # JacobiDavidson.jl's `jdqz` (generalized solver) is broken upstream: it + # references an undefined `verbose` variable that is absent from its + # signature. Until that is fixed, only the standard solver `jdqr` is wired + # up here; point users at the backends that do support generalized problems. + prob.B === nothing || + error("The JacobiDavidson backend currently supports standard eigenvalue problems only. Use `ArpackJL()` or `KrylovKitEigen()` for generalized problems.") + + n = size(prob.A, 2) + nev = LinearSolve.default_nev(prob) + target = _jd_target(prob) + # Search-subspace bounds, capped at the problem size. Users may override + # `subspace_dimensions` (and any other jdqr keyword) via the algorithm. + hi = min(max(2 * nev + 10, 20), n) + lo = min(max(nev + 2, 8), hi) + defaults = (; pairs = nev, target = target, subspace_dimensions = lo:hi) + kw = (; defaults..., prob.kwargs..., alg.kwargs..., kwargs...) + + out = JacobiDavidson.jdqr(prob.A, alg.args...; kw...) + values, vectors = _jd_standard_pairs(prob.A, out[1]) + + values, vectors = LinearSolve._select_eigenpairs( + values, vectors, nev, prob.which, prob.sigma) + retcode = length(values) >= nev ? ReturnCode.Success : ReturnCode.ConvergenceFailure + return LinearSolve.build_eigenvalue_solution( + prob, alg, values, vectors; retcode, stats = out[end]) +end + +# Map the problem's spectral selector onto a JacobiDavidson `Target`. A supplied +# `sigma` is the natural interior target (`Near`), which is Jacobi-Davidson's +# strength; otherwise the `which` symbol selects an extremal target. +function _jd_target(prob) + prob.sigma !== nothing && return JacobiDavidson.Near(ComplexF64(prob.sigma)) + w = prob.which + w === :LM ? JacobiDavidson.LargestMagnitude(0.0 + 0.0im) : + w === :SM ? JacobiDavidson.SmallestMagnitude(0.0 + 0.0im) : + w === :LR ? JacobiDavidson.LargestRealPart(0.0 + 0.0im) : + w === :SR ? JacobiDavidson.SmallestRealPart(0.0 + 0.0im) : + w === :LI ? JacobiDavidson.LargestImaginaryPart(0.0 + 0.0im) : + w === :SI ? JacobiDavidson.SmallestImaginaryPart(0.0 + 0.0im) : + throw(ArgumentError("unsupported `which = $(w)` for JacobiDavidson; expected one of :LM, :SM, :LR, :SR, :LI, :SI, or pass `sigma`")) +end + +# jdqr yields a partial Schur decomposition `A*Q = Q*R`. Eigenpairs are recovered +# from the small projected `R = Q'AQ`: if `R*y = λ*y` then `A*(Q*y) = λ*(Q*y)`. +function _jd_standard_pairs(A, pschur) + T = complex(float(eltype(A))) + k = length(pschur.values) + k == 0 && return (T[], Matrix{T}(undef, size(A, 1), 0)) + Q = pschur.Q[:, 1:k] + F = eigen(Q' * (A * Q)) + return (F.values, Q * F.vectors) +end + +end diff --git a/ext/LinearSolveKrylovKitExt.jl b/ext/LinearSolveKrylovKitExt.jl index 3731b3a8b..9244c5463 100644 --- a/ext/LinearSolveKrylovKitExt.jl +++ b/ext/LinearSolveKrylovKitExt.jl @@ -61,4 +61,49 @@ end LinearSolve.update_tolerances_internal!(cache, alg::KrylovKitJL, atol, rtol) = nothing +function SciMLBase.solve( + prob::LinearSolve.EigenvalueProblem, + alg::LinearSolve.KrylovKitEigen, + args...; kwargs...) + nev = LinearSolve.default_nev(prob) + which = prob.which + kw = (; prob.kwargs..., alg.kwargs..., kwargs...) + values, vectors, info = if prob.sigma !== nothing + _shift_invert_eigsolve(prob, alg, nev, kw) + elseif prob.B === nothing + KrylovKit.eigsolve(prob.A, alg.args..., nev, which; kw...) + else + KrylovKit.geneigsolve((prob.A, prob.B), alg.args..., nev, which; kw...) + end + if prob.sigma !== nothing + values = prob.sigma .+ inv.(values) + end + vecmat = reduce(hcat, vectors) + values, vecmat = LinearSolve._select_eigenpairs(values, vecmat, nev, prob.which, prob.sigma) + retcode = info.converged >= length(values) ? ReturnCode.Success : ReturnCode.ConvergenceFailure + return LinearSolve.build_eigenvalue_solution( + prob, alg, values, vecmat; retcode, resid = info.normres, stats = info) +end + +function _shift_invert_eigsolve(prob, alg, nev, kw) + A, B, sigma = prob.A, prob.B, prob.sigma + T = isnothing(B) ? promote_type(eltype(A), typeof(sigma)) : + promote_type(eltype(A), eltype(B), typeof(sigma)) + if isnothing(B) + F = factorize(A - sigma * I) + op = x -> F \ x + elseif B isa LinearAlgebra.UniformScaling + F = factorize(A - sigma * B) + op = x -> F \ (B.λ * x) + else + F = factorize(A - sigma * B) + op = x -> F \ (B * x) + end + if isempty(alg.args) + return KrylovKit.eigsolve(op, size(A, 2), nev, :LM, T; kw...) + else + return KrylovKit.eigsolve(op, alg.args..., nev, :LM; kw...) + end +end + end diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 292fbd468..dd02ed9ac 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -410,6 +410,7 @@ function defaultalg_symbol end include("verbosity.jl") include("blas_logging.jl") include("generic_lufact.jl") +include("eigenvalue.jl") include("common.jl") include("extension_algs.jl") include("factorization.jl") @@ -573,4 +574,8 @@ export LinearSolveAdjoint export LinearVerbosity +export EigenvalueProblem, EigenvalueSolution, AbstractEigenvalueAlgorithm, + DenseEigen, ArpackJL, ArnoldiMethod, ArnoldiMethodJL, KrylovKitEigen, + JacobiDavidsonJL + end diff --git a/src/eigenvalue.jl b/src/eigenvalue.jl new file mode 100644 index 000000000..14af57c33 --- /dev/null +++ b/src/eigenvalue.jl @@ -0,0 +1,136 @@ +""" + EigenvalueProblem(A[, B]; nev = nothing, which = :LM, sigma = nothing, u0 = nothing, p = SciMLBase.NullParameters()) + +Define a standard or generalized eigenvalue problem. + +The standard problem is ``A v = lambda v``. If `B` is supplied, the generalized +problem is ``A v = lambda B v``. `nev` requests a subset of eigenpairs, `which` +selects the part of the spectrum (`:LM`, `:SM`, `:LR`, `:SR`, `:LI`, `:SI`), +and `sigma` requests eigenvalues nearest to the shift. +""" +struct EigenvalueProblem{AType, BType, NevType, WhichType, SigmaType, U0Type, PType, KType} + A::AType + B::BType + nev::NevType + which::WhichType + sigma::SigmaType + u0::U0Type + p::PType + kwargs::KType +end + +function EigenvalueProblem(A, B = nothing, p = SciMLBase.NullParameters(); + nev = nothing, which = :LM, sigma = nothing, u0 = nothing, kwargs...) + return EigenvalueProblem{ + typeof(A), typeof(B), typeof(nev), typeof(which), typeof(sigma), + typeof(u0), typeof(p), typeof(kwargs), + }(A, B, nev, which, sigma, u0, p, kwargs) +end + +abstract type AbstractEigenvalueAlgorithm <: SciMLBase.AbstractLinearAlgorithm end + +struct DenseEigen <: AbstractEigenvalueAlgorithm end + +struct ArpackJL{A, K} <: AbstractEigenvalueAlgorithm + args::A + kwargs::K +end +ArpackJL(args...; kwargs...) = ArpackJL(args, kwargs) + +struct ArnoldiMethodJL{A, K} <: AbstractEigenvalueAlgorithm + args::A + kwargs::K +end +ArnoldiMethod(args...; kwargs...) = ArnoldiMethodJL(args, kwargs) + +struct KrylovKitEigen{A, K} <: AbstractEigenvalueAlgorithm + args::A + kwargs::K +end +KrylovKitEigen(args...; kwargs...) = KrylovKitEigen(args, kwargs) + +struct JacobiDavidsonJL{A, K} <: AbstractEigenvalueAlgorithm + args::A + kwargs::K +end +JacobiDavidsonJL(args...; kwargs...) = JacobiDavidsonJL(args, kwargs) + +struct EigenvalueSolution{T, N, U, V, P, A, R, S} <: SciMLBase.AbstractNoTimeSolution{T, N} + u::U + vectors::V + prob::P + alg::A + retcode::ReturnCode.T + resid::R + stats::S +end + +function build_eigenvalue_solution(prob, alg, values, vectors; + retcode = ReturnCode.Success, resid = nothing, stats = nothing) + T = eltype(eltype(values)) + N = length((size(values)...,)) + return EigenvalueSolution{ + T, N, typeof(values), typeof(vectors), typeof(prob), typeof(alg), + typeof(resid), typeof(stats), + }(values, vectors, prob, alg, retcode, resid, stats) +end + +SciMLBase.solve(prob::EigenvalueProblem, args...; kwargs...) = + solve(prob, nothing, args...; kwargs...) + +function SciMLBase.solve(prob::EigenvalueProblem, ::Nothing, args...; kwargs...) + return solve(prob, DenseEigen(), args...; kwargs...) +end + +function SciMLBase.solve(prob::EigenvalueProblem, alg::DenseEigen, args...; kwargs...) + kw = (; prob.kwargs..., kwargs...) + F = if isnothing(prob.B) + LinearAlgebra.eigen(prob.A; kw...) + elseif prob.B isa UniformScaling + LinearAlgebra.eigen(prob.A / prob.B.λ; kw...) + else + LinearAlgebra.eigen(prob.A, prob.B; kw...) + end + values, vectors = _select_eigenpairs(F.values, F.vectors, prob.nev, prob.which, prob.sigma) + return build_eigenvalue_solution(prob, alg, values, vectors) +end + +function SciMLBase.solve( + prob::EigenvalueProblem, alg::AbstractEigenvalueAlgorithm, args...; kwargs...) + error("The eigenvalue backend $(typeof(alg)) is not available. Load its package before solving with this algorithm.") +end + +function default_nev(prob::EigenvalueProblem) + n = size(prob.A, 2) + # Only the iterative backends call this; requesting the full dimension `n` + # is invalid/degenerate for them, so default to a small subset. + return prob.nev === nothing ? min(n, 6) : prob.nev +end + +function _select_eigenpairs(values, vectors, nev, which, sigma) + nvals = length(values) + howmany = nev === nothing ? nvals : min(nev, nvals) + ord = _eigenvalue_order(values, which, sigma) + idxs = ord[1:howmany] + return values[idxs], vectors[:, idxs] +end + +function _eigenvalue_order(values, which, sigma) + if sigma !== nothing + return sortperm(abs.(values .- sigma)) + elseif which === :LM + return sortperm(abs.(values); rev = true) + elseif which === :SM + return sortperm(abs.(values)) + elseif which === :LR + return sortperm(real.(values); rev = true) + elseif which === :SR + return sortperm(real.(values)) + elseif which === :LI + return sortperm(imag.(values); rev = true) + elseif which === :SI + return sortperm(imag.(values)) + else + throw(ArgumentError("unsupported eigenvalue selector `which = $which`; expected one of :LM, :SM, :LR, :SR, :LI, or :SI")) + end +end diff --git a/test/Core/eigenvalue.jl b/test/Core/eigenvalue.jl new file mode 100644 index 000000000..ac5a6856f --- /dev/null +++ b/test/Core/eigenvalue.jl @@ -0,0 +1,82 @@ +using LinearSolve +using LinearAlgebra +using Test +using Arpack +using ArnoldiMethod +using KrylovKit +using JacobiDavidson + +A = Diagonal([1.0, 2.0, 3.0, 4.0]) + +prob = EigenvalueProblem(A) +sol = solve(prob) +@test sol.u == [4.0, 3.0, 2.0, 1.0] +@test sol.vectors * Diagonal(sol.u) ≈ A * sol.vectors +@test sol.retcode === ReturnCode.Success + +prob_nev = EigenvalueProblem(A; nev = 2, which = :SM) +sol_nev = solve(prob_nev) +@test sol_nev.u == [1.0, 2.0] + +prob_sigma = EigenvalueProblem(A; nev = 2, sigma = 2.2) +sol_sigma = solve(prob_sigma) +@test sol_sigma.u == [2.0, 3.0] + +B = Diagonal(fill(2.0, 4)) +prob_gen = EigenvalueProblem(A, B; nev = 2, which = :LR) +sol_gen = solve(prob_gen) +@test sol_gen.u == [2.0, 1.5] + +prob_uniform = EigenvalueProblem(A, 2I; nev = 2, which = :LR) +sol_uniform = solve(prob_uniform) +@test sol_uniform.u == [2.0, 1.5] + +A_backend = Diagonal(1.0:8.0) +B_backend = Diagonal(fill(2.0, 8)) + +sol_arpack = solve(EigenvalueProblem(Matrix(A_backend); nev = 2), ArpackJL()) +@test sol_arpack.u ≈ [8.0, 7.0] + +sol_arpack_gen = solve( + EigenvalueProblem(Matrix(A_backend), Matrix(B_backend); nev = 2), ArpackJL()) +@test sol_arpack_gen.u ≈ [4.0, 3.5] + +sol_arnoldi = solve(EigenvalueProblem(A_backend; nev = 2), LinearSolve.ArnoldiMethod()) +@test sol_arnoldi.u ≈ [8.0, 7.0] + +sol_arnoldi_default = solve(EigenvalueProblem(A_backend), LinearSolve.ArnoldiMethod()) +@test length(sol_arnoldi_default.u) == 6 + +sol_arnoldi_shift = solve( + EigenvalueProblem(A_backend; nev = 2, sigma = 3.2), LinearSolve.ArnoldiMethod()) +@test sol_arnoldi_shift.u ≈ [3.0, 4.0] + +sol_krylovkit = solve(EigenvalueProblem(A_backend; nev = 2), KrylovKitEigen()) +@test sol_krylovkit.u ≈ [8.0, 7.0] + +sol_krylovkit_shift = solve( + EigenvalueProblem(Matrix(A_backend); nev = 2, sigma = 3.2), KrylovKitEigen()) +@test sol_krylovkit_shift.u ≈ [3.0, 4.0] + +sol_krylovkit_gen_shift = solve( + EigenvalueProblem(Matrix(A_backend), Matrix(B_backend); nev = 2, sigma = 1.6), + KrylovKitEigen()) +@test sol_krylovkit_gen_shift.u ≈ [1.5, 2.0] + +# Jacobi-Davidson is a target/interior method: it finds the eigenvalues nearest +# the target shift. `sigma` gives the interior target directly, and `:SM`/`:SR` +# work with the implicit zero shift. (`:LM`/`:LR` need a `sigma` guess.) +A_jd = Matrix(Diagonal(Float64.(1:30))) + +sol_jd = solve(EigenvalueProblem(A_jd; nev = 2, sigma = 10.3), JacobiDavidsonJL()) +@test sort(real(sol_jd.u)) ≈ [10.0, 11.0] +@test sol_jd.retcode === ReturnCode.Success +@test norm(A_jd * sol_jd.vectors - sol_jd.vectors * Diagonal(sol_jd.u)) < 1e-6 + +sol_jd_sm = solve(EigenvalueProblem(A_jd; nev = 1, which = :SM), JacobiDavidsonJL()) +@test real(sol_jd_sm.u[1]) ≈ 1.0 + +# Generalized problems are not supported by the JacobiDavidson backend (upstream +# `jdqz` is broken); it should error rather than silently misbehave. +@test_throws ErrorException solve( + EigenvalueProblem(A_jd, Matrix(Diagonal(fill(2.0, 30))); nev = 1), JacobiDavidsonJL()) diff --git a/test/runtests.jl b/test/runtests.jl index 118343b96..2a95218cf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -62,6 +62,7 @@ else default = "All", core = function () @time @safetestset "Basic Tests" include("Core/basictests.jl") + @time @safetestset "EigenvalueProblem" include("Core/eigenvalue.jl") @time @safetestset "Return codes" include("Core/retcodes.jl") @time @safetestset "Re-solve" include("Core/resolve.jl") @time @safetestset "Zero Initialization Tests" include("Core/zeroinittests.jl") From 1ebe83df26ce78dcfd83eaedd4a808f66cf68bae Mon Sep 17 00:00:00 2001 From: Jash Date: Fri, 3 Jul 2026 21:29:34 +0530 Subject: [PATCH 2/7] fixed runic failures --- ext/LinearSolveArnoldiMethodExt.jl | 6 ++++-- ext/LinearSolveArpackExt.jl | 6 ++++-- ext/LinearSolveJacobiDavidsonExt.jl | 23 +++++++++++++---------- ext/LinearSolveKrylovKitExt.jl | 6 ++++-- src/eigenvalue.jl | 15 ++++++++++----- test/Core/eigenvalue.jl | 17 +++++++++++------ 6 files changed, 46 insertions(+), 27 deletions(-) diff --git a/ext/LinearSolveArnoldiMethodExt.jl b/ext/LinearSolveArnoldiMethodExt.jl index 1bc5b3171..f222742a5 100644 --- a/ext/LinearSolveArnoldiMethodExt.jl +++ b/ext/LinearSolveArnoldiMethodExt.jl @@ -8,7 +8,8 @@ using SciMLBase: SciMLBase, ReturnCode function SciMLBase.solve( prob::LinearSolve.EigenvalueProblem, alg::LinearSolve.ArnoldiMethodJL, - args...; kwargs...) + args...; kwargs... + ) prob.B === nothing || error("ArnoldiMethod backend currently supports standard eigenvalue problems only.") nev = LinearSolve.default_nev(prob) @@ -23,7 +24,8 @@ function SciMLBase.solve( values, vectors = LinearSolve._select_eigenpairs(values, vectors, nev, prob.which, prob.sigma) retcode = history.converged ? ReturnCode.Success : ReturnCode.ConvergenceFailure return LinearSolve.build_eigenvalue_solution( - prob, alg, values, vectors; retcode, stats = history) + prob, alg, values, vectors; retcode, stats = history + ) end function _shift_invert_operator(A, sigma) diff --git a/ext/LinearSolveArpackExt.jl b/ext/LinearSolveArpackExt.jl index 31fa5d981..4acfc1e42 100644 --- a/ext/LinearSolveArpackExt.jl +++ b/ext/LinearSolveArpackExt.jl @@ -7,7 +7,8 @@ using SciMLBase: SciMLBase, ReturnCode function SciMLBase.solve( prob::LinearSolve.EigenvalueProblem, alg::LinearSolve.ArpackJL, - args...; kwargs...) + args...; kwargs... + ) nev = LinearSolve.default_nev(prob) base = (; nev, which = prob.which) if prob.sigma !== nothing @@ -24,7 +25,8 @@ function SciMLBase.solve( retcode = nconv >= length(values) ? ReturnCode.Success : ReturnCode.ConvergenceFailure stats = (; nconv, niter, nmult) return LinearSolve.build_eigenvalue_solution( - prob, alg, values, vectors; retcode, resid, stats) + prob, alg, values, vectors; retcode, resid, stats + ) end end diff --git a/ext/LinearSolveJacobiDavidsonExt.jl b/ext/LinearSolveJacobiDavidsonExt.jl index bbc3a202f..9ae23f8af 100644 --- a/ext/LinearSolveJacobiDavidsonExt.jl +++ b/ext/LinearSolveJacobiDavidsonExt.jl @@ -8,7 +8,8 @@ using SciMLBase: SciMLBase, ReturnCode function SciMLBase.solve( prob::LinearSolve.EigenvalueProblem, alg::LinearSolve.JacobiDavidsonJL, - args...; kwargs...) + args...; kwargs... + ) # JacobiDavidson.jl's `jdqz` (generalized solver) is broken upstream: it # references an undefined `verbose` variable that is absent from its # signature. Until that is fixed, only the standard solver `jdqr` is wired @@ -30,10 +31,12 @@ function SciMLBase.solve( values, vectors = _jd_standard_pairs(prob.A, out[1]) values, vectors = LinearSolve._select_eigenpairs( - values, vectors, nev, prob.which, prob.sigma) + values, vectors, nev, prob.which, prob.sigma + ) retcode = length(values) >= nev ? ReturnCode.Success : ReturnCode.ConvergenceFailure return LinearSolve.build_eigenvalue_solution( - prob, alg, values, vectors; retcode, stats = out[end]) + prob, alg, values, vectors; retcode, stats = out[end] + ) end # Map the problem's spectral selector onto a JacobiDavidson `Target`. A supplied @@ -42,13 +45,13 @@ end function _jd_target(prob) prob.sigma !== nothing && return JacobiDavidson.Near(ComplexF64(prob.sigma)) w = prob.which - w === :LM ? JacobiDavidson.LargestMagnitude(0.0 + 0.0im) : - w === :SM ? JacobiDavidson.SmallestMagnitude(0.0 + 0.0im) : - w === :LR ? JacobiDavidson.LargestRealPart(0.0 + 0.0im) : - w === :SR ? JacobiDavidson.SmallestRealPart(0.0 + 0.0im) : - w === :LI ? JacobiDavidson.LargestImaginaryPart(0.0 + 0.0im) : - w === :SI ? JacobiDavidson.SmallestImaginaryPart(0.0 + 0.0im) : - throw(ArgumentError("unsupported `which = $(w)` for JacobiDavidson; expected one of :LM, :SM, :LR, :SR, :LI, :SI, or pass `sigma`")) + return w === :LM ? JacobiDavidson.LargestMagnitude(0.0 + 0.0im) : + w === :SM ? JacobiDavidson.SmallestMagnitude(0.0 + 0.0im) : + w === :LR ? JacobiDavidson.LargestRealPart(0.0 + 0.0im) : + w === :SR ? JacobiDavidson.SmallestRealPart(0.0 + 0.0im) : + w === :LI ? JacobiDavidson.LargestImaginaryPart(0.0 + 0.0im) : + w === :SI ? JacobiDavidson.SmallestImaginaryPart(0.0 + 0.0im) : + throw(ArgumentError("unsupported `which = $(w)` for JacobiDavidson; expected one of :LM, :SM, :LR, :SR, :LI, :SI, or pass `sigma`")) end # jdqr yields a partial Schur decomposition `A*Q = Q*R`. Eigenpairs are recovered diff --git a/ext/LinearSolveKrylovKitExt.jl b/ext/LinearSolveKrylovKitExt.jl index 9244c5463..b81372079 100644 --- a/ext/LinearSolveKrylovKitExt.jl +++ b/ext/LinearSolveKrylovKitExt.jl @@ -64,7 +64,8 @@ LinearSolve.update_tolerances_internal!(cache, alg::KrylovKitJL, atol, rtol) = n function SciMLBase.solve( prob::LinearSolve.EigenvalueProblem, alg::LinearSolve.KrylovKitEigen, - args...; kwargs...) + args...; kwargs... + ) nev = LinearSolve.default_nev(prob) which = prob.which kw = (; prob.kwargs..., alg.kwargs..., kwargs...) @@ -82,7 +83,8 @@ function SciMLBase.solve( values, vecmat = LinearSolve._select_eigenpairs(values, vecmat, nev, prob.which, prob.sigma) retcode = info.converged >= length(values) ? ReturnCode.Success : ReturnCode.ConvergenceFailure return LinearSolve.build_eigenvalue_solution( - prob, alg, values, vecmat; retcode, resid = info.normres, stats = info) + prob, alg, values, vecmat; retcode, resid = info.normres, stats = info + ) end function _shift_invert_eigsolve(prob, alg, nev, kw) diff --git a/src/eigenvalue.jl b/src/eigenvalue.jl index 14af57c33..23843cbc9 100644 --- a/src/eigenvalue.jl +++ b/src/eigenvalue.jl @@ -19,8 +19,10 @@ struct EigenvalueProblem{AType, BType, NevType, WhichType, SigmaType, U0Type, PT kwargs::KType end -function EigenvalueProblem(A, B = nothing, p = SciMLBase.NullParameters(); - nev = nothing, which = :LM, sigma = nothing, u0 = nothing, kwargs...) +function EigenvalueProblem( + A, B = nothing, p = SciMLBase.NullParameters(); + nev = nothing, which = :LM, sigma = nothing, u0 = nothing, kwargs... + ) return EigenvalueProblem{ typeof(A), typeof(B), typeof(nev), typeof(which), typeof(sigma), typeof(u0), typeof(p), typeof(kwargs), @@ -65,8 +67,10 @@ struct EigenvalueSolution{T, N, U, V, P, A, R, S} <: SciMLBase.AbstractNoTimeSol stats::S end -function build_eigenvalue_solution(prob, alg, values, vectors; - retcode = ReturnCode.Success, resid = nothing, stats = nothing) +function build_eigenvalue_solution( + prob, alg, values, vectors; + retcode = ReturnCode.Success, resid = nothing, stats = nothing + ) T = eltype(eltype(values)) N = length((size(values)...,)) return EigenvalueSolution{ @@ -96,7 +100,8 @@ function SciMLBase.solve(prob::EigenvalueProblem, alg::DenseEigen, args...; kwar end function SciMLBase.solve( - prob::EigenvalueProblem, alg::AbstractEigenvalueAlgorithm, args...; kwargs...) + prob::EigenvalueProblem, alg::AbstractEigenvalueAlgorithm, args...; kwargs... + ) error("The eigenvalue backend $(typeof(alg)) is not available. Load its package before solving with this algorithm.") end diff --git a/test/Core/eigenvalue.jl b/test/Core/eigenvalue.jl index ac5a6856f..6739d5c82 100644 --- a/test/Core/eigenvalue.jl +++ b/test/Core/eigenvalue.jl @@ -38,7 +38,8 @@ sol_arpack = solve(EigenvalueProblem(Matrix(A_backend); nev = 2), ArpackJL()) @test sol_arpack.u ≈ [8.0, 7.0] sol_arpack_gen = solve( - EigenvalueProblem(Matrix(A_backend), Matrix(B_backend); nev = 2), ArpackJL()) + EigenvalueProblem(Matrix(A_backend), Matrix(B_backend); nev = 2), ArpackJL() +) @test sol_arpack_gen.u ≈ [4.0, 3.5] sol_arnoldi = solve(EigenvalueProblem(A_backend; nev = 2), LinearSolve.ArnoldiMethod()) @@ -48,19 +49,22 @@ sol_arnoldi_default = solve(EigenvalueProblem(A_backend), LinearSolve.ArnoldiMet @test length(sol_arnoldi_default.u) == 6 sol_arnoldi_shift = solve( - EigenvalueProblem(A_backend; nev = 2, sigma = 3.2), LinearSolve.ArnoldiMethod()) + EigenvalueProblem(A_backend; nev = 2, sigma = 3.2), LinearSolve.ArnoldiMethod() +) @test sol_arnoldi_shift.u ≈ [3.0, 4.0] sol_krylovkit = solve(EigenvalueProblem(A_backend; nev = 2), KrylovKitEigen()) @test sol_krylovkit.u ≈ [8.0, 7.0] sol_krylovkit_shift = solve( - EigenvalueProblem(Matrix(A_backend); nev = 2, sigma = 3.2), KrylovKitEigen()) + EigenvalueProblem(Matrix(A_backend); nev = 2, sigma = 3.2), KrylovKitEigen() +) @test sol_krylovkit_shift.u ≈ [3.0, 4.0] sol_krylovkit_gen_shift = solve( EigenvalueProblem(Matrix(A_backend), Matrix(B_backend); nev = 2, sigma = 1.6), - KrylovKitEigen()) + KrylovKitEigen() +) @test sol_krylovkit_gen_shift.u ≈ [1.5, 2.0] # Jacobi-Davidson is a target/interior method: it finds the eigenvalues nearest @@ -71,7 +75,7 @@ A_jd = Matrix(Diagonal(Float64.(1:30))) sol_jd = solve(EigenvalueProblem(A_jd; nev = 2, sigma = 10.3), JacobiDavidsonJL()) @test sort(real(sol_jd.u)) ≈ [10.0, 11.0] @test sol_jd.retcode === ReturnCode.Success -@test norm(A_jd * sol_jd.vectors - sol_jd.vectors * Diagonal(sol_jd.u)) < 1e-6 +@test norm(A_jd * sol_jd.vectors - sol_jd.vectors * Diagonal(sol_jd.u)) < 1.0e-6 sol_jd_sm = solve(EigenvalueProblem(A_jd; nev = 1, which = :SM), JacobiDavidsonJL()) @test real(sol_jd_sm.u[1]) ≈ 1.0 @@ -79,4 +83,5 @@ sol_jd_sm = solve(EigenvalueProblem(A_jd; nev = 1, which = :SM), JacobiDavidsonJ # Generalized problems are not supported by the JacobiDavidson backend (upstream # `jdqz` is broken); it should error rather than silently misbehave. @test_throws ErrorException solve( - EigenvalueProblem(A_jd, Matrix(Diagonal(fill(2.0, 30))); nev = 1), JacobiDavidsonJL()) + EigenvalueProblem(A_jd, Matrix(Diagonal(fill(2.0, 30))); nev = 1), JacobiDavidsonJL() +) From 7ff3af4c8a197aafac30b987793c03dd04e4eb98 Mon Sep 17 00:00:00 2001 From: Jash Date: Fri, 3 Jul 2026 21:39:25 +0530 Subject: [PATCH 3/7] fixed runic failure --- test/Core/forwarddiff_overloads.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/Core/forwarddiff_overloads.jl b/test/Core/forwarddiff_overloads.jl index 5578b1132..2f68e04f3 100644 --- a/test/Core/forwarddiff_overloads.jl +++ b/test/Core/forwarddiff_overloads.jl @@ -293,8 +293,8 @@ prob = LinearProblem(sparse(plain_A), b) for nchunk in (1, 2, 3) bd = [ ForwardDiff.Dual{Nothing, Float64, nchunk}( - Float64(i), ForwardDiff.Partials(ntuple(k -> sin(i + k), nchunk)) - ) for i in 1:5 + Float64(i), ForwardDiff.Partials(ntuple(k -> sin(i + k), nchunk)) + ) for i in 1:5 ] cache = LinearSolve.__init(LinearProblem(Asp, bd), PureKLUFactorization()) @test eltype(cache.A) == Float64 # A not promoted @@ -303,9 +303,9 @@ prob = LinearProblem(sparse(plain_A), b) @test isapprox(ForwardDiff.value.(u), ForwardDiff.value.(uref); rtol = 1.0e-10) @test all( isapprox( - ForwardDiff.partials(u[i], j), ForwardDiff.partials(uref[i], j); - rtol = 1.0e-8, atol = 1.0e-12 - ) for i in 1:5, j in 1:nchunk + ForwardDiff.partials(u[i], j), ForwardDiff.partials(uref[i], j); + rtol = 1.0e-8, atol = 1.0e-12 + ) for i in 1:5, j in 1:nchunk ) end From 67222619b68b1e52f59fc1437fc1cb7484c72cf3 Mon Sep 17 00:00:00 2001 From: Jash Date: Tue, 7 Jul 2026 10:42:04 +0530 Subject: [PATCH 4/7] addressed reviews --- ext/LinearSolveArnoldiMethodExt.jl | 4 +- ext/LinearSolveArpackExt.jl | 6 +- ext/LinearSolveJacobiDavidsonExt.jl | 16 +-- ext/LinearSolveKrylovKitExt.jl | 16 +-- src/LinearSolve.jl | 6 +- src/eigenvalue.jl | 216 +++++++++++++++++++--------- test/Core/eigenvalue.jl | 15 ++ 7 files changed, 182 insertions(+), 97 deletions(-) diff --git a/ext/LinearSolveArnoldiMethodExt.jl b/ext/LinearSolveArnoldiMethodExt.jl index f222742a5..eb0552337 100644 --- a/ext/LinearSolveArnoldiMethodExt.jl +++ b/ext/LinearSolveArnoldiMethodExt.jl @@ -13,10 +13,10 @@ function SciMLBase.solve( prob.B === nothing || error("ArnoldiMethod backend currently supports standard eigenvalue problems only.") nev = LinearSolve.default_nev(prob) - which = prob.sigma === nothing ? prob.which : :LM + which = prob.sigma === nothing ? LinearSolve._target_symbol(prob.which) : :LM A = prob.sigma === nothing ? prob.A : _shift_invert_operator(prob.A, prob.sigma) kw = (; nev, which, prob.kwargs..., alg.kwargs..., kwargs...) - decomp, history = partialschur(A, alg.args...; kw...) + decomp, history = partialschur(A; kw...) values, vectors = partialeigen(decomp) if prob.sigma !== nothing values = prob.sigma .+ inv.(values) diff --git a/ext/LinearSolveArpackExt.jl b/ext/LinearSolveArpackExt.jl index 4acfc1e42..32eac1891 100644 --- a/ext/LinearSolveArpackExt.jl +++ b/ext/LinearSolveArpackExt.jl @@ -10,7 +10,7 @@ function SciMLBase.solve( args...; kwargs... ) nev = LinearSolve.default_nev(prob) - base = (; nev, which = prob.which) + base = (; nev, which = LinearSolve._target_symbol(prob.which)) if prob.sigma !== nothing base = (; base..., sigma = prob.sigma) end @@ -18,9 +18,9 @@ function SciMLBase.solve( # `Arpack.eigs` takes the generalized-problem matrix `B` positionally, not # as a keyword argument. values, vectors, nconv, niter, nmult, resid = if prob.B === nothing - Arpack.eigs(prob.A, alg.args...; kw...) + Arpack.eigs(prob.A; kw...) else - Arpack.eigs(prob.A, prob.B, alg.args...; kw...) + Arpack.eigs(prob.A, prob.B; kw...) end retcode = nconv >= length(values) ? ReturnCode.Success : ReturnCode.ConvergenceFailure stats = (; nconv, niter, nmult) diff --git a/ext/LinearSolveJacobiDavidsonExt.jl b/ext/LinearSolveJacobiDavidsonExt.jl index 9ae23f8af..3ef4e87d7 100644 --- a/ext/LinearSolveJacobiDavidsonExt.jl +++ b/ext/LinearSolveJacobiDavidsonExt.jl @@ -27,7 +27,7 @@ function SciMLBase.solve( defaults = (; pairs = nev, target = target, subspace_dimensions = lo:hi) kw = (; defaults..., prob.kwargs..., alg.kwargs..., kwargs...) - out = JacobiDavidson.jdqr(prob.A, alg.args...; kw...) + out = JacobiDavidson.jdqr(prob.A; kw...) values, vectors = _jd_standard_pairs(prob.A, out[1]) values, vectors = LinearSolve._select_eigenpairs( @@ -44,14 +44,14 @@ end # strength; otherwise the `which` symbol selects an extremal target. function _jd_target(prob) prob.sigma !== nothing && return JacobiDavidson.Near(ComplexF64(prob.sigma)) + T = LinearSolve.EigenvalueTarget w = prob.which - return w === :LM ? JacobiDavidson.LargestMagnitude(0.0 + 0.0im) : - w === :SM ? JacobiDavidson.SmallestMagnitude(0.0 + 0.0im) : - w === :LR ? JacobiDavidson.LargestRealPart(0.0 + 0.0im) : - w === :SR ? JacobiDavidson.SmallestRealPart(0.0 + 0.0im) : - w === :LI ? JacobiDavidson.LargestImaginaryPart(0.0 + 0.0im) : - w === :SI ? JacobiDavidson.SmallestImaginaryPart(0.0 + 0.0im) : - throw(ArgumentError("unsupported `which = $(w)` for JacobiDavidson; expected one of :LM, :SM, :LR, :SR, :LI, :SI, or pass `sigma`")) + return w == T.LargestMagnitude ? JacobiDavidson.LargestMagnitude(0.0 + 0.0im) : + w == T.SmallestMagnitude ? JacobiDavidson.SmallestMagnitude(0.0 + 0.0im) : + w == T.LargestRealPart ? JacobiDavidson.LargestRealPart(0.0 + 0.0im) : + w == T.SmallestRealPart ? JacobiDavidson.SmallestRealPart(0.0 + 0.0im) : + w == T.LargestImaginaryPart ? JacobiDavidson.LargestImaginaryPart(0.0 + 0.0im) : + JacobiDavidson.SmallestImaginaryPart(0.0 + 0.0im) end # jdqr yields a partial Schur decomposition `A*Q = Q*R`. Eigenpairs are recovered diff --git a/ext/LinearSolveKrylovKitExt.jl b/ext/LinearSolveKrylovKitExt.jl index b81372079..c6e86f3bb 100644 --- a/ext/LinearSolveKrylovKitExt.jl +++ b/ext/LinearSolveKrylovKitExt.jl @@ -67,14 +67,14 @@ function SciMLBase.solve( args...; kwargs... ) nev = LinearSolve.default_nev(prob) - which = prob.which + which = LinearSolve._target_symbol(prob.which) kw = (; prob.kwargs..., alg.kwargs..., kwargs...) values, vectors, info = if prob.sigma !== nothing - _shift_invert_eigsolve(prob, alg, nev, kw) + _shift_invert_eigsolve(prob, nev, kw) elseif prob.B === nothing - KrylovKit.eigsolve(prob.A, alg.args..., nev, which; kw...) + KrylovKit.eigsolve(prob.A, nev, which; kw...) else - KrylovKit.geneigsolve((prob.A, prob.B), alg.args..., nev, which; kw...) + KrylovKit.geneigsolve((prob.A, prob.B), nev, which; kw...) end if prob.sigma !== nothing values = prob.sigma .+ inv.(values) @@ -87,7 +87,7 @@ function SciMLBase.solve( ) end -function _shift_invert_eigsolve(prob, alg, nev, kw) +function _shift_invert_eigsolve(prob, nev, kw) A, B, sigma = prob.A, prob.B, prob.sigma T = isnothing(B) ? promote_type(eltype(A), typeof(sigma)) : promote_type(eltype(A), eltype(B), typeof(sigma)) @@ -101,11 +101,7 @@ function _shift_invert_eigsolve(prob, alg, nev, kw) F = factorize(A - sigma * B) op = x -> F \ (B * x) end - if isempty(alg.args) - return KrylovKit.eigsolve(op, size(A, 2), nev, :LM, T; kw...) - else - return KrylovKit.eigsolve(op, alg.args..., nev, :LM; kw...) - end + return KrylovKit.eigsolve(op, size(A, 2), nev, :LM, T; kw...) end end diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index dd02ed9ac..ebe192bca 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -574,8 +574,8 @@ export LinearSolveAdjoint export LinearVerbosity -export EigenvalueProblem, EigenvalueSolution, AbstractEigenvalueAlgorithm, - DenseEigen, ArpackJL, ArnoldiMethod, ArnoldiMethodJL, KrylovKitEigen, - JacobiDavidsonJL +export AbstractEigenvalueAlgorithm, + DenseEigen, ArpackJL, ArnoldiMethod, ArnoldiMethodJL, + KrylovKitEigen, JacobiDavidsonJL end diff --git a/src/eigenvalue.jl b/src/eigenvalue.jl index 23843cbc9..d20e34b59 100644 --- a/src/eigenvalue.jl +++ b/src/eigenvalue.jl @@ -1,83 +1,149 @@ -""" - EigenvalueProblem(A[, B]; nev = nothing, which = :LM, sigma = nothing, u0 = nothing, p = SciMLBase.NullParameters()) - -Define a standard or generalized eigenvalue problem. - -The standard problem is ``A v = lambda v``. If `B` is supplied, the generalized -problem is ``A v = lambda B v``. `nev` requests a subset of eigenpairs, `which` -selects the part of the spectrum (`:LM`, `:SM`, `:LR`, `:SR`, `:LI`, `:SI`), -and `sigma` requests eigenvalues nearest to the shift. -""" -struct EigenvalueProblem{AType, BType, NevType, WhichType, SigmaType, U0Type, PType, KType} - A::AType - B::BType - nev::NevType - which::WhichType - sigma::SigmaType - u0::U0Type - p::PType - kwargs::KType -end +# `EigenvalueProblem`, `EigenvalueSolution`, `EigenvalueTarget`, and +# `build_eigenvalue_solution` are defined natively in SciMLBase (analogous to +# `LinearProblem`/`LinearSolution`) once an upstream release adds them. Older +# SciMLBase versions lack them, so fall back to local definitions here, +# preserving the same public interface either way. This mirrors the existing +# `@static if isdefined(SciMLBase, :DiffEqArrayOperator)` gate further down in +# this module for the same reason: keep a wide SciMLBase compat range instead +# of forcing every LinearSolve user onto an upstream version they may not need. +@static if isdefined(SciMLBase, :EigenvalueProblem) + using SciMLBase: EigenvalueProblem, EigenvalueSolution, EigenvalueTarget, + build_eigenvalue_solution +else + """ + EigenvalueTarget + + Enum selecting which part of the spectrum is returned when only a subset of the + eigenpairs is requested (via `nev`) in an [`EigenvalueProblem`](@ref). The `which` + keyword accepts either an `EigenvalueTarget` value or, for convenience, the + corresponding ARPACK-style `Symbol` noted for each variant below. + """ + EnumX.@enumx EigenvalueTarget begin + "Eigenvalues of largest magnitude, `abs(λ)` largest (symbol `:LM`)." + LargestMagnitude + "Eigenvalues of smallest magnitude, `abs(λ)` smallest (symbol `:SM`)." + SmallestMagnitude + "Eigenvalues with the largest (most positive) real part (symbol `:LR`)." + LargestRealPart + "Eigenvalues with the smallest (most negative) real part (symbol `:SR`)." + SmallestRealPart + "Eigenvalues with the largest (most positive) imaginary part (symbol `:LI`)." + LargestImaginaryPart + "Eigenvalues with the smallest (most negative) imaginary part (symbol `:SI`)." + SmallestImaginaryPart + end -function EigenvalueProblem( - A, B = nothing, p = SciMLBase.NullParameters(); - nev = nothing, which = :LM, sigma = nothing, u0 = nothing, kwargs... - ) - return EigenvalueProblem{ - typeof(A), typeof(B), typeof(nev), typeof(which), typeof(sigma), - typeof(u0), typeof(p), typeof(kwargs), - }(A, B, nev, which, sigma, u0, p, kwargs) + # Normalize a user-supplied `which` (an `EigenvalueTarget` or an ARPACK-style + # `Symbol`) to an `EigenvalueTarget`, throwing on anything else. + _eigenvalue_target(w::EigenvalueTarget.T) = w + function _eigenvalue_target(w::Symbol) + return w === :LM ? EigenvalueTarget.LargestMagnitude : + w === :SM ? EigenvalueTarget.SmallestMagnitude : + w === :LR ? EigenvalueTarget.LargestRealPart : + w === :SR ? EigenvalueTarget.SmallestRealPart : + w === :LI ? EigenvalueTarget.LargestImaginaryPart : + w === :SI ? EigenvalueTarget.SmallestImaginaryPart : + throw(ArgumentError("unsupported eigenvalue selector `which = $(repr(w))`; expected an `EigenvalueTarget` or one of :LM, :SM, :LR, :SR, :LI, :SI")) + end + _eigenvalue_target(w) = throw(ArgumentError("`which` must be an `EigenvalueTarget` or a Symbol (:LM, :SM, :LR, :SR, :LI, :SI), got $(repr(w))")) + + """ + EigenvalueProblem(A[, B], p = SciMLBase.NullParameters(); + nev = nothing, which = EigenvalueTarget.LargestMagnitude, + sigma = nothing, u0 = nothing) + + Define a standard or generalized eigenvalue problem. + + The standard problem is ``A v = λ v``. If `B` is supplied, the generalized problem + is ``A v = λ B v``. + + ## Keyword arguments + + - `nev`: the number of eigenpairs (eigenvalues together with their eigenvectors) to + compute. `nothing` (the default) requests every eigenpair for the dense solver, or + a solver-chosen default for the iterative backends. + - `which`: which part of the spectrum to return, as an [`EigenvalueTarget`](@ref). + An ARPACK-style `Symbol` (`:LM`, `:SM`, `:LR`, `:SR`, `:LI`, `:SI`) is also accepted + and converted to the corresponding `EigenvalueTarget`. Defaults to the eigenvalues of + largest magnitude. + - `sigma`: if supplied, return the eigenvalues nearest this shift (shift-and-invert). + - `u0`: optional initial guess for the iterative backends. + """ + struct EigenvalueProblem{AType, BType, NevType, WhichType, SigmaType, U0Type, PType, KType} + A::AType + B::BType + nev::NevType + which::WhichType + sigma::SigmaType + u0::U0Type + p::PType + kwargs::KType + end + + function EigenvalueProblem( + A, B = nothing, p = SciMLBase.NullParameters(); + nev = nothing, which = EigenvalueTarget.LargestMagnitude, + sigma = nothing, u0 = nothing, kwargs... + ) + target = _eigenvalue_target(which) + return EigenvalueProblem{ + typeof(A), typeof(B), typeof(nev), typeof(target), typeof(sigma), + typeof(u0), typeof(p), typeof(kwargs), + }(A, B, nev, target, sigma, u0, p, kwargs) + end + + struct EigenvalueSolution{T, N, U, V, P, A, R, S} <: SciMLBase.AbstractNoTimeSolution{T, N} + u::U + vectors::V + prob::P + alg::A + retcode::ReturnCode.T + resid::R + stats::S + end + + function build_eigenvalue_solution( + prob, alg, values, vectors; + retcode = ReturnCode.Success, resid = nothing, stats = nothing + ) + T = eltype(eltype(values)) + N = length((size(values)...,)) + return EigenvalueSolution{ + T, N, typeof(values), typeof(vectors), typeof(prob), typeof(alg), + typeof(resid), typeof(stats), + }(values, vectors, prob, alg, retcode, resid, stats) + end + + export EigenvalueProblem, EigenvalueSolution, EigenvalueTarget end abstract type AbstractEigenvalueAlgorithm <: SciMLBase.AbstractLinearAlgorithm end struct DenseEigen <: AbstractEigenvalueAlgorithm end -struct ArpackJL{A, K} <: AbstractEigenvalueAlgorithm - args::A +# The iterative backends forward any extra keyword arguments to the underlying +# solver (`Arpack.eigs`, `ArnoldiMethod.partialschur`, `KrylovKit.eigsolve`, +# `JacobiDavidson.jdqr`). They are keyword-only: passing positional arguments is +# an error, and unrecognized keywords are rejected by the underlying solver. +struct ArpackJL{K <: NamedTuple} <: AbstractEigenvalueAlgorithm kwargs::K end -ArpackJL(args...; kwargs...) = ArpackJL(args, kwargs) +ArpackJL(; kwargs...) = ArpackJL((; kwargs...)) -struct ArnoldiMethodJL{A, K} <: AbstractEigenvalueAlgorithm - args::A +struct ArnoldiMethodJL{K <: NamedTuple} <: AbstractEigenvalueAlgorithm kwargs::K end -ArnoldiMethod(args...; kwargs...) = ArnoldiMethodJL(args, kwargs) +ArnoldiMethod(; kwargs...) = ArnoldiMethodJL((; kwargs...)) -struct KrylovKitEigen{A, K} <: AbstractEigenvalueAlgorithm - args::A +struct KrylovKitEigen{K <: NamedTuple} <: AbstractEigenvalueAlgorithm kwargs::K end -KrylovKitEigen(args...; kwargs...) = KrylovKitEigen(args, kwargs) +KrylovKitEigen(; kwargs...) = KrylovKitEigen((; kwargs...)) -struct JacobiDavidsonJL{A, K} <: AbstractEigenvalueAlgorithm - args::A +struct JacobiDavidsonJL{K <: NamedTuple} <: AbstractEigenvalueAlgorithm kwargs::K end -JacobiDavidsonJL(args...; kwargs...) = JacobiDavidsonJL(args, kwargs) - -struct EigenvalueSolution{T, N, U, V, P, A, R, S} <: SciMLBase.AbstractNoTimeSolution{T, N} - u::U - vectors::V - prob::P - alg::A - retcode::ReturnCode.T - resid::R - stats::S -end - -function build_eigenvalue_solution( - prob, alg, values, vectors; - retcode = ReturnCode.Success, resid = nothing, stats = nothing - ) - T = eltype(eltype(values)) - N = length((size(values)...,)) - return EigenvalueSolution{ - T, N, typeof(values), typeof(vectors), typeof(prob), typeof(alg), - typeof(resid), typeof(stats), - }(values, vectors, prob, alg, retcode, resid, stats) -end +JacobiDavidsonJL(; kwargs...) = JacobiDavidsonJL((; kwargs...)) SciMLBase.solve(prob::EigenvalueProblem, args...; kwargs...) = solve(prob, nothing, args...; kwargs...) @@ -120,22 +186,30 @@ function _select_eigenpairs(values, vectors, nev, which, sigma) return values[idxs], vectors[:, idxs] end -function _eigenvalue_order(values, which, sigma) +# ARPACK-style symbol for the backends that take `which` as a `Symbol`. +function _target_symbol(w::EigenvalueTarget.T) + return w == EigenvalueTarget.LargestMagnitude ? :LM : + w == EigenvalueTarget.SmallestMagnitude ? :SM : + w == EigenvalueTarget.LargestRealPart ? :LR : + w == EigenvalueTarget.SmallestRealPart ? :SR : + w == EigenvalueTarget.LargestImaginaryPart ? :LI : + :SI +end + +function _eigenvalue_order(values, which::EigenvalueTarget.T, sigma) if sigma !== nothing return sortperm(abs.(values .- sigma)) - elseif which === :LM + elseif which == EigenvalueTarget.LargestMagnitude return sortperm(abs.(values); rev = true) - elseif which === :SM + elseif which == EigenvalueTarget.SmallestMagnitude return sortperm(abs.(values)) - elseif which === :LR + elseif which == EigenvalueTarget.LargestRealPart return sortperm(real.(values); rev = true) - elseif which === :SR + elseif which == EigenvalueTarget.SmallestRealPart return sortperm(real.(values)) - elseif which === :LI + elseif which == EigenvalueTarget.LargestImaginaryPart return sortperm(imag.(values); rev = true) - elseif which === :SI + else # EigenvalueTarget.SmallestImaginaryPart return sortperm(imag.(values)) - else - throw(ArgumentError("unsupported eigenvalue selector `which = $which`; expected one of :LM, :SM, :LR, :SR, :LI, or :SI")) end end diff --git a/test/Core/eigenvalue.jl b/test/Core/eigenvalue.jl index 6739d5c82..86a82c210 100644 --- a/test/Core/eigenvalue.jl +++ b/test/Core/eigenvalue.jl @@ -85,3 +85,18 @@ sol_jd_sm = solve(EigenvalueProblem(A_jd; nev = 1, which = :SM), JacobiDavidsonJ @test_throws ErrorException solve( EigenvalueProblem(A_jd, Matrix(Diagonal(fill(2.0, 30))); nev = 1), JacobiDavidsonJL() ) + +# `which` accepts the EigenvalueTarget enum directly, equivalent to its Symbol alias. +prob_enum = EigenvalueProblem(A; nev = 2, which = EigenvalueTarget.SmallestMagnitude) +@test prob_enum.which === EigenvalueTarget.SmallestMagnitude +@test solve(prob_enum).u == [1.0, 2.0] +@test EigenvalueProblem(A; which = :SM).which === EigenvalueTarget.SmallestMagnitude + +# An invalid `which` is rejected at construction. +@test_throws ArgumentError EigenvalueProblem(A; which = :XX) +@test_throws ArgumentError EigenvalueProblem(A; which = 5) + +# The iterative algorithms are keyword-only: positional arguments error, and +# forwarded keywords reach the backend. +@test_throws MethodError ArpackJL(5) +@test ArpackJL(; maxiter = 500).kwargs.maxiter == 500 From 8020b29559b820aec709c1749f207abba4f5bb36 Mon Sep 17 00:00:00 2001 From: Jash Date: Tue, 7 Jul 2026 15:33:23 +0530 Subject: [PATCH 5/7] addressed reviews --- ext/LinearSolveArnoldiMethodExt.jl | 35 ++++++--- ext/LinearSolveArpackExt.jl | 21 +++++- ext/LinearSolveJacobiDavidsonExt.jl | 12 +-- ext/LinearSolveKrylovKitExt.jl | 38 +++++++--- src/eigenvalue.jl | 110 +++++++++++----------------- test/Core/eigenvalue.jl | 68 ++++++++++------- 6 files changed, 160 insertions(+), 124 deletions(-) diff --git a/ext/LinearSolveArnoldiMethodExt.jl b/ext/LinearSolveArnoldiMethodExt.jl index eb0552337..596c1c0d7 100644 --- a/ext/LinearSolveArnoldiMethodExt.jl +++ b/ext/LinearSolveArnoldiMethodExt.jl @@ -2,7 +2,7 @@ module LinearSolveArnoldiMethodExt using LinearAlgebra using LinearSolve -import ArnoldiMethod: partialschur, partialeigen +import ArnoldiMethod: partialschur, partialeigen, LM, LR, SR, LI, SI using SciMLBase: SciMLBase, ReturnCode function SciMLBase.solve( @@ -12,25 +12,40 @@ function SciMLBase.solve( ) prob.B === nothing || error("ArnoldiMethod backend currently supports standard eigenvalue problems only.") - nev = LinearSolve.default_nev(prob) - which = prob.sigma === nothing ? LinearSolve._target_symbol(prob.which) : :LM - A = prob.sigma === nothing ? prob.A : _shift_invert_operator(prob.A, prob.sigma) + nev = LinearSolve.default_num_eigenpairs(prob) + which = prob.shift === nothing ? _arnoldi_target(prob.eigentarget) : LM() + A = prob.shift === nothing ? prob.A : _shift_invert_operator(prob.A, prob.shift) kw = (; nev, which, prob.kwargs..., alg.kwargs..., kwargs...) decomp, history = partialschur(A; kw...) values, vectors = partialeigen(decomp) - if prob.sigma !== nothing - values = prob.sigma .+ inv.(values) + if prob.shift !== nothing + values = prob.shift .+ inv.(values) end - values, vectors = LinearSolve._select_eigenpairs(values, vectors, nev, prob.which, prob.sigma) + values, vectors = LinearSolve._select_eigenpairs( + values, vectors, nev, prob.eigentarget, prob.shift + ) retcode = history.converged ? ReturnCode.Success : ReturnCode.ConvergenceFailure return LinearSolve.build_eigenvalue_solution( prob, alg, values, vectors; retcode, stats = history ) end -function _shift_invert_operator(A, sigma) - F = factorize(A - sigma * I) - T = promote_type(eltype(A), typeof(sigma)) +# ArnoldiMethod exposes its own `Target` types (preferred over ARPACK-style +# symbols per its own documentation) for all but smallest-magnitude, which it +# does not support at all. +function _arnoldi_target(w::LinearSolve.EigenvalueTarget.T) + T = LinearSolve.EigenvalueTarget + return w == T.LargestMagnitude ? LM() : + w == T.LargestRealPart ? LR() : + w == T.SmallestRealPart ? SR() : + w == T.LargestImaginaryPart ? LI() : + w == T.SmallestImaginaryPart ? SI() : + throw(ArgumentError("ArnoldiMethod does not support `eigentarget = EigenvalueTarget.SmallestMagnitude`; use a different backend (e.g. `KrylovKitEigen()` or `ArpackJL()`) or supply `shift` for shift-and-invert.")) +end + +function _shift_invert_operator(A, shift) + F = factorize(A - shift * I) + T = promote_type(eltype(A), typeof(shift)) return ShiftInvertMap{typeof(F), T}(F, size(A, 1)) end diff --git a/ext/LinearSolveArpackExt.jl b/ext/LinearSolveArpackExt.jl index 32eac1891..9d02caa72 100644 --- a/ext/LinearSolveArpackExt.jl +++ b/ext/LinearSolveArpackExt.jl @@ -9,10 +9,10 @@ function SciMLBase.solve( alg::LinearSolve.ArpackJL, args...; kwargs... ) - nev = LinearSolve.default_nev(prob) - base = (; nev, which = LinearSolve._target_symbol(prob.which)) - if prob.sigma !== nothing - base = (; base..., sigma = prob.sigma) + nev = LinearSolve.default_num_eigenpairs(prob) + base = (; nev, which = _arpack_which(prob.eigentarget)) + if prob.shift !== nothing + base = (; base..., sigma = prob.shift) end kw = (; base..., prob.kwargs..., alg.kwargs..., kwargs...) # `Arpack.eigs` takes the generalized-problem matrix `B` positionally, not @@ -29,4 +29,17 @@ function SciMLBase.solve( ) end +# Arpack.eigs requires a raw ARPACK-style Symbol for `which`; this mapping is +# purely a private adapter to that third-party API, not a general LinearSolve +# concept. +function _arpack_which(w::LinearSolve.EigenvalueTarget.T) + T = LinearSolve.EigenvalueTarget + return w == T.LargestMagnitude ? :LM : + w == T.SmallestMagnitude ? :SM : + w == T.LargestRealPart ? :LR : + w == T.SmallestRealPart ? :SR : + w == T.LargestImaginaryPart ? :LI : + :SI +end + end diff --git a/ext/LinearSolveJacobiDavidsonExt.jl b/ext/LinearSolveJacobiDavidsonExt.jl index 3ef4e87d7..7d4f617a9 100644 --- a/ext/LinearSolveJacobiDavidsonExt.jl +++ b/ext/LinearSolveJacobiDavidsonExt.jl @@ -18,7 +18,7 @@ function SciMLBase.solve( error("The JacobiDavidson backend currently supports standard eigenvalue problems only. Use `ArpackJL()` or `KrylovKitEigen()` for generalized problems.") n = size(prob.A, 2) - nev = LinearSolve.default_nev(prob) + nev = LinearSolve.default_num_eigenpairs(prob) target = _jd_target(prob) # Search-subspace bounds, capped at the problem size. Users may override # `subspace_dimensions` (and any other jdqr keyword) via the algorithm. @@ -31,7 +31,7 @@ function SciMLBase.solve( values, vectors = _jd_standard_pairs(prob.A, out[1]) values, vectors = LinearSolve._select_eigenpairs( - values, vectors, nev, prob.which, prob.sigma + values, vectors, nev, prob.eigentarget, prob.shift ) retcode = length(values) >= nev ? ReturnCode.Success : ReturnCode.ConvergenceFailure return LinearSolve.build_eigenvalue_solution( @@ -40,12 +40,12 @@ function SciMLBase.solve( end # Map the problem's spectral selector onto a JacobiDavidson `Target`. A supplied -# `sigma` is the natural interior target (`Near`), which is Jacobi-Davidson's -# strength; otherwise the `which` symbol selects an extremal target. +# `shift` is the natural interior target (`Near`), which is Jacobi-Davidson's +# strength; otherwise `eigentarget` selects an extremal target. function _jd_target(prob) - prob.sigma !== nothing && return JacobiDavidson.Near(ComplexF64(prob.sigma)) + prob.shift !== nothing && return JacobiDavidson.Near(ComplexF64(prob.shift)) T = LinearSolve.EigenvalueTarget - w = prob.which + w = prob.eigentarget return w == T.LargestMagnitude ? JacobiDavidson.LargestMagnitude(0.0 + 0.0im) : w == T.SmallestMagnitude ? JacobiDavidson.SmallestMagnitude(0.0 + 0.0im) : w == T.LargestRealPart ? JacobiDavidson.LargestRealPart(0.0 + 0.0im) : diff --git a/ext/LinearSolveKrylovKitExt.jl b/ext/LinearSolveKrylovKitExt.jl index c6e86f3bb..0ad0d3b23 100644 --- a/ext/LinearSolveKrylovKitExt.jl +++ b/ext/LinearSolveKrylovKitExt.jl @@ -66,39 +66,53 @@ function SciMLBase.solve( alg::LinearSolve.KrylovKitEigen, args...; kwargs... ) - nev = LinearSolve.default_nev(prob) - which = LinearSolve._target_symbol(prob.which) + nev = LinearSolve.default_num_eigenpairs(prob) + which = _krylovkit_which(prob.eigentarget) kw = (; prob.kwargs..., alg.kwargs..., kwargs...) - values, vectors, info = if prob.sigma !== nothing + values, vectors, info = if prob.shift !== nothing _shift_invert_eigsolve(prob, nev, kw) elseif prob.B === nothing KrylovKit.eigsolve(prob.A, nev, which; kw...) else KrylovKit.geneigsolve((prob.A, prob.B), nev, which; kw...) end - if prob.sigma !== nothing - values = prob.sigma .+ inv.(values) + if prob.shift !== nothing + values = prob.shift .+ inv.(values) end vecmat = reduce(hcat, vectors) - values, vecmat = LinearSolve._select_eigenpairs(values, vecmat, nev, prob.which, prob.sigma) + values, vecmat = LinearSolve._select_eigenpairs( + values, vecmat, nev, prob.eigentarget, prob.shift + ) retcode = info.converged >= length(values) ? ReturnCode.Success : ReturnCode.ConvergenceFailure return LinearSolve.build_eigenvalue_solution( prob, alg, values, vecmat; retcode, resid = info.normres, stats = info ) end +# KrylovKit.eigsolve requires a raw ARPACK-style Symbol for `which`; this +# mapping is purely a private adapter to that third-party API. +function _krylovkit_which(w::LinearSolve.EigenvalueTarget.T) + T = LinearSolve.EigenvalueTarget + return w == T.LargestMagnitude ? :LM : + w == T.SmallestMagnitude ? :SM : + w == T.LargestRealPart ? :LR : + w == T.SmallestRealPart ? :SR : + w == T.LargestImaginaryPart ? :LI : + :SI +end + function _shift_invert_eigsolve(prob, nev, kw) - A, B, sigma = prob.A, prob.B, prob.sigma - T = isnothing(B) ? promote_type(eltype(A), typeof(sigma)) : - promote_type(eltype(A), eltype(B), typeof(sigma)) + A, B, shift = prob.A, prob.B, prob.shift + T = isnothing(B) ? promote_type(eltype(A), typeof(shift)) : + promote_type(eltype(A), eltype(B), typeof(shift)) if isnothing(B) - F = factorize(A - sigma * I) + F = factorize(A - shift * I) op = x -> F \ x elseif B isa LinearAlgebra.UniformScaling - F = factorize(A - sigma * B) + F = factorize(A - shift * B) op = x -> F \ (B.λ * x) else - F = factorize(A - sigma * B) + F = factorize(A - shift * B) op = x -> F \ (B * x) end return KrylovKit.eigsolve(op, size(A, 2), nev, :LM, T; kw...) diff --git a/src/eigenvalue.jl b/src/eigenvalue.jl index d20e34b59..993b2461d 100644 --- a/src/eigenvalue.jl +++ b/src/eigenvalue.jl @@ -14,43 +14,27 @@ else EigenvalueTarget Enum selecting which part of the spectrum is returned when only a subset of the - eigenpairs is requested (via `nev`) in an [`EigenvalueProblem`](@ref). The `which` - keyword accepts either an `EigenvalueTarget` value or, for convenience, the - corresponding ARPACK-style `Symbol` noted for each variant below. + eigenpairs is requested (via `num_eigenpairs`) in an [`EigenvalueProblem`](@ref). """ EnumX.@enumx EigenvalueTarget begin - "Eigenvalues of largest magnitude, `abs(λ)` largest (symbol `:LM`)." + "Eigenvalues of largest magnitude, `abs(λ)` largest." LargestMagnitude - "Eigenvalues of smallest magnitude, `abs(λ)` smallest (symbol `:SM`)." + "Eigenvalues of smallest magnitude, `abs(λ)` smallest." SmallestMagnitude - "Eigenvalues with the largest (most positive) real part (symbol `:LR`)." + "Eigenvalues with the largest (most positive) real part." LargestRealPart - "Eigenvalues with the smallest (most negative) real part (symbol `:SR`)." + "Eigenvalues with the smallest (most negative) real part." SmallestRealPart - "Eigenvalues with the largest (most positive) imaginary part (symbol `:LI`)." + "Eigenvalues with the largest (most positive) imaginary part." LargestImaginaryPart - "Eigenvalues with the smallest (most negative) imaginary part (symbol `:SI`)." + "Eigenvalues with the smallest (most negative) imaginary part." SmallestImaginaryPart end - # Normalize a user-supplied `which` (an `EigenvalueTarget` or an ARPACK-style - # `Symbol`) to an `EigenvalueTarget`, throwing on anything else. - _eigenvalue_target(w::EigenvalueTarget.T) = w - function _eigenvalue_target(w::Symbol) - return w === :LM ? EigenvalueTarget.LargestMagnitude : - w === :SM ? EigenvalueTarget.SmallestMagnitude : - w === :LR ? EigenvalueTarget.LargestRealPart : - w === :SR ? EigenvalueTarget.SmallestRealPart : - w === :LI ? EigenvalueTarget.LargestImaginaryPart : - w === :SI ? EigenvalueTarget.SmallestImaginaryPart : - throw(ArgumentError("unsupported eigenvalue selector `which = $(repr(w))`; expected an `EigenvalueTarget` or one of :LM, :SM, :LR, :SR, :LI, :SI")) - end - _eigenvalue_target(w) = throw(ArgumentError("`which` must be an `EigenvalueTarget` or a Symbol (:LM, :SM, :LR, :SR, :LI, :SI), got $(repr(w))")) - """ EigenvalueProblem(A[, B], p = SciMLBase.NullParameters(); - nev = nothing, which = EigenvalueTarget.LargestMagnitude, - sigma = nothing, u0 = nothing) + num_eigenpairs = nothing, eigentarget = EigenvalueTarget.LargestMagnitude, + shift = nothing, u0 = nothing) Define a standard or generalized eigenvalue problem. @@ -59,22 +43,22 @@ else ## Keyword arguments - - `nev`: the number of eigenpairs (eigenvalues together with their eigenvectors) to - compute. `nothing` (the default) requests every eigenpair for the dense solver, or - a solver-chosen default for the iterative backends. - - `which`: which part of the spectrum to return, as an [`EigenvalueTarget`](@ref). - An ARPACK-style `Symbol` (`:LM`, `:SM`, `:LR`, `:SR`, `:LI`, `:SI`) is also accepted - and converted to the corresponding `EigenvalueTarget`. Defaults to the eigenvalues of - largest magnitude. - - `sigma`: if supplied, return the eigenvalues nearest this shift (shift-and-invert). + - `num_eigenpairs`: the number of eigenpairs (eigenvalues together with their + eigenvectors) to compute. `nothing` (the default) requests every eigenpair for the + dense solver, or a solver-chosen default for the iterative backends. + - `eigentarget`: which part of the spectrum to return, as an + [`EigenvalueTarget`](@ref). Defaults to the eigenvalues of largest magnitude. + - `shift`: if supplied, return the eigenvalues nearest this shift (shift-and-invert). - `u0`: optional initial guess for the iterative backends. """ - struct EigenvalueProblem{AType, BType, NevType, WhichType, SigmaType, U0Type, PType, KType} + struct EigenvalueProblem{ + AType, BType, NevType, TargetType, ShiftType, U0Type, PType, KType, + } A::AType B::BType - nev::NevType - which::WhichType - sigma::SigmaType + num_eigenpairs::NevType + eigentarget::TargetType + shift::ShiftType u0::U0Type p::PType kwargs::KType @@ -82,14 +66,14 @@ else function EigenvalueProblem( A, B = nothing, p = SciMLBase.NullParameters(); - nev = nothing, which = EigenvalueTarget.LargestMagnitude, - sigma = nothing, u0 = nothing, kwargs... + num_eigenpairs = nothing, + eigentarget::EigenvalueTarget.T = EigenvalueTarget.LargestMagnitude, + shift = nothing, u0 = nothing, kwargs... ) - target = _eigenvalue_target(which) return EigenvalueProblem{ - typeof(A), typeof(B), typeof(nev), typeof(target), typeof(sigma), - typeof(u0), typeof(p), typeof(kwargs), - }(A, B, nev, target, sigma, u0, p, kwargs) + typeof(A), typeof(B), typeof(num_eigenpairs), typeof(eigentarget), + typeof(shift), typeof(u0), typeof(p), typeof(kwargs), + }(A, B, num_eigenpairs, eigentarget, shift, u0, p, kwargs) end struct EigenvalueSolution{T, N, U, V, P, A, R, S} <: SciMLBase.AbstractNoTimeSolution{T, N} @@ -161,7 +145,9 @@ function SciMLBase.solve(prob::EigenvalueProblem, alg::DenseEigen, args...; kwar else LinearAlgebra.eigen(prob.A, prob.B; kw...) end - values, vectors = _select_eigenpairs(F.values, F.vectors, prob.nev, prob.which, prob.sigma) + values, vectors = _select_eigenpairs( + F.values, F.vectors, prob.num_eigenpairs, prob.eigentarget, prob.shift + ) return build_eigenvalue_solution(prob, alg, values, vectors) end @@ -171,43 +157,33 @@ function SciMLBase.solve( error("The eigenvalue backend $(typeof(alg)) is not available. Load its package before solving with this algorithm.") end -function default_nev(prob::EigenvalueProblem) +function default_num_eigenpairs(prob::EigenvalueProblem) n = size(prob.A, 2) # Only the iterative backends call this; requesting the full dimension `n` # is invalid/degenerate for them, so default to a small subset. - return prob.nev === nothing ? min(n, 6) : prob.nev + return prob.num_eigenpairs === nothing ? min(n, 6) : prob.num_eigenpairs end -function _select_eigenpairs(values, vectors, nev, which, sigma) +function _select_eigenpairs(values, vectors, num_eigenpairs, eigentarget, shift) nvals = length(values) - howmany = nev === nothing ? nvals : min(nev, nvals) - ord = _eigenvalue_order(values, which, sigma) + howmany = num_eigenpairs === nothing ? nvals : min(num_eigenpairs, nvals) + ord = _eigenvalue_order(values, eigentarget, shift) idxs = ord[1:howmany] return values[idxs], vectors[:, idxs] end -# ARPACK-style symbol for the backends that take `which` as a `Symbol`. -function _target_symbol(w::EigenvalueTarget.T) - return w == EigenvalueTarget.LargestMagnitude ? :LM : - w == EigenvalueTarget.SmallestMagnitude ? :SM : - w == EigenvalueTarget.LargestRealPart ? :LR : - w == EigenvalueTarget.SmallestRealPart ? :SR : - w == EigenvalueTarget.LargestImaginaryPart ? :LI : - :SI -end - -function _eigenvalue_order(values, which::EigenvalueTarget.T, sigma) - if sigma !== nothing - return sortperm(abs.(values .- sigma)) - elseif which == EigenvalueTarget.LargestMagnitude +function _eigenvalue_order(values, eigentarget::EigenvalueTarget.T, shift) + if shift !== nothing + return sortperm(abs.(values .- shift)) + elseif eigentarget == EigenvalueTarget.LargestMagnitude return sortperm(abs.(values); rev = true) - elseif which == EigenvalueTarget.SmallestMagnitude + elseif eigentarget == EigenvalueTarget.SmallestMagnitude return sortperm(abs.(values)) - elseif which == EigenvalueTarget.LargestRealPart + elseif eigentarget == EigenvalueTarget.LargestRealPart return sortperm(real.(values); rev = true) - elseif which == EigenvalueTarget.SmallestRealPart + elseif eigentarget == EigenvalueTarget.SmallestRealPart return sortperm(real.(values)) - elseif which == EigenvalueTarget.LargestImaginaryPart + elseif eigentarget == EigenvalueTarget.LargestImaginaryPart return sortperm(imag.(values); rev = true) else # EigenvalueTarget.SmallestImaginaryPart return sortperm(imag.(values)) diff --git a/test/Core/eigenvalue.jl b/test/Core/eigenvalue.jl index 86a82c210..438bb875a 100644 --- a/test/Core/eigenvalue.jl +++ b/test/Core/eigenvalue.jl @@ -14,87 +14,105 @@ sol = solve(prob) @test sol.vectors * Diagonal(sol.u) ≈ A * sol.vectors @test sol.retcode === ReturnCode.Success -prob_nev = EigenvalueProblem(A; nev = 2, which = :SM) +prob_nev = EigenvalueProblem(A; num_eigenpairs = 2, eigentarget = EigenvalueTarget.SmallestMagnitude) sol_nev = solve(prob_nev) @test sol_nev.u == [1.0, 2.0] -prob_sigma = EigenvalueProblem(A; nev = 2, sigma = 2.2) -sol_sigma = solve(prob_sigma) -@test sol_sigma.u == [2.0, 3.0] +prob_shift = EigenvalueProblem(A; num_eigenpairs = 2, shift = 2.2) +sol_shift = solve(prob_shift) +@test sol_shift.u == [2.0, 3.0] B = Diagonal(fill(2.0, 4)) -prob_gen = EigenvalueProblem(A, B; nev = 2, which = :LR) +prob_gen = EigenvalueProblem(A, B; num_eigenpairs = 2, eigentarget = EigenvalueTarget.LargestRealPart) sol_gen = solve(prob_gen) @test sol_gen.u == [2.0, 1.5] -prob_uniform = EigenvalueProblem(A, 2I; nev = 2, which = :LR) +prob_uniform = EigenvalueProblem(A, 2I; num_eigenpairs = 2, eigentarget = EigenvalueTarget.LargestRealPart) sol_uniform = solve(prob_uniform) @test sol_uniform.u == [2.0, 1.5] A_backend = Diagonal(1.0:8.0) B_backend = Diagonal(fill(2.0, 8)) -sol_arpack = solve(EigenvalueProblem(Matrix(A_backend); nev = 2), ArpackJL()) +sol_arpack = solve(EigenvalueProblem(Matrix(A_backend); num_eigenpairs = 2), ArpackJL()) @test sol_arpack.u ≈ [8.0, 7.0] sol_arpack_gen = solve( - EigenvalueProblem(Matrix(A_backend), Matrix(B_backend); nev = 2), ArpackJL() + EigenvalueProblem(Matrix(A_backend), Matrix(B_backend); num_eigenpairs = 2), ArpackJL() ) @test sol_arpack_gen.u ≈ [4.0, 3.5] -sol_arnoldi = solve(EigenvalueProblem(A_backend; nev = 2), LinearSolve.ArnoldiMethod()) +sol_arnoldi = solve( + EigenvalueProblem(A_backend; num_eigenpairs = 2), LinearSolve.ArnoldiMethod() +) @test sol_arnoldi.u ≈ [8.0, 7.0] sol_arnoldi_default = solve(EigenvalueProblem(A_backend), LinearSolve.ArnoldiMethod()) @test length(sol_arnoldi_default.u) == 6 sol_arnoldi_shift = solve( - EigenvalueProblem(A_backend; nev = 2, sigma = 3.2), LinearSolve.ArnoldiMethod() + EigenvalueProblem(A_backend; num_eigenpairs = 2, shift = 3.2), LinearSolve.ArnoldiMethod() ) @test sol_arnoldi_shift.u ≈ [3.0, 4.0] -sol_krylovkit = solve(EigenvalueProblem(A_backend; nev = 2), KrylovKitEigen()) +# ArnoldiMethod has no smallest-magnitude target of its own. +@test_throws ArgumentError solve( + EigenvalueProblem( + A_backend; num_eigenpairs = 2, eigentarget = EigenvalueTarget.SmallestMagnitude + ), + LinearSolve.ArnoldiMethod() +) + +sol_krylovkit = solve(EigenvalueProblem(A_backend; num_eigenpairs = 2), KrylovKitEigen()) @test sol_krylovkit.u ≈ [8.0, 7.0] sol_krylovkit_shift = solve( - EigenvalueProblem(Matrix(A_backend); nev = 2, sigma = 3.2), KrylovKitEigen() + EigenvalueProblem(Matrix(A_backend); num_eigenpairs = 2, shift = 3.2), KrylovKitEigen() ) @test sol_krylovkit_shift.u ≈ [3.0, 4.0] sol_krylovkit_gen_shift = solve( - EigenvalueProblem(Matrix(A_backend), Matrix(B_backend); nev = 2, sigma = 1.6), + EigenvalueProblem(Matrix(A_backend), Matrix(B_backend); num_eigenpairs = 2, shift = 1.6), KrylovKitEigen() ) @test sol_krylovkit_gen_shift.u ≈ [1.5, 2.0] # Jacobi-Davidson is a target/interior method: it finds the eigenvalues nearest -# the target shift. `sigma` gives the interior target directly, and `:SM`/`:SR` -# work with the implicit zero shift. (`:LM`/`:LR` need a `sigma` guess.) +# the target shift. `shift` gives the interior target directly, and +# `SmallestMagnitude`/`SmallestRealPart` work with the implicit zero shift. +# (`LargestMagnitude`/`LargestRealPart` need a `shift` guess.) A_jd = Matrix(Diagonal(Float64.(1:30))) -sol_jd = solve(EigenvalueProblem(A_jd; nev = 2, sigma = 10.3), JacobiDavidsonJL()) +sol_jd = solve(EigenvalueProblem(A_jd; num_eigenpairs = 2, shift = 10.3), JacobiDavidsonJL()) @test sort(real(sol_jd.u)) ≈ [10.0, 11.0] @test sol_jd.retcode === ReturnCode.Success @test norm(A_jd * sol_jd.vectors - sol_jd.vectors * Diagonal(sol_jd.u)) < 1.0e-6 -sol_jd_sm = solve(EigenvalueProblem(A_jd; nev = 1, which = :SM), JacobiDavidsonJL()) +sol_jd_sm = solve( + EigenvalueProblem( + A_jd; num_eigenpairs = 1, eigentarget = EigenvalueTarget.SmallestMagnitude + ), + JacobiDavidsonJL() +) @test real(sol_jd_sm.u[1]) ≈ 1.0 # Generalized problems are not supported by the JacobiDavidson backend (upstream # `jdqz` is broken); it should error rather than silently misbehave. @test_throws ErrorException solve( - EigenvalueProblem(A_jd, Matrix(Diagonal(fill(2.0, 30))); nev = 1), JacobiDavidsonJL() + EigenvalueProblem(A_jd, Matrix(Diagonal(fill(2.0, 30))); num_eigenpairs = 1), + JacobiDavidsonJL() ) -# `which` accepts the EigenvalueTarget enum directly, equivalent to its Symbol alias. -prob_enum = EigenvalueProblem(A; nev = 2, which = EigenvalueTarget.SmallestMagnitude) -@test prob_enum.which === EigenvalueTarget.SmallestMagnitude +# `eigentarget` only accepts `EigenvalueTarget` values, not ARPACK-style symbols. +prob_enum = EigenvalueProblem( + A; num_eigenpairs = 2, eigentarget = EigenvalueTarget.SmallestMagnitude +) +@test prob_enum.eigentarget === EigenvalueTarget.SmallestMagnitude @test solve(prob_enum).u == [1.0, 2.0] -@test EigenvalueProblem(A; which = :SM).which === EigenvalueTarget.SmallestMagnitude -# An invalid `which` is rejected at construction. -@test_throws ArgumentError EigenvalueProblem(A; which = :XX) -@test_throws ArgumentError EigenvalueProblem(A; which = 5) +# An invalid `eigentarget` is rejected at construction with a clear type error. +@test_throws TypeError EigenvalueProblem(A; eigentarget = :LM) +@test_throws TypeError EigenvalueProblem(A; eigentarget = 5) # The iterative algorithms are keyword-only: positional arguments error, and # forwarded keywords reach the backend. From c9f0beadf71ebfeca714c9d42f679af94ee9deb7 Mon Sep 17 00:00:00 2001 From: Jash Date: Wed, 8 Jul 2026 11:45:24 +0530 Subject: [PATCH 6/7] =?UTF-8?q?fix=20LTS:=20rename=20EigenvalueSolution=20?= =?UTF-8?q?type=20param=20T=E2=86=92Tv=20to=20avoid=20ambiguity=20with=20R?= =?UTF-8?q?eturnCode.T=20field=20annotation?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/eigenvalue.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/eigenvalue.jl b/src/eigenvalue.jl index 993b2461d..1c4f75f47 100644 --- a/src/eigenvalue.jl +++ b/src/eigenvalue.jl @@ -76,7 +76,7 @@ else }(A, B, num_eigenpairs, eigentarget, shift, u0, p, kwargs) end - struct EigenvalueSolution{T, N, U, V, P, A, R, S} <: SciMLBase.AbstractNoTimeSolution{T, N} + struct EigenvalueSolution{Tv, N, U, V, P, A, R, S} <: SciMLBase.AbstractNoTimeSolution{Tv, N} u::U vectors::V prob::P @@ -90,10 +90,10 @@ else prob, alg, values, vectors; retcode = ReturnCode.Success, resid = nothing, stats = nothing ) - T = eltype(eltype(values)) + Tv = eltype(eltype(values)) N = length((size(values)...,)) return EigenvalueSolution{ - T, N, typeof(values), typeof(vectors), typeof(prob), typeof(alg), + Tv, N, typeof(values), typeof(vectors), typeof(prob), typeof(alg), typeof(resid), typeof(stats), }(values, vectors, prob, alg, retcode, resid, stats) end From 1426a702e333c9b744fbe8f5b25d9d5606dd75c9 Mon Sep 17 00:00:00 2001 From: Jash Date: Wed, 8 Jul 2026 13:02:45 +0530 Subject: [PATCH 7/7] updated compat --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 4a5e87410..1eccf9bd3 100644 --- a/Project.toml +++ b/Project.toml @@ -121,6 +121,8 @@ AMDGPU = "1.2, 2" # the ext API is identical and the AMG testset passes on 0.5.1. AlgebraicMultigrid = "0.5, 1, 2" ArrayInterface = "7.19" +ArnoldiMethod = "0.4" +Arpack = "0.5" BandedMatrices = "1.8" BlockDiagonals = "0.2" cuSOLVER = "6"