diff --git a/Project.toml b/Project.toml
index 604b0a187..29045d5a4 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"
@@ -169,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"
@@ -178,8 +187,8 @@ SparseArrays = "1.10"
 SparseColumnPivotedQR = "2"
 SparseMatricesCSR = "0.6"
 Sparspak = "0.3.9"
-StableRNGs = "1"
 SpecializingFactorizations = "0.1"
+StableRNGs = "1"
 StaticArrays = "1.9"
 StaticArraysCore = "1.4.3"
 SuperLUDIST = "1"
@@ -190,6 +199,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"
@@ -201,6 +212,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"
@@ -220,4 +232,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", "StableRNGs", "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"]
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/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 224531b17..a402e159a 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")
@@ -574,4 +575,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..f8373088f
--- /dev/null
+++ b/src/eigenvalue.jl
@@ -0,0 +1,162 @@
+using SciMLBase: EigenvalueProblem, EigenvalueSolution, EigenvalueTarget,
+    build_eigenvalue_solution
+
+"""
+    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
+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/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,
             ),
         ),
     ),
diff --git a/test/runtests.jl b/test/runtests.jl
index c44370b6a..50fa3629e 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 "GESV Factorization" include("Core/gesv.jl")
             @time @safetestset "Return codes" include("Core/retcodes.jl")