Move eigs from ApproxFun (#458)

* Move eigs from ApproxFun

* update docstring

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunBase"
 uuid = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
-version = "0.8.19"
+version = "0.8.20"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

src/ApproxFunBase.jl

@@ -43,7 +43,7 @@ import LinearAlgebra: BlasInt, BlasFloat, norm, ldiv!, mul!, det, cross,
               qr, qr!, rank, isdiag, istril, istriu, issymmetric,
               Tridiagonal, diagm, diagm_container, factorize,
               nullspace, Hermitian, Symmetric, adjoint, transpose, char_uplo,
-              axpy!
+              axpy!, eigvals
 
 import SparseArrays: blockdiag
 
@@ -151,6 +151,7 @@ include("Operators/Operator.jl")
 include("Caching/caching.jl")
 include("PDE/PDE.jl")
 include("Spaces/Spaces.jl")
+include("eigen.jl")
 include("hacks.jl")
 include("testing.jl")
 include("specialfunctions.jl")

src/eigen.jl

@@ -0,0 +1,93 @@
+export eigs
+
+
+eigvals(A::Operator,n::Integer;tolerance::Float64=100eps()) =
+    eigs(A,n;tolerance=tolerance)[1]
+
+"""
+    λ, V = eigs(A::Operator, n::Integer; tolerance::Float64=100eps())
+
+Compute the eigenvalues and eigenvectors of the operator `A`. This is done in the following way:
+
+* Truncate `A` into an `n×n` matrix `A₁`.
+* Compute eigenvalues and eigenvectors of `A₁`.
+* Filter for those eigenvectors of `A₁` that are approximately eigenvectors
+of `A` as well. The `tolerance` argument controls, which eigenvectors of the approximation are kept.
+"""
+function eigs(A::Operator,n::Integer;tolerance::Float64=100eps())
+    typ = eltype(A)
+
+    ds=domainspace(A)
+    C = Conversion(ds,rangespace(A))
+
+    A1,C1 = zeros(typ,n,n),zeros(typ,n,n)
+    A1[1:end,1:n] = A[1:n,1:n]
+    C1[1:end,1:n] = C[1:n,1:n]
+
+    λ,V=eigen(A1,C1)
+
+    pruneeigs(λ,V,ds,tolerance)
+end
+
+eigvals(Bcs::Operator,A::Operator,n::Integer;tolerance::Float64=100eps()) =
+    eigs(Bcs,A,n;tolerance=tolerance)[1]
+
+"""
+    λ, V = eigs(BC::Operator, A::Operator, n::Integer; tolerance::Float64=100eps())
+
+Compute `n` eigenvalues and eigenvectors of the operator `A`,
+subject to the boundary conditions `BC`.
+
+# Examples
+```jldoctest
+julia> #= We compute the spectrum of the second derivative,
+          subject to zero boundary conditions.
+          We solve this eigenvalue problem in the Chebyshev basis =#
+
+julia> S = Chebyshev();
+
+julia> D = Derivative(S, 2);
+
+julia> BC = Dirichlet(S);
+
+julia> λ, v = ApproxFun.eigs(BC, D, 100);
+
+julia> λ[1:10] ≈ [-(n*pi/2)^2 for n in 1:10] # compare with the analytical result
+true
+```
+"""
+function eigs(Bcs_in::Operator,A_in::Operator,n::Integer;tolerance::Float64=100eps())
+    Bcs, A = promotedomainspace([Bcs_in, A_in])
+
+    nf = size(Bcs,1)
+    @assert isfinite(nf)
+
+    typ = promote_type(eltype(Bcs),eltype(A))
+
+    ds=domainspace(A)
+    C = Conversion(ds,rangespace(A))
+
+    A1,C1 = zeros(typ,n,n),zeros(typ,n,n)
+    A1[1:nf,1:n] = Bcs[1:nf,1:n]
+    A1[nf+1:end,1:n] = A[1:n-nf,1:n]
+    C1[nf+1:end,1:n] = C[1:n-nf,1:n]
+
+    λ,V = eigen(A1,C1)
+
+    λ, V = pruneeigs(λ,V,ds,tolerance)
+    p = sortperm(λ; lt=(x,y) -> isless(abs(x),abs(y)))
+    λ[p], V[p]
+end
+
+function pruneeigs(λ,V,ds,tolerance)
+    retλ=eltype(λ)[]
+    retV=VFun{typeof(ds),eltype(V)}[]
+    n=length(λ)
+    for k=1:n
+        if slnorm(V,n-3:n,k)≤tolerance
+            push!(retλ,λ[k])
+            push!(retV,Fun(ds,V[:,k]))
+        end
+    end
+    retλ,retV
+end