diff --git a/Project.toml b/Project.toml
index d9d899663..1eccf9bd3 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"
@@ -115,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"
@@ -139,6 +147,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"
@@ -189,6 +198,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"
@@ -200,6 +211,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"
@@ -218,4 +230,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..596c1c0d7
--- /dev/null
+++ b/ext/LinearSolveArnoldiMethodExt.jl
@@ -0,0 +1,66 @@
+module LinearSolveArnoldiMethodExt
+
+using LinearAlgebra
+using LinearSolve
+import ArnoldiMethod: partialschur, partialeigen, LM, LR, SR, LI, SI
+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_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.shift !== nothing
+        values = prob.shift .+ inv.(values)
+    end
+    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
+
+# 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
+
+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..9d02caa72
--- /dev/null
+++ b/ext/LinearSolveArpackExt.jl
@@ -0,0 +1,45 @@
+module LinearSolveArpackExt
+
+using LinearSolve
+using Arpack
+using SciMLBase: SciMLBase, ReturnCode
+
+function SciMLBase.solve(
+        prob::LinearSolve.EigenvalueProblem,
+        alg::LinearSolve.ArpackJL,
+        args...; kwargs...
+    )
+    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
+    # as a keyword argument.
+    values, vectors, nconv, niter, nmult, resid = if prob.B === nothing
+        Arpack.eigs(prob.A; kw...)
+    else
+        Arpack.eigs(prob.A, prob.B; 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
+
+# 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
new file mode 100644
index 000000000..7d4f617a9
--- /dev/null
+++ b/ext/LinearSolveJacobiDavidsonExt.jl
@@ -0,0 +1,68 @@
+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_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.
+    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; kw...)
+    values, vectors = _jd_standard_pairs(prob.A, out[1])
+
+    values, vectors = LinearSolve._select_eigenpairs(
+        values, vectors, nev, prob.eigentarget, prob.shift
+    )
+    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
+# `shift` is the natural interior target (`Near`), which is Jacobi-Davidson's
+# strength; otherwise `eigentarget` selects an extremal target.
+function _jd_target(prob)
+    prob.shift !== nothing && return JacobiDavidson.Near(ComplexF64(prob.shift))
+    T = LinearSolve.EigenvalueTarget
+    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) :
+        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
+# 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..0ad0d3b23 100644
--- a/ext/LinearSolveKrylovKitExt.jl
+++ b/ext/LinearSolveKrylovKitExt.jl
@@ -61,4 +61,61 @@ 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_num_eigenpairs(prob)
+    which = _krylovkit_which(prob.eigentarget)
+    kw = (; prob.kwargs..., alg.kwargs..., kwargs...)
+    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.shift !== nothing
+        values = prob.shift .+ inv.(values)
+    end
+    vecmat = reduce(hcat, vectors)
+    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, 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 - shift * I)
+        op = x -> F \ x
+    elseif B isa LinearAlgebra.UniformScaling
+        F = factorize(A - shift * B)
+        op = x -> F \ (B.λ * x)
+    else
+        F = factorize(A - shift * B)
+        op = x -> F \ (B * x)
+    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 292fbd468..ebe192bca 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 AbstractEigenvalueAlgorithm,
+    DenseEigen, ArpackJL, ArnoldiMethod, ArnoldiMethodJL,
+    KrylovKitEigen, JacobiDavidsonJL
+
 end
diff --git a/src/eigenvalue.jl b/src/eigenvalue.jl
new file mode 100644
index 000000000..1c4f75f47
--- /dev/null
+++ b/src/eigenvalue.jl
@@ -0,0 +1,191 @@
+# `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
+
+    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
+
+abstract type AbstractEigenvalueAlgorithm <: SciMLBase.AbstractLinearAlgorithm end
+
+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.
+struct ArpackJL{K <: NamedTuple} <: AbstractEigenvalueAlgorithm
+    kwargs::K
+end
+ArpackJL(; kwargs...) = ArpackJL((; kwargs...))
+
+struct ArnoldiMethodJL{K <: NamedTuple} <: AbstractEigenvalueAlgorithm
+    kwargs::K
+end
+ArnoldiMethod(; kwargs...) = ArnoldiMethodJL((; kwargs...))
+
+struct KrylovKitEigen{K <: NamedTuple} <: AbstractEigenvalueAlgorithm
+    kwargs::K
+end
+KrylovKitEigen(; kwargs...) = KrylovKitEigen((; kwargs...))
+
+struct JacobiDavidsonJL{K <: NamedTuple} <: AbstractEigenvalueAlgorithm
+    kwargs::K
+end
+JacobiDavidsonJL(; kwargs...) = JacobiDavidsonJL((; kwargs...))
+
+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.num_eigenpairs, prob.eigentarget, prob.shift
+    )
+    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_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.num_eigenpairs === nothing ? min(n, 6) : prob.num_eigenpairs
+end
+
+function _select_eigenpairs(values, vectors, num_eigenpairs, eigentarget, shift)
+    nvals = length(values)
+    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
+
+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 eigentarget == EigenvalueTarget.SmallestMagnitude
+        return sortperm(abs.(values))
+    elseif eigentarget == EigenvalueTarget.LargestRealPart
+        return sortperm(real.(values); rev = true)
+    elseif eigentarget == EigenvalueTarget.SmallestRealPart
+        return sortperm(real.(values))
+    elseif eigentarget == EigenvalueTarget.LargestImaginaryPart
+        return sortperm(imag.(values); rev = true)
+    else # EigenvalueTarget.SmallestImaginaryPart
+        return sortperm(imag.(values))
+    end
+end
diff --git a/test/Core/eigenvalue.jl b/test/Core/eigenvalue.jl
new file mode 100644
index 000000000..438bb875a
--- /dev/null
+++ b/test/Core/eigenvalue.jl
@@ -0,0 +1,120 @@
+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; num_eigenpairs = 2, eigentarget = EigenvalueTarget.SmallestMagnitude)
+sol_nev = solve(prob_nev)
+@test sol_nev.u == [1.0, 2.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; num_eigenpairs = 2, eigentarget = EigenvalueTarget.LargestRealPart)
+sol_gen = solve(prob_gen)
+@test sol_gen.u == [2.0, 1.5]
+
+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); num_eigenpairs = 2), ArpackJL())
+@test sol_arpack.u ≈ [8.0, 7.0]
+
+sol_arpack_gen = solve(
+    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; 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; num_eigenpairs = 2, shift = 3.2), LinearSolve.ArnoldiMethod()
+)
+@test sol_arnoldi_shift.u ≈ [3.0, 4.0]
+
+# 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); 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); 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. `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; 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; 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))); num_eigenpairs = 1),
+    JacobiDavidsonJL()
+)
+
+# `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]
+
+# 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.
+@test_throws MethodError ArpackJL(5)
+@test ArpackJL(; maxiter = 500).kwargs.maxiter == 500
diff --git a/test/runtests.jl b/test/runtests.jl
index cd7817b9d..3b8c9be4a 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 "Batched RHS" include("Core/batch.jl")
             @time @safetestset "Return codes" include("Core/retcodes.jl")
             @time @safetestset "Re-solve" include("Core/resolve.jl")