Add extension for ArnoldiMethod and tests for cluster using it

Author: RalphAS

Project.toml

@@ -5,19 +5,27 @@ version = "0.2.2-DEV"
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
+[weakdeps]
+ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"
+
+[extensions]
+ArnoldiExt = "ArnoldiMethod"
+
 [compat]
 julia = "1"
 Aqua = "0.8"
+ArnoldiMethod = "0.4"
 LinearAlgebra = "1"
 Quadmath = "0.5"
 Random = "1"
 Test = "1"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"
 Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test","Random","Quadmath","Aqua"]
+test = ["Test","Random","Quadmath","ArnoldiMethod","Aqua"]

ext/ArnoldiExt.jl

@@ -0,0 +1,28 @@
+module ArnoldiExt
+
+using LinearAlgebra
+using IterativeRefinement
+const IR = IterativeRefinement
+using ArnoldiMethod
+
+"""
+    rfeigen(A, S::partialschur, idxλ, DT, maxiter=5) -> vals, vecs, status
+
+Uses a partial Schur decomposition from ArnoldiMethod for cluster refinement.
+"""
+function IR.rfeigen(A::AbstractMatrix{T}, ps::ArnoldiMethod.PartialSchur{TQ,TR}, idxλ,
+                    DT = widen(real(T)), maxiter=5; kwargs...
+) where {T, TQ, TR <: AbstractMatrix{TRe}} where {TRe <: Complex}
+    S = LinearAlgebra.Schur(Matrix(ps.R), Matrix(ps.Q), ps.eigenvalues)
+    IterativeRefinement.rfeigen(A, S, idxλ, DT, maxiter; kwargs...)
+end
+
+function IR.rfeigen(A::AbstractMatrix{T}, ps::ArnoldiMethod.PartialSchur{TQ,TR}, idxλ,
+                    DT = widen(real(T)), maxiter=5; kwargs...
+) where {T, TQ, TR <: AbstractMatrix{TRe}} where {TRe <: Real}
+    Sr = LinearAlgebra.Schur(Matrix(ps.R), Matrix(ps.Q), ps.eigenvalues)
+    S = LinearAlgebra.Schur{complex(TRe)}(Sr)
+    IterativeRefinement.rfeigen(A, S, idxλ, DT, maxiter; kwargs...)
+end
+
+end

src/eigen.jl

@@ -39,7 +39,7 @@ function rfeigen(A::AbstractMatrix{T},
                  factor = lu,
                  scale = true,
                  verbose = false
-                 ) where {T,Tλ,Tx}
+                 ) where {T,Tλ <: Union{Real,Complex},Tx}
     Tr = promote_type(promote_type(Tx,DT),Tλ)
     res = _rfeigen(A, x, λ, Tr, factor, maxiter, tol, scale, verbose)
     return res
@@ -59,7 +59,7 @@ function rfeigen(A::AbstractMatrix{T},
                  factor = lu,
                  scale = true,
                  verbose = false
-                 ) where {T,Tλ}
+                 ) where {T,Tλ <: Union{Real,Complex}}
     # CHECKME: is this condition adequate?
     if issymmetric(A) && (Tλ <: Real)
         Tx = Tλ
@@ -148,20 +148,25 @@ end
 """
     rfeigen(A, S::Schur, idxλ, DT, maxiter=5) -> vals, vecs, status
 
-Improves the precision of a cluster of eigenvalues of matrix `A`
-via multi-precision iterative refinement, using more-precise real type `DT`.
+Improves the precision of a cluster of eigenvalues of square matrix `A`
+via multi-precision iterative refinement, using more-precise real type `DT`,
+using a pre-computed Schur decomposition `S`.
 Returns improved estimates of eigenvalues and vectors generating
 the corresponding invariant subspace.
 
 This method works on the set of eigenvalues in `S.values` indexed by `idxλ`.
 It is designed to handle (nearly) defective cases, but will fail if
 the matrix is extremely non-normal or the initial estimates are poor.
-Note that `S` must be a true Schur decomposition, not a "real Schur".
+
+If `S` is a quasi-triangular "real Schur", it will be converted to
+a complex upper-triangular Schur decomposition. `S` may be a partial
+decomposition (i.e. fewer than `size(A,1)` vectors).
 """
 function rfeigen(A::AbstractMatrix{T}, S::Schur{TS}, idxλ,
                  DT = widen(real(T)), maxiter=5;
                  tol = 1, verbose = false) where {T, TS <: Complex}
     n = size(A,1)
+    ns = size(S.T,1)
     m = length(idxλ)
     λ = [S.values[idxλ[i]] for i in 1:m]
     Tw = promote_type(T,eltype(λ))
