eigs for EigenSystem (#247)

* eigs for EigenSystem

* don't import iterate

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.6.27"
+version = "0.6.28"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"

src/ApproxFunOrthogonalPolynomials.jl

@@ -44,7 +44,7 @@ import ApproxFunBase: Fun, SubSpace, WeightSpace, NoSpace, HeavisideSpace,
                     ℓ⁰, recα, recβ, recγ, ℵ₀, ∞, RectDomain,
                     assert_integer, supportsinplacetransform, ContinuousSpace, SpecialEvalPtType,
                     isleftendpoint, isrightendpoint, evaluation_point, boundaryfn, ldiffbc, rdiffbc,
-                    LeftEndPoint, RightEndPoint, normalizedspace, promotedomainspace
+                    LeftEndPoint, RightEndPoint, normalizedspace, promotedomainspace, eigs
 
 import DomainSets: Domain, indomain, UnionDomain, FullSpace, Point,
             Interval, ChebyshevInterval, boundary, rightendpoint, leftendpoint,
@@ -82,7 +82,7 @@ import SpecialFunctions: erfcx, dawson,
 
 using StaticArrays: SVector
 
-import LinearAlgebra: isdiag, eigvals
+import LinearAlgebra: isdiag, eigvals, eigen
 
 points(d::IntervalOrSegmentDomain{T},n::Integer) where {T} =
     _maybefromcanonical(d, chebyshevpoints(float(real(eltype(T))), n))  # eltype to handle point

src/symeigen.jl

@@ -46,9 +46,9 @@ SkewSymmetricEigensystem
 
 for SET in (:SymmetricEigensystem, :SkewSymmetricEigensystem)
     @eval begin
-        struct $SET{LT,QT} <: EigenSystem
+        struct $SET{LT,QST} <: EigenSystem
             L :: LT
-            Q :: QT
+            QS :: QST
 
             function $SET(L, B, ::Type{QST} = QuotientSpace) where {QST}
                 L2, B2 = promotedomainspace((L, B))
@@ -57,19 +57,18 @@ for SET in (:SymmetricEigensystem, :SkewSymmetricEigensystem)
                 end
 
                 QS = QST(B2)
-                S = domainspace(L2)
-                Q = Conversion(QS, S)
-                new{typeof(L2),typeof(Q)}(L2, Q)
+                new{typeof(L2),typeof(QS)}(L2, QS)
             end
         end
     end
 end
 
 function basis_recombination(SE::EigenSystem)
-    L, Q = SE.L, SE.Q
+    L, QS = SE.L, SE.QS
     S = domainspace(L)
     D1 = Conversion_normalizedspace(S)
     D2 = Conversion_normalizedspace(S, Val(:backward))
+    Q = Conversion(QS, S)
     R = D1*Q;
     C = Conversion(S, rangespace(L))
     P = cache(PartialInverseOperator(C, (0, bandwidth(L, 1) + bandwidth(R, 1) + bandwidth(C, 2))));
@@ -116,3 +115,15 @@ function eigvals(S::EigenSystem, n::Integer)
     SA, SB = bandmatrices_eigen(S, n)
     eigvals(SA, SB)
 end
+
+function eigs(Seig::EigenSystem, n::Integer; tolerance::Float64=100eps())
+    SA, SB = bandmatrices_eigen(Seig, n)
+    λ, v = eigen(SA, SB)
+    vm = Matrix(v)
+    L, QS = Seig.L, Seig.QS
+    S = domainspace(L)
+    Q = Conversion(QS, S)
+    QM = Q[ApproxFunBase.FiniteRange, axes(vm, 1)]
+    V = QM * vm
+    ApproxFunBase.pruneeigs(λ,V,S,tolerance)
+end

test/EigTest.jl

@@ -109,6 +109,10 @@ using Test
         @test u2'(0.25) - u1'(0.25) ≈ 100*u(0.25)
         @test -u1'' + component(V, 1)*u1 ≈ λ[k]*u1
         @test -u2'' + component(V, 2)*u2 ≈ λ[k]*u2
+
+        λ2, f = eigs(Seig, n, tolerance=1e-8);
+        @test λ2[1:10] ≈ λ[1:10]
+        @test f[k] ≈ u || f[k] ≈ -u
     end
 
     @testset "BigFloat negative Laplacian with Dirichlet boundary conditions" begin