diff --git a/Project.toml b/Project.toml index 13a71f94..ea17e7d2 100644 --- a/Project.toml +++ b/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" diff --git a/src/lapack.jl b/src/lapack.jl index d84698dd..5cb8eb13 100644 --- a/src/lapack.jl +++ b/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)) diff --git a/src/symbanded/bandedeigen.jl b/src/symbanded/bandedeigen.jl index f9c12c80..b14a8e57 100644 --- a/src/symbanded/bandedeigen.jl +++ b/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) diff --git a/src/symbanded/symbanded.jl b/src/symbanded/symbanded.jl index f455afee..1d0bda14 100644 --- a/src/symbanded/symbanded.jl +++ 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")) diff --git a/src/symbanded/tridiagonalize.jl b/src/symbanded/tridiagonalize.jl index 8dd908a8..f748b390 100644 --- a/src/symbanded/tridiagonalize.jl +++ 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))) diff --git a/test/test_symbanded.jl b/test/test_symbanded.jl index ac245d75..26aee6a8 100644 --- a/test/test_symbanded.jl +++ b/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