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