dhseqr gives wrong result

[andreasnoack]

This program

program bug
    implicit none

    double precision :: H1(4, 4), H2(4, 4), wr1(4), wi1(4), wr2(4), wi2(4), Z, work(10), dble, epsilon
    integer*8 :: info

    H1(:,:) = 0

    H1(2, 1) = dble(z'bfdf916d32df0e1d')
    H1(3, 2) = dble(z'bf782807624514d9')
    H1(4, 3) = dble(z'bf80d94d89578784')
    H1(1, 2) = dble(z'3fdf916d32df0e1d')
    H1(2, 3) = dble(z'3f782807624514da')
    H1(3, 4) = dble(z'3f80d94d89578784')

    H2(:,:) = H1
    H2(4,4) = epsilon(H2(4,4))

    call dhseqr('E', 'N', 4_8, 1_8, 4_8, H1, 4_8, wr1, wi1, Z, 1_8, work, 10, info)
    call dhseqr('E', 'N', 4_8, 1_8, 4_8, H2, 4_8, wr2, wi2, Z, 1_8, work, 10, info)

    write(*,*) 'wr1: ', wr1
    write(*,*) 'wi1: ', wi1
    write(*,*) 'wr2: ', wr2
    write(*,*) 'wi2: ', wi2

end program bug

produces

 wr1:    0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000
 wi1:   0.73992959728054897      -0.73992959728054897        8.2263841908860064E-003  -8.2263841908860064E-003
 wr2:    0.0000000000000000        0.0000000000000000        1.1102230246251565E-016   1.1102230246251565E-016
 wi2:   0.49328639818703252      -0.49328639818703252        8.2263841908860116E-003  -8.2263841908860116E-003

on my machine with LAPACK 3.8.

I don't think there should be a conditioning issue here. I guess the convergence criterion might be confused by the zero diagonal.

[thijssteel]

i've had a crack at this bug. I'm not sure but it looks like it's caused by underflow/rounding errors in dlanv2.

The following happens during the iterations:

• If the shift is purely imaginary, the structure will be preserved
• because we compare with the diagonal elements, we would need H(3,2) to be under sfmin to get a deflation
• after 10 iterations, we use the ad hoc shift (which is not purely imaginary) and get a deflation

We end up with the following subblock for the first 2 eigenvalues:

  0.0000E+00      4.9329E-01
 -4.9329E-01     -9.8813-324

DLANV2 turns this into:

  0.0000E+00      7.3993E-01
 -7.3993E-01      0.0000E+00

The cause is a rounding error in line 218 of DLANV2, CS is sqrt(0.5), so SN should also be sqrt(0.5), but instead it becomes -1. Not even a valid rotation, let alone one that correctly normalizes the 2x2 block.

[andreasnoack]

Indeed a tiny non-zero value in the diagonal seems to completely ruin the calculation

julia> eigvals([0.0 1.0; -1.0 0.0])
2-element Array{Complex{Float64},1}:
 0.0 - 1.0im
 0.0 + 1.0im

julia> eigvals([0.0 1.0; -1.0 nextfloat(0.0)])
2-element Array{Complex{Float64},1}:
 0.0 - 0.5000000000000001im
 0.0 + 0.5000000000000001im

julia> eigvals([0.0 1.0; -1.0 1e-310])
2-element Array{Complex{Float64},1}:
 0.0 - 1.000000000000023im
 0.0 + 1.000000000000023im