Add more coverage to general eigenvalue problem

Author: andreasnoack

src/eigenGeneral.jl

@@ -1,17 +1,8 @@
 using Printf
 using LinearAlgebra
-using LinearAlgebra: Givens, Rotation
-
-# Auxiliary
-function adiagmax(A::StridedMatrix)
-    adm = zero(typeof(real(A[1])))
-    @inbounds begin
-        for i in size(A, 1)
-            adm = max(adm, abs(A[i, i]))
-        end
-    end
-    return adm
-end
+using LinearAlgebra: Givens, Rotation, givens
+
+import Base: \
 
 # Hessenberg Matrix
 struct HessenbergMatrix{T,S<:StridedMatrix} <: AbstractMatrix{T}
@@ -30,19 +21,23 @@ function LinearAlgebra.ldiv!(H::HessenbergMatrix, B::AbstractVecOrMat)
     n = size(H, 1)
     Hd = H.data
     for i = 1:n-1
-        G, _ = givens!(Hd, i, i + 1, i)
-        lmul!(G, view(Hd, 1:n, i+1:n))
+        G, _ = givens(Hd, i, i + 1, i)
+        lmul!(G, view(Hd, 1:n, i:n))
         lmul!(G, B)
     end
-    ldiv!(Triangular(Hd, :U), B)
+    ldiv!(UpperTriangular(Hd), B)
 end
+(\)(H::HessenbergMatrix, B::AbstractVecOrMat) = ldiv!(copy(H), copy(B))
 
 # Hessenberg factorization
 struct HessenbergFactorization{T,S<:StridedMatrix,U} <: Factorization{T}
     data::S
     τ::Vector{U}
 end
 
+Base.copy(HF::HessenbergFactorization{T,S,U}) where {T,S,U} =
+    HessenbergFactorization{T,S,U}(copy(HF.data), copy(HF.τ))
+
 function _hessenberg!(A::StridedMatrix{T}) where {T}
     n = LinearAlgebra.checksquare(A)
     τ = Vector{Householder{T}}(undef, n - 1)
@@ -60,6 +55,14 @@ LinearAlgebra.hessenberg!(A::StridedMatrix) = _hessenberg!(A)
 
 Base.size(H::HessenbergFactorization, args...) = size(H.data, args...)
 
+function Base.getproperty(F::HessenbergFactorization, s::Symbol)
+    if s === :H
+        return HessenbergMatrix{eltype(F.data),typeof(F.data)}(F.data)
+    else
+        return getfield(F, s)
+    end
+end
+
 # Schur
 struct Schur{T,S<:StridedMatrix} <: Factorization{T}
     data::S
@@ -74,19 +77,12 @@ function Base.getproperty(F::Schur, s::Symbol)
     end
 end
 
-function wilkinson(Hmm, t, d)
-    λ1 = (t + sqrt(t * t - 4d)) / 2
-    λ2 = (t - sqrt(t * t - 4d)) / 2
-    return ifelse(abs(Hmm - λ1) < abs(Hmm - λ2), λ1, λ2)
-end
-
 # We currently absorb extra unsupported keywords in kwargs. These could e.g. be scale and permute. Do we want to check that these are false?
 function _schur!(
     H::HessenbergFactorization{T};
     tol = eps(real(T)),
     shiftmethod = :Francis,
     maxiter = 30 * size(H, 1),
-    kwargs...,
 ) where {T}
 
     n = size(H, 1)
