In-place eigvals for complex Hermitian BandedMatrix

Project.toml

@@ -1,6 +1,6 @@
 name = "BandedMatrices"
 uuid = "aae01518-5342-5314-be14-df237901396f"
-version = "0.17.23"
+version = "0.17.24"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/lapack.jl

@@ -115,6 +115,57 @@ for (fname, elty) in ((:dsbtrd_,:Float64),
     end
 end
 
+for (fname, elty) in ((:zhbtrd_,:ComplexF64),
+                      (:chbtrd_,:ComplexF32))
+    @eval begin
+        local Relty = real($elty)
+        function hbtrd!(vect::Char, uplo::Char,
+                        m::Int, k::Int, A::AbstractMatrix{$elty},
+                        d::AbstractVector{Relty}, e::AbstractVector{Relty}, Q::AbstractMatrix{$elty},
+                        work::AbstractVector{$elty})
+            require_one_based_indexing(A)
+            chkstride1(A)
+            chkuplo(uplo)
+            chkvect(vect)
+            info = Ref{BlasInt}()
+            n    = size(A,2)
+            n ≠ m && throw(ArgumentError("Matrix must be square"))
+            size(A,1) < k+1 && throw(ArgumentError("Not enough bands"))
+            info  = Ref{BlasInt}()
+            ccall((@blasfunc($fname), liblapack), Nothing,
+                (
+                    Ref{UInt8},
+                    Ref{UInt8},
+                    Ref{BlasInt},
+                    Ref{BlasInt},
+                    Ptr{$elty},
+                    Ref{BlasInt},
+                    Ptr{Relty},
+                    Ptr{Relty},
+                    Ptr{$elty},
+                    Ref{BlasInt},
+                    Ptr{$elty},
+                    Ptr{BlasInt},
+                ),
+                 vect,
+                 uplo,
+                 n,
+                 k,
+                 A,
+                 max(1,stride(A,2)),
+                 d,
+                 e,
+                 Q,
+                 max(1,stride(Q,2)),
+                 work,
+                 info
+            )
+            LAPACK.chklapackerror(info[])
+            d, e, Q
+        end
+    end
+end
+
 ## Bidiagonalization
 for (fname, elty) in ((:dgbbrd_,:Float64),
                       (:sgbbrd_,:Float32))

src/symbanded/bandedeigen.jl

@@ -1,3 +1,34 @@
+## eigvals routine
+
+# the symmetric case uses lapack throughout
+eigvals(A::Symmetric{T,<:BandedMatrix{T}}) where T<:Union{Float32, Float64} = eigvals!(copy(A))
+eigvals(A::Symmetric{T,<:BandedMatrix{T}}, irange::UnitRange) where T<:Union{Float32, Float64} = eigvals!(copy(A), irange)
+eigvals(A::Symmetric{T,<:BandedMatrix{T}}, vl::Real, vu::Real) where T<:Union{Float32, Float64} = eigvals!(copy(A), vl, vu)
+
+eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}) = eigvals!(tridiagonalize(A))
+eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, irange::UnitRange) = eigvals!(tridiagonalize(A), irange)
+eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, vl::Real, vu::Real) = eigvals!(tridiagonalize(A), vl, vu)
+
+eigvals!(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, args...) = eigvals!(tridiagonalize!(A), args...)
+
+function eigvals!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T<:Real
+    n = size(A, 1)
+    @assert n == size(B, 1)
+    @assert A.uplo == B.uplo
+    # compute split-Cholesky factorization of B.
+    kb = bandwidth(B)
+    B_data = symbandeddata(B)
+    pbstf!(B.uplo, n, kb, B_data)
+    # convert to a regular symmetric eigenvalue problem.
+    ka = bandwidth(A)
+    A_data = symbandeddata(A)
+    X = Array{T}(undef,0,0)
+    work = Vector{T}(undef,2n)
+    sbgst!('N', A.uplo, n, ka, kb, A_data, B_data, X, work)
+    # compute eigenvalues of symmetric eigenvalue problem.
+    eigvals!(A)
+end
+
 # V = G Q
 struct BandedEigenvectors{T} <: AbstractMatrix{T}
     G::Vector{Givens{T}}
