Add EigenvalueProblem support with dense, Arpack, ArnoldiMethod, KrylovKit, and JacobiDavidson backends

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"]

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",

docs/src/basics/EigenvalueProblem.md

@@ -0,0 +1,6 @@
+# Eigenvalue Problems
+
+```@docs
+EigenvalueProblem
+EigenvalueTarget
+```

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
+```

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).

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

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