From c5606b5e1c6b2455541b68e26db755d15337ef0e Mon Sep 17 00:00:00 2001 From: Jash Date: Fri, 3 Jul 2026 21:09:11 +0530 Subject: [PATCH 01/11] 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 02/11] 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 03/11] 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 04/11] 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 05/11] 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 06/11] =?UTF-8?q?fix=20LTS:=20rename=20EigenvalueSolution?= =?UTF-8?q?=20type=20param=20T=E2=86=92Tv=20to=20avoid=20ambiguity=20with?= =?UTF-8?q?=20ReturnCode.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 07/11] 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" From 1455ca5023a090a5fce96dbd5c1d089308c0c416 Mon Sep 17 00:00:00 2001 From: Jash Date: Sat, 11 Jul 2026 09:43:54 +0530 Subject: [PATCH 08/11] added documentation and fixed QA failure --- Project.toml | 2 +- docs/pages.jl | 4 +- docs/src/basics/EigenvalueProblem.md | 6 + docs/src/solvers/eigenvalue_solvers.md | 77 +++++++++++ docs/src/tutorials/eigenvalue.md | 72 ++++++++++ src/eigenvalue.jl | 173 ++++++++++--------------- test/qa/qa.jl | 6 +- 7 files changed, 234 insertions(+), 106 deletions(-) create mode 100644 docs/src/basics/EigenvalueProblem.md create mode 100644 docs/src/solvers/eigenvalue_solvers.md create mode 100644 docs/src/tutorials/eigenvalue.md diff --git a/Project.toml b/Project.toml index 1eccf9bd3..b73aca7f4 100644 --- a/Project.toml +++ b/Project.toml @@ -178,7 +178,7 @@ RecursiveFactorization = "0.2.26" Reexport = "1.2.2" STRUMPACK_jll = "8" SafeTestsets = "0.1" -SciMLBase = "2.148, 3" +SciMLBase = "3.32" SciMLLogging = "1.10.1, 2.0" SciMLOperators = "1.15" SciMLTesting = "1" diff --git a/docs/pages.jl b/docs/pages.jl index c39e6250b..8ee13f0b8 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -11,16 +11,18 @@ pages = [ "tutorials/accelerating_choices.md", "tutorials/gpu.md", "tutorials/autotune.md", + "tutorials/eigenvalue.md", ], "Basics" => Any[ "basics/LinearProblem.md", + "basics/EigenvalueProblem.md", "basics/algorithm_selection.md", "basics/common_solver_opts.md", "basics/OperatorAssumptions.md", "basics/Preconditioners.md", "basics/FAQ.md", ], - "Solvers" => Any["solvers/solvers.md"], + "Solvers" => Any["solvers/solvers.md", "solvers/eigenvalue_solvers.md"], "Advanced" => Any[ "advanced/developing.md", "advanced/custom.md", diff --git a/docs/src/basics/EigenvalueProblem.md b/docs/src/basics/EigenvalueProblem.md new file mode 100644 index 000000000..c6a965ce3 --- /dev/null +++ b/docs/src/basics/EigenvalueProblem.md @@ -0,0 +1,6 @@ +# Eigenvalue Problems + +```@docs +EigenvalueProblem +EigenvalueTarget +``` diff --git a/docs/src/solvers/eigenvalue_solvers.md b/docs/src/solvers/eigenvalue_solvers.md new file mode 100644 index 000000000..0f91dd40c --- /dev/null +++ b/docs/src/solvers/eigenvalue_solvers.md @@ -0,0 +1,77 @@ +# [Eigenvalue Problem Solvers](@id eigenvaluesolvers) + +`LS.solve(prob::LS.EigenvalueProblem, alg; kwargs)` + +Solves for the eigenpairs $$(\lambda, v)$$ of $$Av = \lambda v$$ (or the generalized +problem $$Av = \lambda Bv$$ when a second operator is supplied) defined by `prob` using +the algorithm `alg`. If no algorithm is given, `DenseEigen()` is used. + +## Recommended Methods + +The default `DenseEigen()` wraps `LinearAlgebra.eigen` and is the right choice for +small to moderately sized dense matrices, or whenever most/all of the spectrum is +needed. It always computes the full dense eigendecomposition before selecting the +requested eigenpairs, so its cost does not depend on `num_eigenpairs`. + +For large sparse or structured matrices where only a handful of eigenpairs are needed, +an iterative Krylov-based backend is preferred, since these only ever factorize/apply +`A` (and `B`) and never form the dense decomposition: + + - `KrylovKitEigen()` is a solid default iterative choice: it supports both standard + and generalized problems, extremal targets, and interior (`shift`) targets. + - `ArpackJL()` is a mature, widely-used implicitly restarted Arnoldi/Lanczos solver, + also supporting extremal and shifted targets on standard and generalized problems. + - `LS.ArnoldiMethod()` is a pure-Julia implicitly restarted Arnoldi method. It has no + direct `EigenvalueTarget.SmallestMagnitude` support; use `shift` (shift-and-invert) + or another backend for interior/smallest eigenvalues instead. + - `JacobiDavidsonJL()` is a target/interior method: it is most effective when a + `shift` close to the desired eigenvalues is known. It does not support generalized + eigenvalue problems. + +## Full List of Methods + +### LinearAlgebra (built-in) + +```@docs +DenseEigen +``` + +### Arpack.jl + +!!! note + + Using this solver requires adding the package Arpack.jl, i.e. `using Arpack` + +```@docs +ArpackJL +``` + +### ArnoldiMethod.jl + +!!! note + + Using this solver requires adding the package ArnoldiMethod.jl, i.e. `using ArnoldiMethod` + +```@docs +ArnoldiMethod +``` + +### KrylovKit.jl + +!!! note + + Using this solver requires adding the package KrylovKit.jl, i.e. `using KrylovKit` + +```@docs +KrylovKitEigen +``` + +### JacobiDavidson.jl + +!!! note + + Using this solver requires adding the package JacobiDavidson.jl, i.e. `using JacobiDavidson` + +```@docs +JacobiDavidsonJL +``` diff --git a/docs/src/tutorials/eigenvalue.md b/docs/src/tutorials/eigenvalue.md new file mode 100644 index 000000000..f1eff4521 --- /dev/null +++ b/docs/src/tutorials/eigenvalue.md @@ -0,0 +1,72 @@ +# Getting Started with Solving Eigenvalue Problems in Julia + +An eigenvalue problem $$Av = \lambda v$$ (or, generalized, $$Av = \lambda Bv$$) is +specified by defining an `EigenvalueProblem`. As with `LinearProblem`, the same +`EigenvalueProblem` can be solved with different algorithms just by swapping out the +`alg` argument to `solve`. + +```@example eigsys1 +import LinearSolve as LS +import LinearAlgebra + +A = LinearAlgebra.Diagonal([1.0, 2.0, 3.0, 4.0]) +prob = LS.EigenvalueProblem(A) +sol = LS.solve(prob) +sol.u +``` + +`sol.u` holds the eigenvalues and `sol.vectors` holds the corresponding eigenvectors as +the columns of a matrix, satisfying `A * sol.vectors ≈ sol.vectors * Diagonal(sol.u)`. +With no algorithm specified, `LS.DenseEigen()` is used, which wraps +`LinearAlgebra.eigen` and then selects the requested eigenpairs out of the full +decomposition. + +## Requesting a Subset of Eigenpairs + +For large problems, computing every eigenpair is wasteful when only a handful are +needed. `num_eigenpairs` controls how many eigenpairs are returned, and `eigentarget` +(an [`EigenvalueTarget`](@ref)) controls which part of the spectrum they come from: + +```@example eigsys1 +prob = LS.EigenvalueProblem( + A; num_eigenpairs = 2, eigentarget = LS.EigenvalueTarget.SmallestMagnitude +) +sol = LS.solve(prob) +sol.u +``` + +Alternatively, `shift` requests the eigenvalues nearest a given value +(shift-and-invert): + +```@example eigsys1 +prob = LS.EigenvalueProblem(A; num_eigenpairs = 2, shift = 2.2) +sol = LS.solve(prob) +sol.u +``` + +## Generalized Eigenvalue Problems + +Passing a second operator `B` solves the generalized problem $$Av = \lambda Bv$$ +instead: + +```@example eigsys1 +B = LinearAlgebra.Diagonal(fill(2.0, 4)) +prob = LS.EigenvalueProblem(A, B; num_eigenpairs = 2, eigentarget = LS.EigenvalueTarget.LargestRealPart) +sol = LS.solve(prob) +sol.u +``` + +## Switching to an Iterative Backend + +For large sparse or structured matrices, an iterative Krylov-based backend is often +preferable to the dense decomposition. Just as with `LinearProblem`, changing the +solver is a matter of swapping out the algorithm passed to `solve`: + +```julia +using KrylovKit +sol = LS.solve(LS.EigenvalueProblem(A; num_eigenpairs = 2), LS.KrylovKitEigen()) +sol.u +``` + +For more information on the available eigenvalue backends and when to choose each one, +see [the eigenvalue solvers page](@ref eigenvaluesolvers). diff --git a/src/eigenvalue.jl b/src/eigenvalue.jl index 1c4f75f47..f8373088f 100644 --- a/src/eigenvalue.jl +++ b/src/eigenvalue.jl @@ -1,129 +1,100 @@ -# `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 `num_eigenpairs`) in an [`EigenvalueProblem`](@ref). - """ - EnumX.@enumx EigenvalueTarget begin - "Eigenvalues of largest magnitude, `abs(λ)` largest." - LargestMagnitude - "Eigenvalues of smallest magnitude, `abs(λ)` smallest." - SmallestMagnitude - "Eigenvalues with the largest (most positive) real part." - LargestRealPart - "Eigenvalues with the smallest (most negative) real part." - SmallestRealPart - "Eigenvalues with the largest (most positive) imaginary part." - LargestImaginaryPart - "Eigenvalues with the smallest (most negative) imaginary part." - SmallestImaginaryPart - end - - """ - EigenvalueProblem(A[, B], p = SciMLBase.NullParameters(); - num_eigenpairs = nothing, eigentarget = EigenvalueTarget.LargestMagnitude, - shift = 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 - - - `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, TargetType, ShiftType, U0Type, PType, KType, - } - A::AType - B::BType - num_eigenpairs::NevType - eigentarget::TargetType - shift::ShiftType - u0::U0Type - p::PType - kwargs::KType - end +using SciMLBase: EigenvalueProblem, EigenvalueSolution, EigenvalueTarget, + build_eigenvalue_solution - function EigenvalueProblem( - A, B = nothing, p = SciMLBase.NullParameters(); - num_eigenpairs = nothing, - eigentarget::EigenvalueTarget.T = EigenvalueTarget.LargestMagnitude, - shift = nothing, u0 = nothing, kwargs... - ) - return EigenvalueProblem{ - 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{Tv, N, U, V, P, A, R, S} <: SciMLBase.AbstractNoTimeSolution{Tv, 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 - ) - Tv = eltype(eltype(values)) - N = length((size(values)...,)) - return EigenvalueSolution{ - Tv, 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 +""" + AbstractEigenvalueAlgorithm +Base type for algorithms that solve an [`EigenvalueProblem`](@ref). +""" abstract type AbstractEigenvalueAlgorithm <: SciMLBase.AbstractLinearAlgorithm end +""" + DenseEigen() + +Solve the `EigenvalueProblem` with `LinearAlgebra.eigen`. This is the default +algorithm: it computes the full dense eigendecomposition and then selects the +requested eigenpairs (via `num_eigenpairs`, `eigentarget`, or `shift`) from it. +Best for small to moderately sized dense matrices where every eigenpair (or a +sizable fraction of them) is needed. +""" struct DenseEigen <: AbstractEigenvalueAlgorithm end # 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. +""" + ArpackJL(; kwargs...) + +Solve the `EigenvalueProblem` with [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl)'s +`eigs`, an iterative Krylov (implicitly restarted Arnoldi/Lanczos) solver well suited to +computing a handful of extremal eigenpairs of a large sparse or structured matrix. +Extra `kwargs` are forwarded to `Arpack.eigs`. + +!!! note + + Using this solver requires loading Arpack.jl, i.e. `using Arpack`. +""" struct ArpackJL{K <: NamedTuple} <: AbstractEigenvalueAlgorithm kwargs::K end ArpackJL(; kwargs...) = ArpackJL((; kwargs...)) +""" + ArnoldiMethodJL + +Algorithm type constructed by [`ArnoldiMethod`](@ref); see its docstring for details. +""" struct ArnoldiMethodJL{K <: NamedTuple} <: AbstractEigenvalueAlgorithm kwargs::K end + +""" + ArnoldiMethod(; kwargs...) + +Solve the `EigenvalueProblem` with +[ArnoldiMethod.jl](https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl)'s +`partialschur`, a pure-Julia implicitly restarted Arnoldi method for large sparse or +structured matrices. Does not support `eigentarget = EigenvalueTarget.SmallestMagnitude` +directly; use `shift` for shift-and-invert instead, or another backend. Extra `kwargs` +are forwarded to `ArnoldiMethod.partialschur`. + +!!! note + + Using this solver requires loading ArnoldiMethod.jl, i.e. `using ArnoldiMethod`. +""" ArnoldiMethod(; kwargs...) = ArnoldiMethodJL((; kwargs...)) +""" + KrylovKitEigen(; kwargs...) + +Solve the `EigenvalueProblem` with +[KrylovKit.jl](https://github.com/Jutho/KrylovKit.jl)'s `eigsolve`, a Krylov solver +supporting both standard and generalized eigenvalue problems, extremal and interior +(shifted) targets. Extra `kwargs` are forwarded to `KrylovKit.eigsolve`. + +!!! note + + Using this solver requires loading KrylovKit.jl, i.e. `using KrylovKit`. +""" struct KrylovKitEigen{K <: NamedTuple} <: AbstractEigenvalueAlgorithm kwargs::K end KrylovKitEigen(; kwargs...) = KrylovKitEigen((; kwargs...)) +""" + JacobiDavidsonJL(; kwargs...) + +Solve the `EigenvalueProblem` with +[JacobiDavidson.jl](https://github.com/haampie/JacobiDavidson.jl)'s `jdqr`, a +target/interior method that finds the eigenvalues nearest a given `shift`. Does not +support generalized eigenvalue problems (upstream `jdqz` is broken). Extra `kwargs` are +forwarded to `JacobiDavidson.jdqr`. + +!!! note + + Using this solver requires loading JacobiDavidson.jl, i.e. `using JacobiDavidson`. +""" struct JacobiDavidsonJL{K <: NamedTuple} <: AbstractEigenvalueAlgorithm kwargs::K end diff --git a/test/qa/qa.jl b/test/qa/qa.jl index e355aa79b..edf7877f7 100644 --- a/test/qa/qa.jl +++ b/test/qa/qa.jl @@ -32,7 +32,7 @@ run_qa( aqua_kwargs = (; deps_compat = (; ignore = [:MKL_jll]), stale_deps = (; ignore = [:MKL_jll]), - piracies = (; treat_as_own = [LinearProblem]), + piracies = (; treat_as_own = [LinearProblem, EigenvalueProblem]), ), ei_kwargs = (; no_implicit_imports = (; @@ -59,8 +59,8 @@ run_qa( ignore = ( Symbol("@blasfunc"), :AbstractSciMLOperator, :AbstractSparseMatrixCSC, :ArrayInterface, :BLASELTYPES, :BlasInt, :StaticArray, :UMFPACK_OK, - :chkargsok, :chkfinite, :chkstride1, :getcolptr, :inv, - :pattern_changed, :require_one_based_indexing, + :build_eigenvalue_solution, :chkargsok, :chkfinite, :chkstride1, + :getcolptr, :inv, :pattern_changed, :require_one_based_indexing, ), ), ), From a6200b29e5bb99b929501dd36059f4d884abe173 Mon Sep 17 00:00:00 2001 From: Jash Date: Sat, 11 Jul 2026 10:07:57 +0530 Subject: [PATCH 09/11] added StableRNGs in test targets --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index b73aca7f4..5457168da 100644 --- a/Project.toml +++ b/Project.toml @@ -230,4 +230,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -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"] +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", "StableRNGs", "StaticArrays", "Test", "Zygote"] From 2794af30db3c6a4b3a104b80ddd98d6aecf80021 Mon Sep 17 00:00:00 2001 From: Jash Date: Sat, 11 Jul 2026 11:11:47 +0530 Subject: [PATCH 10/11] resolve merge conflict --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 5457168da..03296799d 100644 --- a/Project.toml +++ b/Project.toml @@ -188,6 +188,7 @@ SparseColumnPivotedQR = "2" SparseMatricesCSR = "0.6" Sparspak = "0.3.9" SpecializingFactorizations = "0.1" +StableRNGs = "1" StaticArrays = "1.9" StaticArraysCore = "1.4.3" SuperLUDIST = "1" @@ -225,6 +226,7 @@ SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac" SpecializingFactorizations = "fa08b7a1-13d3-4faf-875d-5cbc1520e3f3" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" From 6b8c126b09905d56ac78bad59aff2ffca0170433 Mon Sep 17 00:00:00 2001 From: Jash Date: Sat, 11 Jul 2026 12:24:04 +0530 Subject: [PATCH 11/11] fix --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7639bca60..29045d5a4 100644 --- a/Project.toml +++ b/Project.toml @@ -187,7 +187,6 @@ SparseArrays = "1.10" SparseColumnPivotedQR = "2" SparseMatricesCSR = "0.6" Sparspak = "0.3.9" -StableRNGs = "1" SpecializingFactorizations = "0.1" StableRNGs = "1" StaticArrays = "1.9"