JuliaLinearAlgebra · smataigne · Dec 22, 2023 · Dec 21, 2023 · Dec 21, 2023 · Dec 21, 2023
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "SkewLinearAlgebra"
 uuid = "5c889d49-8c60-4500-9d10-5d3a22e2f4b9"
 authors = ["smataigne <[email protected]> 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
@@ -67,23 +67,21 @@
 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 @@
             ζ = 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 @@
     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 @@
             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
 
-@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 @@
             ζ = 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 @@
             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 @@
         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 @@
             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
 
 @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 @@
     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 @@
             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 @@
         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
diff --git 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