[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 #22, bug #53700 (project octave):

@Dmitri: Could you look at the results for the failing case for all those
versions and see if the number of converged eigenvalues is the same?  The code
is:


A = magic (100) / 10 + eye (100);
opts.v0 = (1:100)';
opts.maxit = 10;
d = eigs (A, 10, "sm", opts);


I've experimented with things like

opts.maxit = 100;
opts.tol = 1e-13;

to see at what point I begin to see all the eigenvalues converging.  I'm
surprised that I have to set the tolerance to 1e-13 to 1e-12 to get all ten
eigenvalues converging.  That's four orders of magnitude above the floating
point resolution.

Again, picking the number of eigenvalues to be large seems to make convergence
better, so it has me wondering if there is some error in the loop.  I've
experimented a bit, but can't uncover anything.  The dnaupd function:

http://www.caam.rice.edu/software/ARPACK/UG/node137.html


      F77_FUNC (dnaupd, DNAUPD)
        (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
         F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
         k, tol, presid, p, v, n, iparam,
         ipntr, workd, workl, lwork, tmp_info
         F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));


wants the programmer to supply the following when ido=1:


c          IDO =  1: compute  Y = OP * Z  and Z = B * X where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y,
c                    IPNTR(3) is the pointer into WORKD for Z.


Currently, Y = OP * Z is done using a L U decomposition approach that appears
to be a sparse routine with P[] and Q[] indexing.  The P Q pairs are pretty
much 1,1 2,2 3,3 etc. but along the way the indices get out of sync by 1, like
98 takes the place of 6.  I really don't know the sparse code well enough to
say why it is choosing the indices in that way.

I also noticed that currently there is no explicit code that sets Z = B * X
where in this case B is the identity matrix.  So, I added a little hunk of
code that will do that, and it didn't seem to affect the result any.  I
suppose the DNAUPD is setting vector Z = X by default.

Otherwise, all the code associated with EigsRealNonSymmetricMatrixShift
appears to be reasonable in terms of how the library function should be used. 
If anything, my concern would be that the L-U sparse decomposition carries out
Y = OP * Z precisely.

I think Marco's suggestion is the best at this point.  I'm sure I won't be the
only person who's system produces a single convergent eigenvalue that doesn't
come in a conjugate pair.

    _______________________________________________________

Reply to this item at:

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

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