diff --git a/Project.toml b/Project.toml index 2dba0f6..c6ef3bf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SkewLinearAlgebra" uuid = "5c889d49-8c60-4500-9d10-5d3a22e2f4b9" authors = ["smataigne and contributors"] -version = "0.1.0" +version = "0.1.1" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -13,6 +13,7 @@ LinearAlgebra = "1.6" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [targets] -test = ["Test","Random"] +test = ["Test","Random", "SparseArrays"] diff --git a/src/skeweigen.jl b/src/skeweigen.jl index 1a897ee..9bfa804 100644 --- a/src/skeweigen.jl +++ b/src/skeweigen.jl @@ -67,23 +67,21 @@ end end function getshift(ev::AbstractVector{T}) where T - if abs(ev[2]) ≈ abs(ev[1]) - return 0 - end return ev[2]^2 end -@views function implicitstep_novec(ev::AbstractVector{T} , n::Integer ) where T +@views function implicitstep_novec(ev::AbstractVector{T} , n::Integer, start::Integer) where T bulge = zero(T) shift = getshift(ev[n-1:n]) - activation = 0 - @inbounds(for i = 1:n-1 - α = (i > 1 ? ev[i-1] : zero(ev[i])) + tol = T(1) * eps(T) + @inbounds(for i = start:n-1 + α = (i > start ? ev[i-1] : zero(ev[i])) β = ev[i] γ = ev[i+1] x1 = - α * α - β * β + shift x2 = - α * bulge + β * γ - c, s = ((iszero(x1) && iszero(x2)) ? getgivens(α, bulge) : getgivens(x1, x2)) + c, s = (i > start ? getgivens(α, bulge) : getgivens(x1, x2)) + if i > 1 ev[i-1] = c * α + s * bulge end @@ -95,17 +93,13 @@ end ζ = ev[i+2] ev[i+2] *= c bulge = s * ζ - end - #Make it compulsory to initiate a bulge before stopping - if !iszero(bulge) - activation += 1 - else - if activation > 0 - return + if abs(bulge) < tol && abs(ev[i]) < tol + start = i + 1 + return start end end end) - return + return start end @views function skewtrieigvals!(A::SkewHermTridiagonal{T,V,Vim}) where {T<:Real,V<:AbstractVector{T},Vim<:Nothing} @@ -120,9 +114,10 @@ end tol = eps(T) * T(10) max_iter = 30 * n iter = 0 ; - N = n + mem = T(1); count_static = 0 #mem and count_static allow to detect the failure of Wilkinson shifts. + start = 1 #start remembers if a zero eigenvalue appeared in the middle of ev. while n > 2 && iter < max_iter - implicitstep_novec(ev, n - 1) + start = implicitstep_novec(ev, n - 1, start) while n > 2 && abs(ev[n-1]) <= tol values[n] = 0 n -= 1 @@ -131,6 +126,20 @@ end eigofblock(ev[n - 1], values[n-1:n] ) n -= 2 end + if start > n-2 + start = 1 + end + if n>1 && abs(mem - ev[n-1]) < T(0.0001) * abs(ev[n-1]) + count_static += 1 + if count_static > 4 + #Wilkinson shifts have failed, change strategy using LAPACK tridiagonal symmetric solver. + values[1:n] .= complex.(0, skewtrieigvals_backup!(SkewHermTridiagonal(ev[1:n-1]))) + return values + end + else + count_static = 0 + end + mem = (n>1 ? ev[n-1] : T(0)) iter += 1 end if n == 2 @@ -146,18 +155,18 @@ end end end -@views function implicitstep_vec!(ev::AbstractVector{T}, Qeven::AbstractMatrix{T}, Qodd::AbstractMatrix{T}, n::Integer, N::Integer) where T +@views function implicitstep_vec!(ev::AbstractVector{T}, Qeven::AbstractMatrix{T}, Qodd::AbstractMatrix{T}, n::Integer, N::Integer, start::Integer) where T bulge = zero(T) shift = getshift(ev[n-1:n]) - activation = 0 - @inbounds(for i=1:n-1 - α = (i > 1 ? ev[i-1] : zero(ev[i])) + tol = 10 * eps(T) + @inbounds(for i = start:n-1 + α = (i > start ? ev[i-1] : zero(ev[i])) β = ev[i] γ = ev[i+1] x1 = - α * α - β * β + shift x2 = - α * bulge + β * γ - c, s = ((iszero(x1) && iszero(x2)) ? getgivens(α,bulge) : getgivens(x1, x2)) + c, s = (i > start ? getgivens(α,bulge) : getgivens(x1, x2)) if i > 1 ev[i-1] = c * α + s * bulge end @@ -167,6 +176,9 @@ end ζ = ev[i+2] ev[i+2] *= c bulge = s * ζ + if abs(bulge) < tol && abs(ev[i]) < tol + start = i + 1 + end end Q = (isodd(i) ? Qodd : Qeven) k = div(i+1, 2) @@ -176,19 +188,12 @@ end Q[j, k] = c*σ + s*ω Q[j, k+1] = -s*σ + c*ω end - - if !iszero(bulge) - activation += 1 - else - if activation > 0 - return - end - end end) - return + return start end @views function skewtrieigen_merged!(A::SkewHermTridiagonal{T}) where {T<:Real} + Backup = copy(A) n = size(A, 1) values = complex(zeros(T, n)) vectors = similar(A, Complex{T}, n, n) @@ -202,13 +207,14 @@ end reducetozero(ev, Ginit, n) end - tol = eps(T)*T(10) + tol = eps(T) * T(10) max_iter = 30 * n iter = 0 ; halfN = div(n, 2) - + mem = T(1); count_static = 0 #mem and count_static allow to detect the failure of Wilkinson shifts. + start = 1 #start remembers if a zero eigenvalue appeared in the middle of ev. while n > 2 && iter < max_iter - implicitstep_vec!(ev, Qeven, Qodd, n - 1, halfN) + start = implicitstep_vec!(ev, Qeven, Qodd, n - 1, halfN, start) while n > 2 && abs(ev[n-1]) <= tol values[n] = 0 n -= 1 @@ -217,8 +223,24 @@ end eigofblock(ev[n - 1], values[n-1:n] ) n -= 2 end + if start > n-2 + start = 1 + end + + if n>1 && abs(mem - ev[n-1]) < T(0.0001) * abs(ev[n-1]) + count_static += 1 + if count_static > 4 + #Wilkinson shifts have failed, change strategy using LAPACK tridiagonal symmetric solver. + values, Q = skewtrieigen_backup!(Backup) + return Eigen(values, Q) + end + else + count_static = 0 + end + mem = (n>1 ? ev[n-1] : T(0)) iter += 1 end + if n > 0 if n == 2 eigofblock(ev[1], values[1:2]) @@ -298,6 +320,7 @@ end end @views function skewtrieigen_divided!(A::SkewHermTridiagonal{T}) where {T<:Real} + Backup = copy(A) n = size(A, 1) values = complex(zeros(T, n)) vectorsreal = similar(A, n, n) @@ -317,8 +340,10 @@ end max_iter = 30 * n iter = 0 ; halfN = div(n, 2) + mem = T(1); count_static = 0 #mem and count_static allow to detect the failure of Wilkinson shifts. + start = 1 #start remembers if a zero eigenvalue appeared in the middle of ev. while n > 2 && iter < max_iter - implicitstep_vec!(ev, Qeven, Qodd, n - 1, halfN) + start = implicitstep_vec!(ev, Qeven, Qodd, n - 1, halfN, start) while n > 2 && abs(ev[n-1]) <= tol values[n] = 0 n -= 1 @@ -327,6 +352,27 @@ end eigofblock(ev[n - 1], values[n-1:n] ) n -= 2 end + if n > 2 && abs(ev[n-1]-mem) < tol + eigofblock(ev[n - 1], values[n-1:n] ) + n -= 2 + end + if start > n-2 + start = 1 + end + if n>1 && abs(mem - ev[n-1]) < T(0.0001) * abs(ev[n-1]) + count_static += 1 + if count_static > 4 + #Wilkinson shifts have failed, change strategy using LAPACK tridiagonal symmetric solver. + Q = complex.(vectorsreal, vectorsim) + values, Q = skewtrieigen_backup!(Backup) + vectorsreal .= real.(Q) + vectorsim .= imag.(Q) + return values, vectorsreal, vectorsim + end + else + count_static = 0 + end + mem = (n>1 ? ev[n-1] : T(0)) iter += 1 end if n > 0 @@ -408,4 +454,37 @@ end error("Maximum number of iterations reached, the algorithm didn't converge") end +end + +#The Wilkinson shifts have some pathological cases. +#In these cases, the skew-symmetric eigenvalue problem is solved as detailed in +#C. Penke, A. Marek, C. Vorwerk, C. Draxl, P. Benner, High Performance Solution of Skew-symmetric Eigenvalue Problems with Applications in Solving the Bethe-Salpeter Eigenvalue Problem, Parallel Computing, Volume 96, 2020. + +@views function skewtrieigvals_backup!(S::SkewHermTridiagonal{T,V,Vim}) where {T<:Real,V<:AbstractVector{T},Vim<:Nothing} + n = size(S,1) + H = SymTridiagonal(zeros(eltype(S.ev), n), copy(S.ev)) + vals = eigvals!(H) + return vals .= .-vals +end + +@views function skewtrieigen_backup!(S::SkewHermTridiagonal{T,V,Vim}) where {T<:Real,V<:AbstractVector{T},Vim<:Nothing} + + n = size(S, 1) + H = SymTridiagonal(zeros(T, n), S.ev) + trisol = eigen!(H) + vals = complex.(0, -trisol.values) + Qdiag = complex(zeros(T,n,n)) + c = 1 + @inbounds for j=1:n + c = 1 + @simd for i=1:2:n-1 + Qdiag[i,j] = trisol.vectors[i,j] * c + Qdiag[i+1,j] = complex(0, trisol.vectors[i+1,j] * c) + c *= (-1) + end + end + if n % 2 == 1 + Qdiag[n,:] = trisol.vectors[n,:] * c + end + return vals, Qdiag end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index fb27209..af50c03 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,4 @@ -using LinearAlgebra, Random +using LinearAlgebra, Random, SparseArrays using SkewLinearAlgebra using Test @@ -480,3 +480,24 @@ end end end + +@testset "issue#118 and issue#130" begin + #issue #130 + sp = sparse([2, 7, 1, 3, 6, 8, 2, 4, 7, 9, 3, 5, 8, 10, 4, 9, 2, 7, 1, 3, 6, 8, 2, 4, 7, 9, 3, 5, 8, 10, 4, 9], [1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10], [-0.8414709848, 1.5403023059, 0.8414709848, -0.8414709848, 0.4596976941, 1.5403023059, 0.8414709848, -0.8414709848, 0.4596976941, 1.5403023059, 0.8414709848, -0.8414709848, 0.4596976941, 1.5403023059, 0.8414709848, 0.4596976941, -0.4596976941, 0.8414709848, -1.5403023059, -0.4596976941, -0.8414709848, 0.8414709848, -1.5403023059, -0.4596976941, -0.8414709848, 0.8414709848, -1.5403023059, -0.4596976941, -0.8414709848, 0.8414709848, -1.5403023059, -0.8414709848], 10, 10) + A = SkewHermitian(Matrix(sp)) + E = eigen(A) + @test E.vectors*Diagonal(E.values)*E.vectors' ≈ A + sp = sparse([26, 50, 51, 52, 27, 51, 52, 53, 28, 52, 53, 54, 29, 53, 54, 55, 30, 54, 55, 56, 31, 55, 56, 32, 33, 57, 58, 32, 33, 34, 57, 58, 59, 33, 34, 35, 58, 59, 60, 34, 35, 36, 59, 60, 61, 35, 36, 37, 60, 61, 62, 36, 37, 38, 61, 62, 63, 37, 38, 39, 62, 63, 64, 38, 39, 40, 63, 64, 65, 39, 40, 41, 64, 65, 66, 40, 41, 42, 65, 66, 41, 42, 43, 66, 42, 43, 44, 43, 44, 45, 44, 45, 46, 45, 46, 47, 46, 47, 48, 47, 48, 49, 48, 49, 50, 49, 50, 51, 1, 50, 51, 52, 2, 51, 52, 53, 3, 52, 53, 54, 4, 53, 54, 55, 5, 54, 55, 56, 6, 55, 56, 7, 8, 7, 8, 9, 58, 8, 9, 10, 59, 9, 10, 11, 60, 10, 11, 12, 61, 11, 12, 13, 62, 12, 13, 14, 63, 13, 14, 15, 64, 14, 15, 16, 65, 15, 16, 17, 66, 16, 17, 18, 17, 18, 19, 18, 19, 20, 19, 20, 21, 20, 21, 22, 21, 22, 23, 22, 23, 24, 23, 24, 25, 1, 24, 25, 26, 1, 2, 25, 26, 27, 1, 2, 3, 26, 27, 28, 2, 3, 4, 27, 28, 29, 3, 4, 5, 28, 29, 30, 4, 5, 6, 29, 30, 31, 5, 6, 30, 31, 7, 8, 7, 8, 9, 33, 8, 9, 10, 34, 9, 10, 11, 35, 10, 11, 12, 36, 11, 12, 13, 37, 12, 13, 14, 38, 13, 14, 15, 39, 14, 15, 16, 40, 15, 16, 17, 41], [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 26, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 28, 29, 29, 29, 29, 30, 30, 30, 30, 31, 31, 31, 32, 32, 33, 33, 33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 36, 36, 36, 36, 37, 37, 37, 37, 38, 38, 38, 38, 39, 39, 39, 39, 40, 40, 40, 40, 41, 41, 41, 41, 42, 42, 42, 43, 43, 43, 44, 44, 44, 45, 45, 45, 46, 46, 46, 47, 47, 47, 48, 48, 48, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 58, 58, 58, 58, 59, 59, 59, 59, 60, 60, 60, 60, 61, 61, 61, 61, 62, 62, 62, 62, 63, 63, 63, 63, 64, 64, 64, 64, 65, 65, 65, 65, 66, 66, 66, 66], [-1.176e-13, -0.0980580675690928, 0.987744983947536, -0.0980580675690911, -4.9e-14, -0.0980580675690911, 0.987744983947536, -0.0980580675690927, 1.96e-14, -0.0980580675690927, 0.987744983947536, -0.098058067569091, -1.96e-13, -0.098058067569091, 0.987744983947536, -0.0980580675690928, -1.274e-13, -0.0980580675690928, 0.987744983947536, -0.0980580675690928, -5.88e-14, -0.0980580675690928, 0.987744983947536, -10.0, -0.4902903378454601, -100.97143868952186, -0.098058067569092, 0.4902903378454601, -10.0, -0.4902903378454601, -0.098058067569092, -100.97143868952186, -0.098058067569092, 0.4902903378454601, -10.0, -0.4902903378454601, -0.098058067569092, -100.97143868952186, -0.0980580675690921, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690921, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.49029033784546, -0.0980580675690919, -100.97143868952186, -0.0980580675690923, 0.49029033784546, -10.0, -0.4902903378454601, -0.0980580675690923, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454603, 0.4902903378454603, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454602, 0.4902903378454603, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454603, 0.4902903378454603, -10.0, -0.4902903378454599, 1.176e-13, 0.4902903378454599, -10.0, -0.4902903378454602, 4.9e-14, 0.4902903378454602, -10.0, -0.4902903378454599, -1.96e-14, 0.4902903378454599, -10.0, -0.4902903378454603, 1.96e-13, 0.4902903378454603, -10.0, -0.4902903378454599, 1.274e-13, 0.4902903378454599, -10.0, -0.4902903378454599, 5.88e-14, 0.4902903378454599, -10.0, 10.0, -0.4902903378454601, 0.4902903378454601, 10.0, -0.4902903378454601, 2.4e-15, 0.4902903378454601, 10.0, -0.4902903378454601, 4.9e-15, 0.4902903378454601, 10.0, -0.4902903378454601, 7.3e-15, 0.4902903378454601, 10.0, -0.4902903378454601, 9.8e-15, 0.4902903378454601, 10.0, -0.49029033784546, 1.22e-14, 0.49029033784546, 10.0, -0.4902903378454601, 1.47e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 1.71e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 1.96e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 2.2e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 0.4902903378454601, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454603, 0.4902903378454603, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454603, 0.4902903378454602, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454603, 0.0980580675690928, 0.4902903378454603, 10.0, -0.4902903378454599, -0.987744983947536, 0.0980580675690911, 0.4902903378454599, 10.0, -0.4902903378454602, 0.0980580675690911, -0.987744983947536, 0.0980580675690927, 0.4902903378454602, 10.0, -0.4902903378454599, 0.0980580675690927, -0.987744983947536, 0.098058067569091, 0.4902903378454599, 10.0, -0.4902903378454603, 0.098058067569091, -0.987744983947536, 0.0980580675690928, 0.4902903378454603, 10.0, -0.4902903378454599, 0.0980580675690928, -0.987744983947536, 0.0980580675690928, 0.4902903378454599, 10.0, -0.4902903378454599, 0.0980580675690928, -0.987744983947536, 0.4902903378454599, 10.0, 100.97143868952186, 0.098058067569092, 0.098058067569092, 100.97143868952186, 0.098058067569092, -2.4e-15, 0.098058067569092, 100.97143868952186, 0.0980580675690921, -4.9e-15, 0.0980580675690921, 100.97143868952186, 0.0980580675690919, -7.3e-15, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -9.8e-15, 0.0980580675690919, 100.97143868952186, 0.0980580675690923, -1.22e-14, 0.0980580675690923, 100.97143868952186, 0.0980580675690919, -1.47e-14, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -1.71e-14, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -1.96e-14, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -2.2e-14], 66, 66) + A = SkewHermitian(Matrix(sp)) + E = eigen(A) + @test E.vectors*Diagonal(E.values)*E.vectors' ≈ A + #issue #118 + for v ∈ ([1.0, 0.001, 1.0, 0.0001, 1.0], [2.0, 1e-11, 2.0, 1e-11, 2.0]) + A = SkewHermTridiagonal(v) + E = eigen(A) + @test E.vectors*Diagonal(E.values)*E.vectors' ≈ A + B = SkewHermitian(Matrix(A)) + E = eigen(B) + @test E.vectors*Diagonal(E.values)*E.vectors' ≈ B + end +end