From 8cfb93bf0203c5d3d66048447ffac81fe4beb5de Mon Sep 17 00:00:00 2001 From: smataigne Date: Thu, 21 Dec 2023 14:15:15 -0800 Subject: [PATCH 1/7] solve issues #118 and #130 --- .vscode/settings.json | 3 + src/skeweigen.jl | 150 ++++++++++++++++++++++++++++++++---------- test/runtests.jl | 21 +++++- 3 files changed, 137 insertions(+), 37 deletions(-) create mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..1d3a6e6 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,3 @@ +{ + "julia.environmentPath": "c:\\Users\\simon\\Documents\\Mes documents\\Doctorat\\SkewLinearAlgebra.jl" +} \ No newline at end of file diff --git a/src/skeweigen.jl b/src/skeweigen.jl index 1a897ee..9470c25 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 abs(mem - ev[n-1]) < T(1e-4) * 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 = ev[n-1] 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,10 @@ end ζ = ev[i+2] ev[i+2] *= c bulge = s * ζ + if abs(bulge) < tol && abs(ev[i]) < tol + start = i + 1 + return start + end end Q = (isodd(i) ? Qodd : Qeven) k = div(i+1, 2) @@ -176,19 +189,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 +208,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,6 +224,20 @@ end eigofblock(ev[n - 1], values[n-1:n] ) n -= 2 end + if start > n-2 + start = 1 + end + if abs(mem - ev[n-1]) < T(1e-4) * 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 = ev[n-1] iter += 1 end if n > 0 @@ -298,6 +319,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 +339,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 +351,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 abs(mem - ev[n-1]) < T(1e-4) * 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 = ev[n-1] iter += 1 end if n > 0 @@ -408,4 +453,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..710a969 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,22 @@ 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 From abb71859a9d2d446fb855c4117a09442f7997870 Mon Sep 17 00:00:00 2001 From: smataigne Date: Thu, 21 Dec 2023 14:40:46 -0800 Subject: [PATCH 2/7] solve issues #118 and #130 --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 710a969..b74b319 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,4 @@ +import Pkg; Pkg.add("SparseArrays") using LinearAlgebra, Random, SparseArrays using SkewLinearAlgebra using Test From 05964c8d60bc29813660ea6f7ce33a00bd9f1a9e Mon Sep 17 00:00:00 2001 From: smataigne Date: Thu, 21 Dec 2023 14:56:50 -0800 Subject: [PATCH 3/7] solve issues #118 and #130 --- Project.toml | 3 ++- test/runtests.jl | 1 - 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 2dba0f6..269108c 100644 --- a/Project.toml +++ b/Project.toml @@ -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/test/runtests.jl b/test/runtests.jl index b74b319..710a969 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,3 @@ -import Pkg; Pkg.add("SparseArrays") using LinearAlgebra, Random, SparseArrays using SkewLinearAlgebra using Test From 7b915409b4cfea3354bbe3fee18396ca840cd6ab Mon Sep 17 00:00:00 2001 From: smataigne Date: Thu, 21 Dec 2023 16:37:51 -0800 Subject: [PATCH 4/7] solve typo --- src/skeweigen.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/skeweigen.jl b/src/skeweigen.jl index 9470c25..9bfa804 100644 --- a/src/skeweigen.jl +++ b/src/skeweigen.jl @@ -129,7 +129,7 @@ end if start > n-2 start = 1 end - if abs(mem - ev[n-1]) < T(1e-4) * abs(ev[n-1]) + 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. @@ -139,7 +139,7 @@ end else count_static = 0 end - mem = ev[n-1] + mem = (n>1 ? ev[n-1] : T(0)) iter += 1 end if n == 2 @@ -178,7 +178,6 @@ end bulge = s * ζ if abs(bulge) < tol && abs(ev[i]) < tol start = i + 1 - return start end end Q = (isodd(i) ? Qodd : Qeven) @@ -227,7 +226,8 @@ end if start > n-2 start = 1 end - if abs(mem - ev[n-1]) < T(1e-4) * abs(ev[n-1]) + + 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. @@ -237,9 +237,10 @@ end else count_static = 0 end - mem = ev[n-1] + mem = (n>1 ? ev[n-1] : T(0)) iter += 1 end + if n > 0 if n == 2 eigofblock(ev[1], values[1:2]) @@ -358,7 +359,7 @@ end if start > n-2 start = 1 end - if abs(mem - ev[n-1]) < T(1e-4) * abs(ev[n-1]) + 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. @@ -371,7 +372,7 @@ end else count_static = 0 end - mem = ev[n-1] + mem = (n>1 ? ev[n-1] : T(0)) iter += 1 end if n > 0 From f091d8b5cf564906a1e6566cb2cf510cb72c42c3 Mon Sep 17 00:00:00 2001 From: smataigne Date: Thu, 21 Dec 2023 16:43:54 -0800 Subject: [PATCH 5/7] solve typo in new tests --- test/runtests.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 710a969..af50c03 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -480,6 +480,7 @@ 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) @@ -492,10 +493,11 @@ end @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 + 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 From 15fb40f03b1854be139b69bfeb32ea5bb1c7f574 Mon Sep 17 00:00:00 2001 From: smataigne Date: Thu, 21 Dec 2023 17:41:56 -0800 Subject: [PATCH 6/7] New version 0.1.1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 269108c..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" From 6287b37f2e8e5880ef7f59cd52101a80748c6c15 Mon Sep 17 00:00:00 2001 From: smataigne Date: Thu, 21 Dec 2023 17:49:13 -0800 Subject: [PATCH 7/7] New version 0.1.1 --- .vscode/settings.json | 3 --- 1 file changed, 3 deletions(-) delete mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json deleted file mode 100644 index 1d3a6e6..0000000 --- a/.vscode/settings.json +++ /dev/null @@ -1,3 +0,0 @@ -{ - "julia.environmentPath": "c:\\Users\\simon\\Documents\\Mes documents\\Doctorat\\SkewLinearAlgebra.jl" -} \ No newline at end of file