Added hseqr! LAPACK function

stdlib/LinearAlgebra/docs/src/index.md

@@ -745,6 +745,7 @@ LinearAlgebra.LAPACK.trexc!
 LinearAlgebra.LAPACK.trsen!
 LinearAlgebra.LAPACK.tgsen!
 LinearAlgebra.LAPACK.trsyl!
+LinearAlgebra.LAPACK.hseqr!
 ```
 
 ```@meta

stdlib/LinearAlgebra/src/lapack.jl

@@ -5741,6 +5741,104 @@ for (ormhr, elty) in
     end
 end
 
+for (hseqr, elty) in
+    ((:zhseqr_,:ComplexF64),
+     (:chseqr_,:ComplexF32))
+    @eval begin
+        # *     .. Scalar Arguments ..
+        #       CHARACTER          JOB, COMPZ
+        #       INTEGER            N, ILO, IHI, LWORK, LDH, LDZ, INFO
+        # *     ..
+        # *     .. Array Arguments ..
+        #       COMPLEX*16         H( LDH, * ), Z( LDZ, * ), WORK( * )
+        function hseqr!(job::AbstractChar, compz::AbstractChar, ilo::Integer, ihi::Integer,
+                        H::AbstractMatrix{$elty}, Z::AbstractMatrix{$elty})
+            require_one_based_indexing(H, Z)
+            chkstride1(H)
+            n = checksquare(H)
+            checksquare(Z) == n || throw(DimensionMismatch())
+            ldh = max(1, stride(H, 2))
+            ldz = max(1, stride(Z, 2))
+            w = similar(H, $elty, n)
+            work = Vector{$elty}(undef, 1)
+            lwork = BlasInt(-1)
+            info = Ref{BlasInt}()
+            for i = 1:2  # first call returns lwork as work[1]
+                ccall((@blasfunc($hseqr), libblastrampoline), Cvoid,
+                    (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
+                    Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
+                    Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
+                    Ptr{BlasInt}),
+                    job, compz, n, ilo, ihi,
+                    H, ldh, w, Z, ldz, work,
+                    lwork, info)
+                chklapackerror(info[])
+                if i == 1
+                    lwork = BlasInt(real(work[1]))
+                    resize!(work, lwork)
+                end
+            end
+            H, Z, w
+        end
+    end
+end
+
+for (hseqr, elty) in
+    ((:dhseqr_,:Float64),
+     (:shseqr_,:Float32))
+    @eval begin
+        # *     .. Scalar Arguments ..
+        #       CHARACTER          JOB, COMPZ
+        #       INTEGER            N, ILO, IHI, LWORK, LDH, LDZ, INFO
+        # *     ..
+        # *     .. Array Arguments ..
+        #       COMPLEX*16         H( LDH, * ), Z( LDZ, * ), WORK( * )
+        function hseqr!(job::AbstractChar, compz::AbstractChar, ilo::Integer, ihi::Integer,
+                        H::AbstractMatrix{$elty}, Z::AbstractMatrix{$elty})
+            require_one_based_indexing(H, Z)
+            chkstride1(H)
+            n = checksquare(H)
+            checksquare(Z) == n || throw(DimensionMismatch())
+            ldh = max(1, stride(H, 2))
+            ldz = max(1, stride(Z, 2))
+            wr = similar(H, $elty, n)
+            wi = similar(H, $elty, n)
+            work = Vector{$elty}(undef, 1)
+            lwork = BlasInt(-1)
+            info = Ref{BlasInt}()
+            for i = 1:2  # first call returns lwork as work[1]
+                ccall((@blasfunc($hseqr), libblastrampoline), Cvoid,
+                    (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
+                    Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ptr{$elty},
+                    Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
+                    Ptr{BlasInt}),
+                    job, compz, n, ilo, ihi,
+                    H, ldh, wr, wi, Z, ldz, work,
+                    lwork, info)
+                chklapackerror(info[])
+                if i == 1
+                    lwork = BlasInt(real(work[1]))
+                    resize!(work, lwork)
+                end
+            end
+            H, Z, complex.(wr, wi)
+        end
+    end
+end
+hseqr!(H::StridedMatrix{T}, Z::StridedMatrix{T}) where {T<:BlasFloat} = hseqr!('S', 'V', 1, size(H, 1), H, Z)
+hseqr!(H::StridedMatrix{T}) where {T<:BlasFloat} = hseqr!('S', 'I', 1, size(H, 1), H, similar(H))
+
+"""
+    hseqr!(job, compz, ilo, ihi, H, Z) -> (H, Z, w)
+
+Computes all eigenvalues and (optionally) the Schur factorization of a matrix
+reduced to Hessenberg form. If `H` is balanced with `gebal!`
+then `ilo` and `ihi` are the outputs of `gebal!`. Otherwise they should be
+`ilo = 1` and `ihi = size(H,2)`. `tau` contains the elementary reflectors of
+the factorization.
+"""
+hseqr!(job::AbstractChar, compz::AbstractChar, ilo::Integer, ihi::Integer, H::AbstractMatrix, Z::AbstractMatrix)
+
 for (hetrd, elty) in
     ((:dsytrd_,Float64),
      (:ssytrd_,Float32),

stdlib/LinearAlgebra/src/schur.jl

@@ -102,6 +102,8 @@ julia> A
 """
 schur!(A::StridedMatrix{<:BlasFloat}) = Schur(LinearAlgebra.LAPACK.gees!('V', A)...)
 
+schur!(A::UpperHessenberg{T}) where {T<:BlasFloat} = Schur(LinearAlgebra.LAPACK.hseqr!(parent(A))...)
+
 """
     schur(A) -> F::Schur
 
@@ -153,6 +155,7 @@ true
 ```
 """
 schur(A::AbstractMatrix{T}) where {T} = schur!(copy_similar(A, eigtype(T)))
+schur(A::UpperHessenberg{T}) where {T} = schur!(copy_similar(A, eigtype(T)))
 function schur(A::RealHermSymComplexHerm)
     F = eigen(A; sortby=nothing)
     return Schur(typeof(F.vectors)(Diagonal(F.values)), F.vectors, F.values)

stdlib/LinearAlgebra/test/schur.jl

@@ -202,4 +202,20 @@ end
     @test A' ≈ C ≈ E
 end
 
+@testset "UpperHessenberg schur" begin
+    A = UpperHessenberg(rand(ComplexF64, 100, 100))
+    B = Array(A)
+    fact1 = schur(A)
+    fact2 = schur(B)
+    @test fact1.values ≈ fact2.values
+    @test fact1.Z * fact1.T * fact1.Z' ≈ B
+
+    A = UpperHessenberg(rand(Int32, 50, 50))
+    B = Array(A)
+    fact1 = schur(A)
+    fact2 = schur(B)
+    @test fact1.values ≈ fact2.values
+    @test fact1.Z * fact1.T * fact1.Z' ≈ B
+end
+
 end # module TestSchur