From 2c335d14a622f879d41f49b39b436656f7023ff6 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 18 May 2023 14:45:59 +0530 Subject: [PATCH 1/6] in-place eigvals for complex hermitian BandedMatrix --- src/lapack.jl | 30 ++++++++++++++++++++++++++++++ src/symbanded/symbanded.jl | 18 +++--------------- src/symbanded/tridiagonalize.jl | 21 +++++++++++++++++++++ test/test_symbanded.jl | 29 ++++++++++++++++++++--------- 4 files changed, 74 insertions(+), 24 deletions(-) diff --git a/src/lapack.jl b/src/lapack.jl index d84698dd..83ef2e96 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -115,6 +115,36 @@ for (fname, elty) in ((:dsbtrd_,:Float64), end end +for (fname, elty) in ((:zhbtrd_,:ComplexF64), + (:chbtrd_,:ComplexF32)) + @eval begin + 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{$elty}, Ptr{$elty}, 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/symbanded.jl b/src/symbanded/symbanded.jl index f455afee..5cde39fb 100644 --- a/src/symbanded/symbanded.jl +++ b/src/symbanded/symbanded.jl @@ -117,25 +117,13 @@ 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<:Union{Float32, Float64} = 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)) +eigvals!(A::HermOrSym{T,<:BandedMatrix{T}}) where T <: Real = eigvals!(tridiagonalize!(A)) +eigvals!(A::Hermitian{T,<:BandedMatrix{T}}) where T <: Complex = 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) 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..0703973a 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 + for T in [Float32, Float64] + A = brand(T, 10, 10, 2, 4) + std = eigvals(Symmetric(Matrix(A))) + @test eigvals(Symmetric(A)) ≈ std + @test eigvals!(copy(Symmetric(A))) ≈ std + @test eigvals(Hermitian(A)) ≈ std + @test eigvals!(copy(Hermitian(A))) ≈ std + @test eigvals(Hermitian(big.(A))) ≈ std + + @test eigvals(Symmetric(A, :L)) ≈ eigvals(Symmetric(Matrix(A), :L)) + @test eigvals(Hermitian(A, :L)) ≈ eigvals(Hermitian(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,9 +69,13 @@ 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 + for T in [ComplexF64, ComplexF32] + A = brand(T, 20, 20, 4, 2) + std = eigvals(Hermitian(Matrix(A), :L)) + @test eigvals(Hermitian(A, :L)) ≈ std + @test eigvals(Hermitian(big.(A), :L)) ≈ std + @test eigvals!(copy(Hermitian(A, :L))) ≈ std + end A = Symmetric(brand(Float64, 100, 100, 2, 4)) F = eigen(A) From 6e3d0427ecdaf95d0d22232bb787547fa9ca9530 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 18 May 2023 14:46:23 +0530 Subject: [PATCH 2/6] bump version to 0.17.24 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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" From aa26da48dc7db26cab207ae4cb088ec3de6907df Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 18 May 2023 16:06:50 +0530 Subject: [PATCH 3/6] fix pointer types --- src/lapack.jl | 35 ++++++++++++++++++++++++++++------- 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/src/lapack.jl b/src/lapack.jl index 83ef2e96..5cb8eb13 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -118,7 +118,7 @@ end for (fname, elty) in ((:zhbtrd_,:ComplexF64), (:chbtrd_,:ComplexF32)) @eval begin - Relty = real($elty) + 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}, @@ -133,12 +133,33 @@ for (fname, elty) in ((:zhbtrd_,:ComplexF64), 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{$elty}, Ptr{$elty}, 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) + ( + 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 From b0fe27131710f8aa7d6a2343a8a458a89a191def Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 18 May 2023 17:28:34 +0530 Subject: [PATCH 4/6] complex hermitian inplace eigen --- src/symbanded/bandedeigen.jl | 15 ++++++++++++++- test/test_symbanded.jl | 32 ++++++++++++++++++++------------ 2 files changed, 34 insertions(+), 13 deletions(-) diff --git a/src/symbanded/bandedeigen.jl b/src/symbanded/bandedeigen.jl index f9c12c80..f6d97b42 100644 --- a/src/symbanded/bandedeigen.jl +++ b/src/symbanded/bandedeigen.jl @@ -50,7 +50,7 @@ function eigen(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{ eigen!(AA, copy(B)) end -function eigen!(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real +function eigen!(A::HermOrSym{T,<:BandedMatrix{T}}) where T <: Real N = size(A, 1) KD = bandwidth(A) D = Vector{T}(undef, N) @@ -63,6 +63,19 @@ function eigen!(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real Eigen(Λ, BandedEigenvectors(G, Q, similar(Q, size(Q,1)))) end +function eigen!(A::Hermitian{T,<:BandedMatrix{T}}) 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)) + 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/test/test_symbanded.jl b/test/test_symbanded.jl index 0703973a..0f531ff7 100644 --- a/test/test_symbanded.jl +++ b/test/test_symbanded.jl @@ -47,17 +47,14 @@ import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedCol # (generalized) eigen & eigvals Random.seed!(0) - for T in [Float32, Float64] + @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!(copy(Symmetric(A))) ≈ std - @test eigvals(Hermitian(A)) ≈ std @test eigvals!(copy(Hermitian(A))) ≈ std - @test eigvals(Hermitian(big.(A))) ≈ std @test eigvals(Symmetric(A, :L)) ≈ eigvals(Symmetric(Matrix(A), :L)) - @test eigvals(Hermitian(A, :L)) ≈ eigvals(Hermitian(Matrix(A), :L)) end A = brand(ComplexF64, 20, 20, 4, 0) @@ -69,14 +66,6 @@ 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] - for T in [ComplexF64, ComplexF32] - A = brand(T, 20, 20, 4, 2) - std = eigvals(Hermitian(Matrix(A), :L)) - @test eigvals(Hermitian(A, :L)) ≈ std - @test eigvals(Hermitian(big.(A), :L)) ≈ std - @test eigvals!(copy(Hermitian(A, :L))) ≈ std - end - A = Symmetric(brand(Float64, 100, 100, 2, 4)) F = eigen(A) Λ, Q = F @@ -188,6 +177,25 @@ 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(big.(A), uplo)) ≈ std + @test eigvals!(copy(Hermitian(A, uplo))) ≈ std + 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!(A) + λ2, v2 = eigen!(copy(MA)) + @test λ1 ≈ λ2 + @test v1' * MA * v1 ≈ Diagonal(λ1) + end + end end @testset "LDLᵀ" begin From 4391c515e695ed6e85aefbae8c31d72d3381d2e4 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 18 May 2023 23:59:36 +0530 Subject: [PATCH 5/6] accept ranges in eigen --- src/symbanded/bandedeigen.jl | 44 ++++++++++++++++++++++++++++++++---- src/symbanded/symbanded.jl | 28 ----------------------- test/test_symbanded.jl | 39 +++++++++++++++++++++++++++++--- 3 files changed, 75 insertions(+), 36 deletions(-) diff --git a/src/symbanded/bandedeigen.jl b/src/symbanded/bandedeigen.jl index f6d97b42..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::HermOrSym{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,11 +93,11 @@ function eigen!(A::HermOrSym{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}}) where T <: Complex +function eigen!(A::Hermitian{T,<:BandedMatrix{T}}, args...) where T <: Complex N = size(A, 1) KD = bandwidth(A) D = Vector{real(T)}(undef, N) @@ -72,7 +106,7 @@ function eigen!(A::Hermitian{T,<:BandedMatrix{T}}) where T <: Complex WORK = Vector{T}(undef, N) AB = hermbandeddata(A) hbtrd!('V', A.uplo, N, KD, AB, D, E, Q, WORK) - Λ, W = eigen(SymTridiagonal(D, E)) + Λ, W = eigen!(SymTridiagonal(D, E), args...) Eigen(Λ, Q * W) end diff --git a/src/symbanded/symbanded.jl b/src/symbanded/symbanded.jl index 5cde39fb..1d0bda14 100644 --- a/src/symbanded/symbanded.jl +++ b/src/symbanded/symbanded.jl @@ -114,34 +114,6 @@ function materialize!(M::BlasMatMulVecAdd{<:HermitianLayout{<:BandedColumnMajor} _banded_hbmv!(symmetricuplo(A), α, A, x, β, y) end - -## eigvals routine - -eigvals(A::Symmetric{T,<:BandedMatrix{T}}) where T<:Union{Float32, Float64} = 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::HermOrSym{T,<:BandedMatrix{T}}) where T <: Real = eigvals!(tridiagonalize!(A)) -eigvals!(A::Hermitian{T,<:BandedMatrix{T}}) where T <: Complex = 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/test/test_symbanded.jl b/test/test_symbanded.jl index 0f531ff7..a95af942 100644 --- a/test/test_symbanded.jl +++ b/test/test_symbanded.jl @@ -51,8 +51,11 @@ import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedCol 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(Hermitian(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 @@ -66,13 +69,23 @@ 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] - 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(Λ) FD = convert(Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}, F) @test FD.vectors'Matrix(A)*FD.vectors ≈ Diagonal(F.values) + 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(Λ) + MQ = Matrix(Q) @test Q[axes(Q,1),1:3] ≈ MQ[axes(Q,1),1:3] @test Q[:,1:3] ≈ MQ[:,1:3] @@ -183,17 +196,37 @@ end 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!(A) + λ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 From 11da89caecc04b5a7aed6435a9a555b882d69254 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Fri, 19 May 2023 00:19:20 +0530 Subject: [PATCH 6/6] reorder tests --- test/test_symbanded.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/test_symbanded.jl b/test/test_symbanded.jl index a95af942..26aee6a8 100644 --- a/test/test_symbanded.jl +++ b/test/test_symbanded.jl @@ -76,6 +76,12 @@ import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedCol FD = convert(Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}, F) @test FD.vectors'Matrix(A)*FD.vectors ≈ Diagonal(F.values) + MQ = Matrix(Q) + @test Q[axes(Q,1),1:3] ≈ MQ[axes(Q,1),1:3] + @test Q[:,1:3] ≈ MQ[:,1:3] + @test Q[:,3] ≈ MQ[:,3] + @test Q[1,2] ≈ MQ[1,2] + F = eigen(A, 2:4) Λ, Q = F QM = Matrix(Q) @@ -86,12 +92,6 @@ import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedCol QM = Matrix(Q) @test QM' * (Matrix(A)*QM) ≈ Diagonal(Λ) - MQ = Matrix(Q) - @test Q[axes(Q,1),1:3] ≈ MQ[axes(Q,1),1:3] - @test Q[:,1:3] ≈ MQ[:,1:3] - @test Q[:,3] ≈ MQ[:,3] - @test Q[1,2] ≈ MQ[1,2] - function An(::Type{T}, N::Int) where {T} A = Symmetric(BandedMatrix(Zeros{T}(N,N), (0, 2))) for n = 0:N-1