bugGNU 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

Submitter:  Patrick Alken <psa>
Submitted:  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: 

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

fixed in commit cf8f4650b8aef8ac8172656939d25080dcbb33ef

Patrick Alken <psa>
Group administrator
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>
Group administrator
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>
Group administrator

 

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

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)

 

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)
  •  

    Follow 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-03-27 psa StatusNone Fixed
        Open/ClosedOpen Closed
    2022-03-18 mdunlap1 Attached File- Added francis.c, #52998
    2020-12-29 psa Attached File- Added eig.c, #50588

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code