@@ -176,6 +172,8 @@ function _schur!(
     return Schur{T,typeof(HH)}(HH, τ)
 end
 _schur!(A::StridedMatrix; kwargs...) = _schur!(_hessenberg!(A); kwargs...)
+
+# FIXME! Move this method to piracy extension
 LinearAlgebra.schur!(A::StridedMatrix; kwargs...) = _schur!(A; kwargs...)
 
 function singleShiftQR!(

test/eigengeneral.jl

@@ -14,6 +14,17 @@ using Test, GenericLinearAlgebra, LinearAlgebra
         @test sort(vGLA, by = cplxord) ≈ sort(vLAPACK, by = cplxord)
         @test sort(vGLA, by = cplxord) ≈ sort(complex(eltype(A)).(vBig), by = cplxord)
         @test issorted(vBig, by = cplxord)
+
+        if T <: Complex
+            @testset "Rayleigh shifts" begin
+                @test sort(
+                    GenericLinearAlgebra._eigvals!(
+                        GenericLinearAlgebra._schur!(copy(A), shiftmethod = :Rayleigh),
+                    ),
+                    by = t -> (real(t), imag(t)),
+                ) ≈ sort(eigvals(A), by = t -> (real(t), imag(t)))
+            end
+        end
     end
 
     @testset "make sure that solver doesn't hang" begin
@@ -235,6 +246,9 @@ using Test, GenericLinearAlgebra, LinearAlgebra
             A[:, (i+1):end] = A[:, (i+1):end] * HM'
         end
         @test tril(A, -2) ≈ zeros(n, n) atol = 1e-14
-    end
 
+        @test eigvals(HF.H) ≈ eigvals(A)
+        @test eigvals(HF.H) ≈ eigvals!(copy(HF))
+        @test HF.H \ ones(n) ≈ Matrix(HF.H) \ ones(n)
+    end
 end

test/lapack.jl

@@ -69,34 +69,41 @@ using GenericLinearAlgebra.LAPACK2
         @test !(_vals ≈ sort(eigvals(T)))
     end
 
-    @testset "syevd: eltype=$eltype, uplo=$uplo" for eltype in (Float32, Float64, ComplexF32, ComplexF64), uplo in ('U', 'L')
+    @testset "syevd: eltype=$eltype, uplo=$uplo" for eltype in (
+            Float32,
+            Float64,
+            ComplexF32,
+            ComplexF64,
+        ),
+        uplo in ('U', 'L')
+
         A = randn(eltype, n, n)
         A = A + A'
         if eltype <: Real
             vals, vecs = LAPACK2.syevd!('V', uplo, copy(A))
         else
             vals, vecs = LAPACK2.heevd!('V', uplo, copy(A))
         end
-        @test diag(vecs'*A*vecs) ≈ eigvals(A)
+        @test diag(vecs' * A * vecs) ≈ eigvals(A)
     end
 
-    @testset "tgevc: eltype=$eltype, side=$side, howmny=$howmny" for eltype in (Float32, Float64), side in ('L', 'R', 'B'), howmny in ('A', #='B', =#'S')
+    @testset "tgevc: eltype=$eltype, side=$side, howmny=$howmny" for eltype in
+                                                                     (Float32, Float64),
+        side in ('L', 'R', 'B'),
+        howmny in ('A', 'S')
+ #='B', =#
         select = ones(Int, n)
         S, P = triu(randn(eltype, n, n)), triu(randn(eltype, n, n))
-        VL, VR, m = LAPACK2.tgevc!(
-            side,
-            howmny,
-            select,
-            copy(S),
-            copy(P),
-        )
+        VL, VR, m = LAPACK2.tgevc!(side, howmny, select, copy(S), copy(P))
         if side ∈ ('R', 'B')
-            w = diag(S*VR) ./ diag(P*VR)
-            @test S*VR ≈ P*VR*Diagonal(w) rtol=sqrt(eps(eltype)) atol=sqrt(eps(eltype))
+            w = diag(S * VR) ./ diag(P * VR)
+            @test S * VR ≈ P * VR * Diagonal(w) rtol = sqrt(eps(eltype)) atol =
+                sqrt(eps(eltype))
         end
         if side ∈ ('L', 'B')
-            w = w = diag(VL'*S) ./ diag(VL'*P)
-            @test VL'*S ≈ Diagonal(w)*VL'*P rtol=sqrt(eps(eltype)) atol=sqrt(eps(eltype))
+            w = w = diag(VL' * S) ./ diag(VL' * P)
+            @test VL' * S ≈ Diagonal(w) * VL' * P rtol = sqrt(eps(eltype)) atol =
+                sqrt(eps(eltype))
         end
     end
 end