@@ -174,8 +179,8 @@ function rfeigen(A::AbstractMatrix{T}, S::Schur{TS}, idxλ,
     # M is an upper-triangular matrix of mixing coefficients
 
     # Most of the work is in the space of Schur vectors
-    Z = zeros(Tw, n, m)
-    z = zeros(Tw, n)
+    Z = zeros(Tw, ns, m)
+    z = zeros(Tw, ns)
     idxz = Vector{Int}(undef, m)
 
     k = idxλ[1]
@@ -195,7 +200,9 @@ function rfeigen(A::AbstractMatrix{T}, S::Schur{TS}, idxλ,
     for l=2:m
         kp = k
         k = idxλ[l]
-        @assert k > kp
+        if k <= kp
+            throw(ArgumentError("idxλ must be an increasing sequence"))
+        end
         x0 = (S.T[1:k-1,1:k-1] - λ[l]*I) \ S.T[1:k-1,k]
         X1 = (S.T[1:k-1,1:k-1] - λ[l]*I) \ Z[1:k-1,1:l-1]
         z[1:k-1] .= -x0
@@ -329,3 +336,9 @@ function rfeigen(A::AbstractMatrix{T}, S::Schur{TS}, idxλ,
     λnew = λbar .+ dλ
     λnew, Xnew, status
 end
+
+function rfeigen(A::AbstractMatrix{T}, S::Schur{TS}, idxλ, DT = widen(real(T)), maxiter=5;
+                 kwargs...) where {T, TS <: Real}
+    Sc = Schur{Complex{TS}}(S)
+    rfeigen(A, Sc, idxλ, DT, maxiter; kwargs...)
+end

test/eigen.jl

@@ -74,6 +74,7 @@ end
     tol = 20.0
     for n in [5,10,32]
         for k in [2,3]
+            # testset seems to screw with scope, so revert to ancient practice
             @label retry
             A = mkmat_multiple(n,k,etarget,dmin,T)
             # we need a true Schur here
@@ -99,3 +100,41 @@ end
         end
     end
 end
+
+using ArnoldiMethod: partialschur
+
+@testset "eigenvalue cluster $T" for T in (Float32, ComplexF32)
+    DT = widen(T)
+    RT = real(T)
+    etarget = T(2)
+    dmin = RT(1e3) * eps(RT)
+    maxit = 5
+    tol = 20.0
+    bcond = RT(100)
+    for n in [200]
+        for k in [2,3]
+            etargets = etarget * (T(1) .+ dmin * collect(0:k-1))
+            @label retry
+            A = mkmat_cond(n,etargets,bcond,T; fac22=0.1)
+            ps,_ = partialschur(A, nev=3*k)
+            Ad = convert.(DT,A)
+            ew = eigvals(Ad)
+            idx = findall(abs.(ps.eigenvalues .- etarget) .< 0.2)
+            # don't try if A is so nonnormal that initial estimates are bad
+            if length(idx) != k
+                @goto retry
+            end
+            e = eps(real(T))
+            if verbose
+                ewerrs = [minimum(abs.(ps.eigenvalues[j] .- ew)) for j in idx]
+                println("initial errors ", ewerrs / (e * n))
+            end
+            newew, newvecs = rfeigen(A, ps, idx, DT, maxit)
+            ewerrs = [minimum(abs.(newew[j] .- ew)) for j in 1:k]
+            if verbose
+                println("final errors ", ewerrs / (e * n))
+            end
+            @test maximum(ewerrs) / abs(etarget) < tol * e * n
+        end
+    end
+end

test/utils.jl

@@ -98,19 +98,20 @@ construct a matrix with a cluster of eigenvalues with specified condition.
 to have worse condition than the cluster of interest, which may be
 undesirable.)
 """
-function mkmat_cond(n, targets, cond, ::Type{T}; lbdiag=false, verbose=false) where T
+function mkmat_cond(n, targets, cond, ::Type{T};
+                    lbdiag=false, verbose=false, fac22=1) where T
     if (cond < 1)
         throw(ArgumentError("condition cannot be < 1"))
     end
     k = length(targets)
     verbose && println("Matrix{$T} N=$n $k-cluster w/ cond $cond")
     Tw = (T <: Real) ? Float64 : ComplexF64
+    roffset = (T <: Real) ? 1/T(2) : (1+im)/T(2)
     A11 = diagm(0=>Tw.(targets)) + triu(rand(Tw,k,k),1)
-    ews = rand(n-k)
     if lbdiag
-        A22 = diagm(0=>rand(Tw,n-k))
+        A22 = fac22 * diagm(0=>(rand(Tw,n-k) .- roffset))
     else
-        A22 = triu(rand(Tw,n-k,n-k))
+        A22 = fac22 * triu((rand(Tw,n-k,n-k) .- roffset))
     end
     R = rand(Tw,k,n-k)
     condr = sqrt(cond^2 - 1.0)