Hello, please forgive the briefness of my previous email, I was rather
excited.
Allow me to introduce myself.
My name is Michael Dunlap. I graduated from the University of Arizona in
2018 with a
bachelor's degree. I majored in math and minored in computer science.
Having made use of
the GNU Science Library, I thought it would be worth checking to see if
there was a way I
could make a contribution to the project. The website says that anyone
interested in
contributing should start with working on the issues mentioned in the bug
tracker. Being
acquainted with the QR algorithm I decided to work on bug #59781.
Thank you very much to Patrick Alken for including well written commentary
in the code.
Without the guidance to Matrix Computations and DLAHQR, working with the
code would have
been substantially more difficult.
Part of algorithm 7.5.2 from Matrix Computations involves looking for
subdiagonal entries
that are sufficiently close to zero. If any are found we can break the
problem up into two
smaller subproblems ("deflation").
In the book Matrix Computations the criterion by which a subdiagonal entry
is deemed small
enough is:
|h_i,i-1| <= tol * ( |h_i,i| + |h_i-1,i-1| )
where |.| refers to absolute value, and tol is a suitably small number. In
the code the
value GSL_DBL_EPSILON is used for tol.
Here is a snippet of output generated by adding print statements to the
code that shows the
problem quite clearly:
[image: fig1.png]
(fig 1)
In the above we are checking the subdiagonal for sufficiently small
entries. The program
fails to find any and proceeds to perform a francis qrstep using Wilkinson
shift.
However it is quite clear that the subdiagonal entry -2.47033e-323 is small
enough
according to the aforementioned criterion. The program should have set it
to zero.
The precise nature of the bug can be stated succinctly as:
In the function francis_search_subdiag_small_elements the variables del and
dpel are
backwards.
Here is the original code:
static inline size_t
francis_search_subdiag_small_elements(gsl_matrix * A)
{
const size_t N = A->size1;
size_t i;
double dpel = gsl_matrix_get(A, N - 2, N - 2);
for (i = N - 1; i > 0; --i)
{
double sel = gsl_matrix_get(A, i, i - 1);
double del = gsl_matrix_get(A, i, i);
if ((sel == 0.0) ||
(fabs(sel) < GSL_DBL_EPSILON * (fabs(del) + fabs(dpel))))
{
gsl_matrix_set(A, i, i - 1, 0.0);
return (i);
}
dpel = del;
}
return (0);
} /* francis_search_subdiag_small_elements() */
At first del is set to h_N-1,N-1 and dpel is set to h_N-2,N-2. In the
second pass of the
loop when the comparison is made del = h_N-2,N-2 and dpel is set to
h_N-1,N-1.
This does not match algorithm 7.5.2 from Matrix Computations.
Going back to the snippet we see that the program is comparing
-2.47033e-323 to the two
zero entries on the diagonal. Hence "deflation" of the problem fails to
occur. The
algorithm should be comparing -2.47033e-323 to 1.21667e-16 and 0. Instead
it is comparing
-2.47033e-323 to the two zero valued diagonal entries. "Circling" the entry
for sel with
<,> and the entries for del and dpel with [,] we have the following:
[image: fig2.png]
(fig 2)
Now, the determination of "close enough to zero" is not necessarily unique.
So perhaps it could be argued that this
"bug" is not really a bug but an alternative way of making such
determination.
As has already been mentioned, gsl_eigen_nonsymm works on the transpose of
the matrix in
the bug submission. Additionally, balancing the matrix prior to handing it
off to gsl_eigen_nonsymm also results in
getting the correct eigenvalues.
Sincerely,
Michael Dunlap
---------- Forwarded message ---------
From: Michael Dunlap <mdunlap1@email.arizona.edu>
Date: Fri, Mar 18, 2022 at 12:23 AM
Subject: bug #59781
To: <bug-gsl@gnu.org>
Hello!
I believe that I have found the bug and fixed it.
The bug was pertaining to the variable dpel in the function:
francis_search_subdiag_small_elements.
I have attached a corrected version with comments to show where the issue
was.
I have tested this version of francis.c against the matrix submitted by
Dmitry Cheshkov and also checked it against test.c, both were successful.
Thanks,
Michael Dunlap