failing test cases

[KlausC]

I prepared some (further simple) test cases, which should all succeed - but ...
some work remaining. The failing cases are now all related to eigenvalues of 0 and zeros in the sub-diagonal.

It would be good to insert them into runtests.jl

julia> using Test

julia> @testset "eigvals" begin
       @testset "eigvals robustness $ev" for ev in ([0,0,0], [1,2,3,2,1], [0,1,0], [1,0,0], [1,0,1], [1,1,0], [0,1,1], [0,0,1], [1,1,1],[1,1,0,1,1])
           A = SkewHermTridiagonal(float.(ev))
           a = sort(eigvals(A), by = imag)
           b = sort(eigvals(im * Matrix(A)) / im, by = imag)
           @test a ≈ b
       end
       end
[...]
Test Summary:                        | Pass  Error  Total  Time
eigvals                              |    6      4     10  1.6s
  eigvals robustness [0, 0, 0]       |    1             1  0.5s
  eigvals robustness [1, 2, 3, 2, 1] |    1             1  0.0s
  eigvals robustness [0, 1, 0]       |           1      1  1.0s
  eigvals robustness [1, 0, 0]       |    1             1  0.0s
  eigvals robustness [1, 0, 1]       |    1             1  0.0s
  eigvals robustness [1, 1, 0]       |           1      1  0.0s
  eigvals robustness [0, 1, 1]       |           1      1  0.0s
  eigvals robustness [0, 0, 1]       |    1             1  0.0s
  eigvals robustness [1, 1, 1]       |    1             1  0.0s
  eigvals robustness [1, 1, 0, 1, 1] |           1      1  0.0s
ERROR: Some tests did not pass: 6 passed, 0 failed, 4 errored, 0 broken.

[KlausC]

In the meanwhile I have the impression, that the algorithm for skew symmetric real tridiagonal matrices is not reliable enough for production use. I propose to fall back to the proven algorithms for tridiagonal symmetric matrices, even if they do not exploit the inverse signed eigenvalues, which appear with a zero diagonal.
The examples that made me think so are attached below. I have the impression, that the effort to get the current algorithm robust is not worth while.

julia> using SkewLinearAlgebra
[ Info: Precompiling SkewLinearAlgebra [5c889d49-8c60-4500-9d10-5d3a22e2f4b9]

julia> A = SkewHermTridiagonal([1, 1e-4, 1, 1e-4, 1.0]);

julia> eigvals(A)
6-element Vector{ComplexF64}:
 0.0 + 1.000000005im
 0.0 - 1.000000005im
 0.0 + 1.1441228068479692im
 0.0 - 1.1441228068479692im
 0.0 + 0.43701602350395685im
 0.0 - 0.43701602350395685im

julia> Random.seed!(0)
TaskLocalRNG()

julia> A = SkewHermTridiagonal(rand(2000));

julia> eigvals(A);
ERROR: Maximum number of iterations reached, the algorithm didn't converge
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] skewtrieigvals!(A::SkewHermTridiagonal{Float64, Vector{Float64}, Nothing})
   @ SkewLinearAlgebra ~/dev/SkewLinearAlgebra/src/skeweigen.jl:129
 [3] eigvals!(A::SkewHermTridiagonal{Float64, Vector{Float64}, Nothing}, sortby::Nothing)
   @ SkewLinearAlgebra ~/dev/SkewLinearAlgebra/src/tridiag.jl:447
 [4] eigvals (repeats 2 times)
   @ ~/dev/SkewLinearAlgebra/src/tridiag.jl:484 [inlined]
 [5] top-level scope
   @ REPL[4]:1

[smataigne]

At the beginning, this package was using the technic described in https://pure.mpg.de/rest/items/item_3181584/component/file_3309065/content to solve antisymmetric eigenvalue problems. But it appeared that symmetric tridiagonal eigensolver from LAPACK itself could encounter issues with LAPACKException(22). Therefore came up the idea at the end of my internship in august to implement a specialized algorithm for skew-symmetric matrices that would be able to solve problems that LAPACK could not (not a little challenge).

For this QR algorithm, zeros are indeed difficult to handle, the main reason being that these eigenvalues do not appear by pairs during the iterations. However, for dense matrices, the performances obtained are quite nice as it can be faster than the symmetric solver as long as the non-availability of a skew-symmetric matrix-vector product for the hessenberg reduction does not become a too large disadvantage. I have only had 1 week to really work on that eigensolver, reason for which some caveats may remain for very particular cases.

I think it is worthwhile to try to make it work for these difficult cases, although it is not easy. Meanwhile, the question I am asking myself is that this eigensolver was more intended to solve dense skew-symmetric eigenvalue problems, for which it is much more robust since the caveats existing for special tridiagonal cases have a very little probability to appear. So if the efforts invested don't solve all issues with tridiagonal matrices, maybe it is worth mentioning it in the documentation and live with that caveat since tridiagonal skew-symmetric matrices don't appear by themselves that much in practice. The other possibility is to come back to the method proposed in the paper above, but caveats will remain anyway with #64.

[KlausC]

While it might be acceptable, that the algorithm diverges for some "corner cases" it is not acceptable, that it converges to false results. The second example above is such a case. It has no zeros in the subdiagonal.

[smataigne]

I could not agree more with you, reason for which I will use my next free time to find what creates these issues.

[smataigne]

#117 solves the main problem of handling zeros. I stopped trying to make zeros appear in pairs. Eigenvectors are in consequence a bit more complicated to build. For the case A = SkewHermTridiagonal([1, 1e-4, 1, 1e-4, 1.0]), an accuracy issue remains. Also I stopped shifting if the last two elements are equal.

A pathological problem of the double shift QR algorithm remains: if several zeros are consecutive in the middle of the matrix, the bulge disappears and the second part of the matrix is not updated. This should be tackled by performing classical iterations at some occasions.

[smataigne]

I close this issue since the remaining problem is not the same anymore.