DGEEV function gives non-smooth results while finding eigenvalues of tri-diagonal matrix

[N0rbert]

I have old Fortran 90 program for calculating eigenvalues of tri-diagonal matrix. Technically it relates to Mathieu functions, where eigenvalues should be monothonic when plotted for different orders and values of q (see https://dlmf.nist.gov/28.2#F1 ).

I remember the program was working fine in early 2010s with reference LAPACK.
Then something was broken in BLAS or LAPACK, I can't get the really cause because of big codebase of both projects, strace, ltrace, GDB do not help me.
The last two versions where the program was working fine were Debian 6 (with lapack 3.2.1-8) and Ubuntu 12.04 LTS (with lapack 3.3.1-1). It seems to be broken first in Ubuntu 12.10 (with lapack 3.4.1-6).

I was a bit surprised that it is still broken.
The eigenvalues a_n(q) of high even order (from about >20) get unexpected spikes despite values of q.

Below is the fully reproducible way to test the program using two Docker containers with Ubuntu 12.04 LTS (to test ordinary LAPACK) and Debian testing (to test modern OpenBLAS).
The program calculates matrixes of eigenvalues for 51 orders (from 0 to 50) and for 100 values of q (from 0 to 30000) and then saves text file in GNU Octave-friendly format for post-processing. Below is the code, to be placed in file named dgeev_test.f90 in some folder:

cat <<\EOF > dgeev_test.f90
program dgeev_test
    implicit none
    integer n, m, kd, q_i, m_i, a_tab_i
    real t1, t2
    integer file_unit
    integer min_order, max_order, order_step
    doubleprecision max_q, q, a;
    doubleprecision q_array(100)
    integer orders_array(51)
    doubleprecision a_n(51,100)
    character(len=32) fmt_string_orders
    character(len=32) fmt_string_qs
    character(len=32) fmt_string_a_n

    min_order = 0
    order_step = 1
    max_order = 50

    max_q = 3e4;

    do m_i = 1,1+(max_order - min_order)/order_step    
        orders_array(m_i) = min_order + order_step*(m_i-1)
        m = orders_array(m_i)

        print *, 'm = ', m
        call cpu_time ( t1 )

        do q_i = 1,100
            q_array(q_i) = (q_i-1)*(max_q/99)
            q = q_array(q_i)

            call max_n(m, q, n)

            if (mod(m, 2) == 0) then
                kd = 1;
            else
                kd = 2;
            endif

            call fcoef(kd, m, q, n, a)
            a_n(m_i, q_i) = a
        enddo

        call cpu_time ( t2 )
        write ( *, * ) 'Elapsed CPU time = ', t2 - t1
    enddo

    file_unit=10
    open(unit=file_unit, file='fortran.txt', status='replace', action='write')

    write(file_unit, '(A)') "# Created for Octave by Fortran 90 program"
    write(file_unit, '(A)') "# name: orders_array"
    write(file_unit, '(A)') "# type: matrix"
    write(file_unit, '(A)') "# rows: 1"
    write(file_unit, '(A, I0)') "# columns: ", size(orders_array)
    write(fmt_string_orders, '(A, I0, A)') '(', size(orders_array), '(I0, 1X))'
    write(file_unit, fmt_string_orders) orders_array

    write(file_unit, '(A)') "# name: q_array"
    write(file_unit, '(A)') "# type: matrix"
    write(file_unit, '(A)') "# rows: 1"
    write(file_unit, '(A, I0)') "# columns: ", size(q_array)
    write(fmt_string_qs, '(A, I0, A)') '(', size(q_array), '(ES25.15, 1X))'
    write(file_unit, fmt_string_qs) q_array

    write(file_unit, '(A)') "# name: a_n_matrix"
    write(file_unit, '(A)') "# type: matrix"
    write(file_unit, '(A, I0)') "# rows: ", size(orders_array)
    write(file_unit, '(A, I0)') "# columns: ", size(q_array)
    write(fmt_string_a_n, '(A, I0, A)') '(', size(q_array), '(ES25.15, 1X))'
    do a_tab_i = 1, size(orders_array)
        write(file_unit, fmt_string_a_n) a_n(a_tab_i, :size(q_array))
    end do
    close(file_unit)
end program

subroutine max_n(m, q, n)
    integer n, m;
    doubleprecision q, qm;

    if (abs(q) <= 1.0d0) then
        qm=7.5+56.1*sqrt(abs(q))-134.7*abs(q)+90.7*sqrt(abs(q))*abs(q);
    else
        qm=17.0+3.1*sqrt(abs(q))-.126*abs(q)+.0037*sqrt(abs(q))*abs(q);
    endif

    n = int(qm+0.5*int(m));

    if (n > 200) then
        n = 200;
    endif

    return;
end

subroutine fcoef(kd, m, q, n, a)
    integer n, i, j, INFO, m, kd
    doubleprecision e1(n), e2(n), e3(n), s3(n,n), q, WR(n), WI(n), VL(n,n), VR(n,n), WORK(4*n), fc(n), a
    integer nn(n)
    i = 1;

    if (kd == 1) then
        e1(1:n) = q;
        e2(1:n) = 0.0;
        e3(1:n) = q;
        s3(1:n, 1:n) = 0;
        nn = (/(i,i=1,n)/);
        e2(1:n) = (2*(nn-1))**2;

        do i = 1,n
            s3(i,i) = e2(i);
            do j = 1,n
                if (j == i - 1) then
                    s3(i, j) = e3(i);
                endif
                if (j == i + 1) then
                    s3(i, j) = e1(i);
                endif
            enddo
        enddo

        s3(1,1) = epsilon(1.0);
        s3(1,2) = 2*q;

        call DGEEV( 'V', 'V', n, s3, n, WR, WI, VL, n, VR, n, WORK, 4*n, INFO )

        fc=vl(1:n,(m+2)/2) / sqrt(1+abs(vl(1,(m+2)/2))**2);
        a = wr((m+2)/2);
    elseif (kd == 2) then
        e1(1:n) = q;
        e2(1:n) = 0.0;
        e3(1:n) = q;
        s3(1:n, 1:n) = 0;
        nn = (/(i,i=1,n)/);
        e2(1:n) = (2*nn-1)**2;

        do i = 1,n
            s3(i,i) = e2(i);
            do j = 1,n
                if (j == i - 1) then
                    s3(i, j) = e3(i);
                endif
                if (j == i + 1) then
                    s3(i, j) = e1(i);
                endif
            enddo
        enddo

        s3(1,1)=1+q

        call DGEEV( 'V', 'V', n, s3, n, WR, WI, VL, n, VR, n, WORK, 4*n, INFO )

        fc= vl(1:n,(m+1)/2);
        a = wr((m+1)/2);
    endif

    if (real(fc(1)) < 0) then
        fc = -fc;
    endif
    return
end
EOF

Then we need to run Docker containers in the same folder :

docker run -it --rm -v $(pwd):/root ubuntu:precise bash -c " sed -i 's/archive/old-releases/g' /etc/apt/sources.list ; apt-get update ; apt-get dist-upgrade -y ; apt-get install -y gfortran liblapack-dev libatlas3gf-base; cd /root ; gfortran dgeev_test.f90 -llapack -Wall -Werror -o dgeev_test ; ls -l /etc/alternatives | tee ls-etc-alternatives_ubuntu-precise ; ldd ./dgeev_test | tee ldd_ubuntu-precise ; dpkg -l | tee dpkg-l_ubuntu-precise; ./dgeev_test ; mv -v fortran.txt fortran-ubuntu-precise.txt ; rm -v ./dgeev_test "

docker run -it --rm -v $(pwd):/root debian:testing bash -c " apt-get update ; apt-get dist-upgrade -y ; apt-get install -y gfortran libopenblas-dev; cd /root ; gfortran dgeev_test.f90 -lopenblas -Wall -Werror -o dgeev_test ; ls -l /etc/alternatives | tee ls-etc-alternatives_debian-testing ; ldd ./dgeev_test | tee ldd_debian-testing ; dpkg -l | tee dpkg-l_debian-testing ; ./dgeev_test ; mv -v fortran.txt fortran-debian-testing.txt ; rm -v ./dgeev_test "

And finally get the plots by the GNU Octave installed on the host system:

octave --persist --eval="load fortran-ubuntu-precise.txt; subplot(1,2,1); plot(q_array, a_n_matrix); title('Ubuntu 12.04 LTS with LAPACK'); xlabel('q'); ylabel('a_n(q)'); clear; subplot(1,2,2); load fortran-debian-testing.txt; plot(q_array, a_n_matrix); title('Debian testing with OpenBLAS'); xlabel('q'); ylabel('a_n(q)'); "

Other debug information is attached - output of ls -l /etc/alternatives, ldd dgeev_test and dpkg -l - see dgeev_test.zip.

[martin-frbg]

Okay... it would certainly help to know if it is broken with a recent version of both BLAS and LAPACK from the "netlib" Reference-LAPACK - I think the /etc/alternatives mechanism should let you test this easily in your Debian setup ?

OpenBLAS copies the LAPACK part from them and only replaces a small number of its functions, so any behavioural change (or bug) in DGEEV would be present in OpenBLAS as well.
On the other hand, a bug in one of OpenBLAS' BLAS kernels could derail a perfectly fine LAPACK, but most BLAS kernels are cpu-specific so (a) it would be useful to know your cpu type(s) and (b) the error would be expected to vary with cpu type (either by running the test on different computers, or by setting OPENBLAS_CORETYPE to force selection of a specific runtime target in the default Debian/Ubuntu build of the library)

[dasergatskov]

fcoef picks the characteristic value positionally -- a = wr((m+2)/2), which silently assumes DGEEV returns eigenvalues in ascending order. LAPACK's documentation has never promised any ordering for DGEEV. (I assume) LAPACK 3.3 the old double-shift QR just happened to deflate this tridiagonal matrix nearly in ascending order.
Attached is the mod to your program that adds subroutine, kth_smallest_idx, that index-sorts WR after each DGEEV call and selects the k-th smallest eigenvalue.

dgeev_test_fixed.f90.gz

The results (on Fedora 44 x86_64, the same also on Mac/M4 with Homebrew):

[N0rbert]

Thanks, it really works!
Just successfully tested with LAPACK, OpenBLAS and even MKL.