@@ -44,13 +75,16 @@ convert(::Type{GeneralizedEigen{T, T, Matrix{T}, Vector{T}}}, F::GeneralizedEige
 compress(F::Eigen{T, T, BandedEigenvectors{T}, Vector{T}}) where T = convert(Eigen{T, T, Matrix{T}, Vector{T}}, F)
 compress(F::GeneralizedEigen{T, T, BandedGeneralizedEigenvectors{T}, Vector{T}}) where T = convert(GeneralizedEigen{T, T, Matrix{T}, Vector{T}}, F)
 
-eigen(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real = eigen!(copy(A))
+eigen(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}) = eigen!(copy(A))
+eigen(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, irange::UnitRange) = eigen!(copy(A), irange)
+eigen(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, vl::Real, vu::Real) = eigen!(copy(A), vl, vu)
+
 function eigen(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
     AA = _copy_bandedsym(A, B)
     eigen!(AA, copy(B))
 end
 
-function eigen!(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
+function eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, args...) where T <: Real
     N = size(A, 1)
     KD = bandwidth(A)
     D = Vector{T}(undef, N)
@@ -59,10 +93,23 @@ function eigen!(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
     WORK = Vector{T}(undef, N)
     AB = symbandeddata(A)
     sbtrd!('V', A.uplo, N, KD, AB, D, E, G, WORK)
-    Λ, Q = eigen(SymTridiagonal(D, E))
+    Λ, Q = eigen!(SymTridiagonal(D, E), args...)
     Eigen(Λ, BandedEigenvectors(G, Q, similar(Q, size(Q,1))))
 end
 
+function eigen!(A::Hermitian{T,<:BandedMatrix{T}}, args...) where T <: Complex
+    N = size(A, 1)
+    KD = bandwidth(A)
+    D = Vector{real(T)}(undef, N)
+    E = Vector{real(T)}(undef, N-1)
+    Q = Matrix{T}(undef, N, N)
+    WORK = Vector{T}(undef, N)
+    AB = hermbandeddata(A)
+    hbtrd!('V', A.uplo, N, KD, AB, D, E, Q, WORK)
+    Λ, W = eigen!(SymTridiagonal(D, E), args...)
+    Eigen(Λ, Q * W)
+end
+
 function eigen!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
     isdiag(A) || isdiag(B) || symmetricuplo(A) == symmetricuplo(B) || throw(ArgumentError("uplo of matrices do not match"))
     S = splitcholesky!(B)

src/symbanded/symbanded.jl

@@ -114,46 +114,6 @@ function materialize!(M::BlasMatMulVecAdd{<:HermitianLayout{<:BandedColumnMajor}
     _banded_hbmv!(symmetricuplo(A), α, A, x, β, y)
 end
 
-
-## eigvals routine
-
-
-function _tridiagonalize!(A::AbstractMatrix{T}, ::SymmetricLayout{<:BandedColumnMajor}) where T
-    n=size(A, 1)
-    d = Vector{T}(undef,n)
-    e = Vector{T}(undef,n-1)
-    Q = Matrix{T}(undef,0,0)
-    work = Vector{T}(undef,n)
-    sbtrd!('N', symmetricuplo(A), size(A,1), bandwidth(A), symbandeddata(A), d, e, Q, work)
-    SymTridiagonal(d,e)
-end
-
-tridiagonalize!(A::AbstractMatrix) = _tridiagonalize!(A, MemoryLayout(typeof(A)))
-
-eigvals(A::Symmetric{T,<:BandedMatrix{T}}) where T<:Base.IEEEFloat = eigvals!(copy(A))
-eigvals(A::Symmetric{T,<:BandedMatrix{T}}) where T<:Real = eigvals!(tridiagonalize(A))
-eigvals(A::Hermitian{T,<:BandedMatrix{T}}) where T<:Real = eigvals!(tridiagonalize(A))
-eigvals(A::Hermitian{T,<:BandedMatrix{T}}) where T<:Complex = eigvals!(tridiagonalize(A))
-
-eigvals!(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real = eigvals!(tridiagonalize!(A))
-function eigvals!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T<:Real
-    n = size(A, 1)
-    @assert n == size(B, 1)
-    @assert A.uplo == B.uplo
-    # compute split-Cholesky factorization of B.
-    kb = bandwidth(B)
-    B_data = symbandeddata(B)
-    pbstf!(B.uplo, n, kb, B_data)
-    # convert to a regular symmetric eigenvalue problem.
-    ka = bandwidth(A)
-    A_data = symbandeddata(A)
-    X = Array{T}(undef,0,0)
-    work = Vector{T}(undef,2n)
-    sbgst!('N', A.uplo, n, ka, kb, A_data, B_data, X, work)
-    # compute eigenvalues of symmetric eigenvalue problem.
-    eigvals!(A)
-end
-
 function copyto!(A::Symmetric{<:Number,<:BandedMatrix}, B::Symmetric{<:Number,<:BandedMatrix})
     size(A) == size(B) || throw(ArgumentError("sizes of A and B must match"))
     bandwidth(A) >= bandwidth(B) || throw(ArgumentError("bandwidth of A must exceed that of B"))

src/symbanded/tridiagonalize.jl

@@ -84,3 +84,24 @@ function copybands(A::AbstractMatrix{T}, d::Integer) where T
     B
 end
 
+function _tridiagonalize!(A::AbstractMatrix{T}, ::SymmetricLayout{<:BandedColumnMajor}) where T<:Union{Float32,Float64}
+    n=size(A, 1)
+    d = Vector{T}(undef,n)
+    e = Vector{T}(undef,n-1)
+    Q = Matrix{T}(undef,0,0)
+    work = Vector{T}(undef,n)
+    sbtrd!('N', symmetricuplo(A), size(A,1), bandwidth(A), symbandeddata(A), d, e, Q, work)
+    SymTridiagonal(d,e)
+end
+
+function _tridiagonalize!(A::AbstractMatrix{T}, ::HermitianLayout{<:BandedColumnMajor}) where T<:Union{ComplexF32,ComplexF64}
+    n=size(A, 1)
+    d = Vector{real(T)}(undef,n)
+    e = Vector{real(T)}(undef,n-1)
+    Q = Matrix{T}(undef,0,0)
+    work = Vector{T}(undef,n)
+    hbtrd!('N', symmetricuplo(A), size(A,1), bandwidth(A), hermbandeddata(A), d, e, Q, work)
+    SymTridiagonal(d,e)
+end
+
+tridiagonalize!(A::AbstractMatrix) = _tridiagonalize!(A, MemoryLayout(typeof(A)))

test/test_symbanded.jl

@@ -47,13 +47,20 @@ import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedCol
     # (generalized) eigen & eigvals
     Random.seed!(0)
 
-    A = brand(Float64, 100, 100, 2, 4)
-    std = eigvals(Symmetric(Matrix(A)))
-    @test eigvals(Symmetric(A)) ≈ std
-    @test eigvals(Hermitian(A)) ≈ std
-    @test eigvals(Hermitian(big.(A))) ≈ std
+    @testset for T in (Float32, Float64)
+        A = brand(T, 10, 10, 2, 4)
+        std = eigvals(Symmetric(Matrix(A)))
+        @test eigvals(Symmetric(A)) ≈ std
+        @test eigvals(Symmetric(A), 2:4) ≈ std[2:4]
+        @test eigvals(Symmetric(A), 1, 2) ≈ std[1 .<= std .<= 2]
+        @test eigvals!(copy(Symmetric(A))) ≈ std
+        @test eigvals!(copy(Symmetric(A)), 2:4) ≈ std[2:4]
+        @test eigvals!(copy(Symmetric(A)), 1, 2) ≈ std[1 .<= std .<= 2]
+
+        @test eigvals(Symmetric(A, :L)) ≈ eigvals(Symmetric(Matrix(A), :L))
+    end
 
-    A = brand(ComplexF64, 100, 100, 4, 0)
+    A = brand(ComplexF64, 20, 20, 4, 0)
     @test Symmetric(A)[2:10,1:9] isa BandedMatrix
     @test Hermitian(A)[2:10,1:9] isa BandedMatrix
     @test isempty(Hermitian(A)[1:0,1:9])
@@ -62,11 +69,7 @@ import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedCol
     @test [Symmetric(A)[k,j] for k=2:10, j=1:9] == Symmetric(A)[2:10,1:9]
     @test [Hermitian(A)[k,j] for k=2:10, j=1:9] == Hermitian(A)[2:10,1:9]
 
-    std = eigvals(Hermitian(Matrix(A), :L))
-    @test eigvals(Hermitian(A, :L)) ≈ std
-    @test eigvals(Hermitian(big.(A), :L)) ≈ std
-
-    A = Symmetric(brand(Float64, 100, 100, 2, 4))
+    A = Symmetric(brand(Float64, 10, 10, 2, 4))
     F = eigen(A)
     Λ, Q = F
     @test Q'Matrix(A)*Q ≈ Diagonal(Λ)
@@ -79,6 +82,16 @@ import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedCol
     @test Q[:,3] ≈ MQ[:,3]
     @test Q[1,2] ≈ MQ[1,2]
 
+    F = eigen(A, 2:4)
+    Λ, Q = F
+    QM = Matrix(Q)
+    @test QM' * (Matrix(A)*QM) ≈ Diagonal(Λ)
+
+    F = eigen(A, 1, 2)
+    Λ, Q = F
+    QM = Matrix(Q)
+    @test QM' * (Matrix(A)*QM) ≈ Diagonal(Λ)
+
     function An(::Type{T}, N::Int) where {T}
         A = Symmetric(BandedMatrix(Zeros{T}(N,N), (0, 2)))
         for n = 0:N-1
@@ -177,6 +190,45 @@ end
     @test all(A*x .=== muladd!(one(T),A,x,zero(T),copy(x)) .===
                 materialize!(MulAdd(one(T),A,x,zero(T),copy(x))) .===
                 BLAS.hbmv!('L', 1, one(T), view(parent(A).data,3:4,:), x, zero(T), similar(x)))
+
+    @testset "eigen" begin
+        @testset for T in (Float32, Float64, ComplexF64, ComplexF32), uplo in (:L, :U)
+            A = brand(T, 20, 20, 4, 2)
+            std = eigvals(Hermitian(Matrix(A), uplo))
+            @test eigvals(Hermitian(A, uplo)) ≈ std
+            @test eigvals(Hermitian(A, uplo), 2:4) ≈ std[2:4]
+            @test eigvals(Hermitian(A, uplo), 1, 2) ≈ std[1 .<= std .<= 2]
+            @test eigvals(Hermitian(big.(A), uplo)) ≈ std
+            @test eigvals!(copy(Hermitian(A, uplo))) ≈ std
+            @test eigvals!(copy(Hermitian(A, uplo)), 2:4) ≈ std[2:4]
+            @test eigvals!(copy(Hermitian(A, uplo)), 1, 2) ≈ std[1 .<= std .<= 2]
+        end
+
+        @testset for T in (Float32, Float64, ComplexF32, ComplexF64), uplo in (:L, :U)
+            A = Hermitian(brand(T, 20, 20, 2, 4), uplo)
+            MA = Hermitian(Matrix(A), uplo)
+            λ1, v1 = eigen!(copy(A))
+            λ2, v2 = eigen!(copy(MA))
+            @test λ1 ≈ λ2
+            @test v1' * MA * v1 ≈ Diagonal(λ1)
+
+            λ3, v3 = eigen!(copy(A), 2:4)
+            @test λ3 ≈ λ1[2:4]
+            @test v3' * (MA * v3) ≈ Diagonal(λ3)
+
+            λ4, v4 = eigen(A, 2:4)
+            @test λ4 ≈ λ3
+            @test v4' * (MA * v4) ≈ Diagonal(λ4)
+
+            λ3, v3 = eigen!(copy(A), 1, 2)
+            @test λ3 ≈ λ1[1 .<= λ1 .<= 2]
+            @test v3' * (MA * v3) ≈ Diagonal(λ3)
+
+            λ4, v4 = eigen!(copy(A), 1, 2)
+            @test λ4 ≈ λ1[1 .<= λ1 .<= 2]
+            @test v4' * (MA * v4) ≈ Diagonal(λ4)
+        end
+    end
 end
 
 @testset "LDLᵀ" begin