[Octave-bug-tracker] [bug #50846] Unexpected results with ordered qz decomposition

[Rik]

Update of bug #50846 (project octave):

        Operating System:       Microsoft Windows => Any                    

    _______________________________________________________

Follow-up Comment #2:

This looks like it is probably roundoff error, rather than a bug.  The sixth
eigenvalue is only missing being equal to 1 by eps/2.  Once you are below
machine precision, odd things can happen.

Try this


load octave_qz.mat
format bank
[ss,tt,w,eigvalS] = qz(e,d,'S');
[ss,tt,w,eigvalB] = qz(e,d,'B');
abs (eigvalS)
ans =

                   0.00
                   0.43
                   0.80
                   0.90
                   1.26
                   1.00
                   1.12
   11088015243066408.00
                   2.36
                   1.01

size (eigvalS)
ans =

   10.00    1.00

abs (eigvalB)
ans =

   1.26
   1.00
   1.12
   2.36
   1.01
   0.00
   0.43
   0.80
   0.90

size (eigvalB)
ans =

   9.00   1.00



The bank format is just a quick hack to force Octave not to use scientific
notation.  The sixth eigenvalue really looks like it is >=1 and is being
included in the larger block.  Also note that if you reverse the sorting by
using 'B' then there are only 9 eigenvalues returned.  The one that is
essentially infinite was discarded.

The code for this is in libinterp/corefcn/qz.cc if you want to take a look. 
Around line 698 it eventually calls a Fortran function


      F77_XFCN (dsubsp, DSUBSP,
                (nn, nn, aa.fortran_vec (), bb.fortran_vec (),
                 ZZ.fortran_vec (), sort_test, eps, ndim, fail,
                 ind.fortran_vec ()));


Importantly, the precision to solve for (the variable eps) is not necessarily
the same eps as the Octave interpreter uses.  It was calculated earlier as


      double inf_norm;

      F77_XFCN (xdlange, XDLANGE,
                (F77_CONST_CHAR_ARG2 ("I", 1),
                 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm
                 F77_CHAR_ARG_LEN (1)));

      double eps = std::numeric_limits<double>::epsilon () * inf_norm * nn;





    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?50846>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/