[Octave-bug-tracker] [bug #53700] eigs test failure related to ARPACK generating real NaN rather than complex NaN+1i*NaN

[Dan Sebald]

Follow-up Comment #5, bug #53700 (project octave):

OK, yeah 3.5.0+real-2 here as well (64-bit build).

Here's settings for what could be relevant:


  BLAS libraries:                -lopenblas
[snip]
  Build static libraries:               no
  Build shared libraries:               yes
  Dynamic Linking:                      yes (dlopen)
  64-bit array dims and indexing:       yes
  64-bit BLAS array dims and indexing:  no
  OpenMP SMP multithreading:            yes
  Truncate intermediate FP results:     yes


I don't have a debugger set up at the moment and it is inconvenient to jump
back and forth between systems, so let me just print out some stuff as a
start...

It appears to me that the tests jump back and forth between complex versions
of the library call and real versions of the library call.  The one that fails
is real.  I.e., in the following:


  octave_idx_type nconv;
  if (a_is_complex || b_is_complex)
    {
std::cerr << "I AM COMPLEX\n";
      ComplexMatrix eig_vec;
      ComplexColumnVector eig_val;
[snip]
      if (nargout < 2)
        retval(0) = eig_val;
      else
        retval = ovl (eig_vec, ComplexDiagMatrix (eig_val), double (info));
    }
  else if (sigmai != 0.)
    {
std::cerr << "PROMOTED REAL TO COMPLEX\n";
      // Promote real problem to a complex one.
      ComplexMatrix eig_vec;
      ComplexColumnVector eig_val;
[snip]
     if (nargout < 2)
        retval(0) = eig_val;
      else
        retval = ovl (eig_vec, ComplexDiagMatrix (eig_val), double (info));
    }
  else
    {
std::cerr << "I AM REAL\n";
      if (symmetric)
        {
std::cerr << "    SYMMETRIC\n";
          Matrix eig_vec;
          ColumnVector eig_val;
[snip]
         if (nargout < 2)
            retval(0) = eig_val;
          else
            retval = ovl (eig_vec, DiagMatrix (eig_val), double (info));
        }
      else
        {
std::cerr << "    NOT SYMMETRIC\n";
          ComplexMatrix eig_vec;
          ComplexColumnVector eig_val;
[snip]
          if (nargout < 2)
            retval(0) = eig_val;
          else
            retval = ovl (eig_vec, ComplexDiagMatrix (eig_val), double
(info));
        }


Notice in the above that the I AM REAL case is the only one which returns a
real vector, otherwise the rest return complex.  Here is the output of that
failed case:


octave:2> A = magic (100) / 10 + eye (100);
octave:3> opts.v0 = (1:100)';
octave:4> opts.maxit = 10;
octave:5> d = eigs (A, 10, "sm", opts);
ans = eigs 3
CALLED __eigs__
I AM REAL
     NOT SYMMETRIC
warning: eigs: Only 1 of the 10 requested eigenvalues converged
warning: called from
    eigs at line 268 column 18


The NOT SYMMETRIC case has 


          ComplexMatrix eig_vec;
          ComplexColumnVector eig_val;


so one would think that is the octave_value that is returned because I don't
see anywhere in which there is conditional code that would check the size of
the imaginary portion, say, and discard it if the value is negligible.

How Octave is returning a non-complex value has me wondering.  But the
interesting thing is that even though that last case (i.e., REAL and NOT
SYMMETRIC, lines 517 through 555 of __eigs__.cc) uses a ComplexColumnVector
for the return value, the routines of the ARPACK are all of the EigsRealXYZ()
variety, not complex.  (In the I AM COMPLEX and PROMOTED COMPLEX scenarios the
functions are EigsComplexXYZ().)  Could this be a compiler issue with calling
those EigsRealXYZ routines and eig_val being ComplexColumnVector?  That might
also explain why you are seeing three eigenvalues of 1 while my build only
reports one eigenvalue of 1.

    _______________________________________________________

Reply to this item at:

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

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