octave-bug-tracker
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

## [Octave-bug-tracker] [bug #53700] eigs test failure related to ARPACK ge

 From: Dan Sebald Subject: [Octave-bug-tracker] [bug #53700] eigs test failure related to ARPACK generating real NaN rather than complex NaN+1i*NaN Date: Fri, 20 Apr 2018 14:40:37 -0400 (EDT) User-agent: Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:59.0) Gecko/20100101 Firefox/59.0

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

I see there is a difference in algorithms when running

d = eigs (A, 10, "sm", opts);

versus

d = eigs (A, 100, "sm", opts);

The former doesn't use ___eigs___.  Rather it uses the, presumably, non-sparse
eig() routine.  So that is one difference, but more investigation shows the
situation hasn't changed.

d = eigs (A, 100, "sm", opts);
d = eigs (A, 99, "sm", opts);
d = eigs (A, 98, "sm", opts);

use the eig() routine.  For any k <= 97, the routine that is used is
___eigs___().  In eigs.m, here is the hunk of code that makes that decision:

if (rows (a) - k < 3)
call_eig = true;
endif

I suppose that is a crude guess as to whether something might be sparse.  So,
let's try

octave:42> d = eigs (A, 97, "sm", opts);
octave:44> size(d(d==d))
ans =

97    1

octave:45> d(97)
ans =  1.0000e+00 + 2.1304e-13i

The algorithm behaves well when k is large (those last three eigenvalues are
the non 1.0 values).

Where does the algorithm seem to break?  At 50, which happens to be half the
number of columns or eigenvalues, and that probably isn't a coincidence
either.

octave:63> d = eigs (A, 50, "sm", opts);
octave:64> d = eigs (A, 49, "sm", opts);
warning: eigs: Only 43 of the 49 requested eigenvalues converged
warning: called from
eigs at line 268 column 18

and it so happens that the number of iterations through that do{} loop jumps
from 100 to 243 before the ido=99 comes back when transitioning from k=50 to
k=49.

I've got a feeling there is some type of bug in Octave's do{} loop.  In the
following do loop there are several comments about exactly where some
important information lies in Eig ARPACK result:

do
{
std::cerr << "counterval " << counterval++ << "\n";
F77_INT tmp_info = octave::to_f77_int (info);

// On exit, ip(4) <= k + 1 is the number of converged eigenvalues.
// See dnaupd2.
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));
// k is not changed

std::cerr << "ido = " << ido << "\n";
info = tmp_info;

if (disp > 0 && ! octave::math::isnan (workl[iptr (5)-1]))
{
if (iter++)
{
os << "Iteration " << iter - 1 <<
": a few Ritz values of the " << p << "-by-" <<
p << " matrix\n";
if (ido == 99) // convergence
{
os << "    " << workl[iptr(5)+k] << "\n";
for (F77_INT i = 0; i < k; i++)
os << "    " << workl[iptr(5)+i-1] << "\n";
}
else
{
// the wanted Ritz estimates are at the end
for (F77_INT i = p - k - 1; i < p; i++)
os << "    " << workl[iptr(5)+i-1] << "\n";
}
}

// This is a kludge, as ARPACK doesn't give its
// iteration pointer.  But as workl[iptr(5)-1] is
// an output value updated at each iteration, setting
// a value in this array to NaN and testing for it
// is a way of obtaining the iteration counter.
if (ido != 99)
workl[iptr(5)-1] = octave::numeric_limits<double>::NaN ();
}

Could there be something strange going on with the location of some of these
results when k in the ARPACK is below N_cols/2 as opposed to greater than
N_cols/2?  And that information is being used to incorrectly "prime" the next
pass of the loop?

I don't have time to look into this, but it sure seems to me that the Eigs
ARPACK routines shouldn't be failing for this example, i.e., there should be
all convergent eigenvalues and no NaNs.

_______________________________________________________

Reply to this item at:

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

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

```

reply via email to

 [Prev in Thread] Current Thread [Next in Thread]