# GNU Scientific Library - Bugs: bug #59781, Eigenvalue failure

You are not allowed to post comments on this tracker with your current authentication level.

# bug #59781: Eigenvalue failure

 Submitted by: Patrick Alken Submitted on: Tue 29 Dec 2020 11:38:25 PM UTC Category: Runtime error Severity: 3 - Normal Operating System: Status: Fixed Assigned to: None Open/Closed: Closed Release:
* Mandatory Fields

Comment Type & Canned Response:

## Discussion

Sun 27 Mar 2022 07:46:54 PM UTC, comment #2:

fixed in commit cf8f4650b8aef8ac8172656939d25080dcbb33ef

Patrick Alken <psa>
Thu 24 Mar 2022 11:35:32 PM UTC, comment #1:

emails from M. Dunlap mdunlap1 =at= email =dot= arizona =dot= edu

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

===================
original email
===================

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

fig1.png

fig2.png

Patrick Alken <psa>
Tue 29 Dec 2020 11:38:25 PM UTC, original submission:

The attached code demonstrates a failure of the non-symmetric eigenvalue solver. Reported by Dmitry Cheshkov (dcheshkov =at= gmail =dot= com)

Patrick Alken <psa>

## Attached Files

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Comment:

Attached Files
file #52998:  francis.c added by mdunlap1 (29KiB - text/x-csrc - The bug was in francis_search_subdiag_small_elements(). I have attached a corrected version. (With comments to show where the issue was.) This version successfully fixes this issue and passes the tests in test.c.)
file #50588:  eig.c added by psa (2KiB - text/x-csrc)

## Dependencies

Depends on the following items: None found

Items that depend on this one: None found

Carbon-Copy List
• -email is unavailable- added by mdunlap1 (Updated the item)
• -email is unavailable- added by psa (Submitted the item)
•