Marco: I pushed the changeset from comment #32 to stable.
Can this report be closed now, or is more work needed for the 4.4 release?

I would stay with if (isreal (d))... since we already did it. I attach a changeset.
(file #44034)

Is more work on this problem needed for 4.4?

The warning can be turned off with:
So the suggestion of
won't work. How about floor(length(A)/4) rather than floor((length(A)1)/2)? That would put the NEV somewhere in the middle of the range where the Krylov subspace isn't adequate.
It's up to Marco to make the choice. Check for real/imaginary or pick an NEV somewhere that generates lots of 1.0+iZi and NaN+iNaN. Probably doesn't matter too much. We're seeing so much variation in the convergence behavior for all the different libraries and compilers out there.

After the patch on the stable branch, I see the following behavior"
I begin to see nonNaN results at i=49:

Do the two patches pushed today for bug #53719 also resolve this report?

OK, got it.
Note that "lm" ends up being a different case; it doesn't necessarily run the same ARPACK function or the same implementation of that function. So let's not pursue that idea.
Also, an observation here is that on Mint 18.2 I find that ARPACK issues error 14 when selecting NEV to be 8 or less. On Ubuntu 18.04 (the system that produces one convergent real 1.0 eigenvalue) it fails for NEV equal to 5 or less:
So, lots of library/compiler variation going on. (I don't even see a 14 error listed in the documentation I cited previously.)
Nonetheless, I suspect that the real/complex issue is more likely to occur when NEV is just above that transition where it doesn't produce an error. A more assured way of always getting a complex result would be to push the NEV higher:
That is, the NEV=49 case has so many converged eigenvalues of 1.0 that there is bound to be complex conjugate pairs, and also plenty of NaN + NaNi as well. So the following would also be an easy way of fixing this test:

@Dan: for the "lm" version, I have real NaNs both with 9 and 10.
NEV to 50 is p to 100.

Oh, so the reason for the less robust system for NEV=10 is that the Krylov subspace dimensionality is not even near the cardinality of the largest eigenvalue group, which is 97 in this case. I.e., 97 eigenvalues of 1.0. That's what the algorithm is struggling with. On the other hand, had we chosen a system for which there were more variation in the eigenvalues, e.g., 1.0 1.0 3.0 3.0 3.5 22j 2+2j, etc., and used a Krylov subspace of 10 then we'd get good convergence.
Going the other direction, increasing NEV to 50 (or is it p=50?) makes the Krylov subspace nearer to the size of the largest eigenvalue group in this case.
Thanks for the explanation.
Would there be any other type of more assured test? Say, if instead of looking for "sm" eigenvalues, we looked for "lm". The largest eigenvalues have a group cardinality of 1 or 2 and therefore those few really large eigenvalues should converge nicely on all systems, compilers and linear algebra libraries with a choice of NEV=10. I'm seeing
But that too ends up being not complex, so I don't know if that gets us anywhere different.
How about choosing an odd number for NEV? Will that force there to be a conjugate pair somehow?
What do others see for such tests? Is it just as arbitrary?

@Dan: I try to make a cumulative answer to your deep analysis. The main "problem" in ARPACK is a proper choice of nev (number of desired eigenvalues, k in eigs) and ncv (dimension of the Krylov subspace of approximate eigenpairs (p in eigs). As written at the end of section 2.3.3 of Arpack User's Guide, "it is best to avoid setting nev in such a way that will split clusters of eigenvalues". The choice 10 is then to avoid (but it was made in order to produce a failure!). Later "as a rule of thumb, ncv >= 2*nev is reasonable". ncv=2*nev is the choice made in eigs. Therefore, requiring 10 eigenvalues correspond to using 20 as Krylov subspace. Requiring 50 eugenvalues correspond to using 100 as Krylov subspace, which is the dimension of the original matrix: therefore, you are pretty sure that you find all the eigenvalues. You can force a new ncv by setting opts.p="something larger than 20" and keeping k=10. Solution will improve.
Z = B * X (that is Z = X, in our case) is done in dnaup2, which is called be dnaupd.
Matlab can find only few eigenvalues, too.
Therefore, although it is not possible to exclude a hidden bug in the chain eigs.m, _eigs_.cc, eigsbase.cc, arpack, "Why when we request less eigenvalues do we have a less robust situation?" can be explained as above and it is not a clear symptom of a bug.

I got 4 with atlas and 2 with ref blas and openblas.
Dmitri.


@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:
I've experimented with things like
opts.maxit = 100;
opts.tol = 1e13;
to see at what point I begin to see all the eigenvalues converging. I'm surprised that I have to set the tolerance to 1e13 to 1e12 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
wants the programmer to supply the following when ido=1:
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 LU 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.

it works for me wit hany blas:
Notice the times are different wrt blas version  atlas seems faster than ref/openblas. So it is possible that actual solutions are different. (I think the size of the problem is too small to judge the efficiency of particular blas implementation.)
Dmitri.


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, nonsparse 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:
I suppose that is a crude guess as to whether something might be sparse. So, let's try
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.
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:
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.

At first I thought there was something wrong about the system setup, i.e., this large matrix with numbers much larger than eye() caused an instability somewhere. However, experimenting a bit, I don't quite understand what is happening with the number of eigenvalues. Going back to the original test:
However, change the requested number of eigenvalues from 10 to 100, i.e., all the eigenvalues and surprisingly there is no convergence problem:
Doing sum(A), every value is 50006, and it probably isn't a coincidence that the last eigenvalue in the list is 50006. Seems reasonable.
So I'm confused. Why when we request less eigenvalues do we have a less robust situation? Typically when one puts a restriction on an algorithm it leads to improved robustness...say if we know some eigenvalues are problematic we might choose sigma and the number of eigenvalues in a way that squeezes those problematic eigenvalues out of the picture.

I think Marco's fix is the right one. The purpose of the test is to verify that the function was asked for 10 eigenvalues, a vector of length 10 is returned, and the last one is either NaN or (NaN,NaN). It doesn't matter whether 1 or 5 valid eigenvalues were returned, as long as it's less than 10, and it doesn't matter what their values are.

@Marco: I don't see any problems with checking that the eigenvalues are within a certain tolerance. When it comes to applications, users often have to do such things in order to interpret results.

@Mike: I didn't see your post about Ubuntu 18.04 until now. Your result is different than what you reported from a different system, so technically it isn't reproduced. On the older system you have
and on Ubuntu 18.04 you have
So, I suspect that the algorithm is finding eigenvalues by iteratively solving some equation that could result in a negative residual somewhere that results in some imaginary pair, and if said residual is positive the result is a single real value. Just guessing, but these numerical algorithms are prone to this sort of thing with a large space and lots of iterations.
In other words, my setup could be slightly different in terms of compilers and libraries I installed that gives me a real 1.0 eigenvalue in the first location rather than a complex 1.0 + 0.0 pair.
When I think of this test, I don't know what it is supposed to reveal. It's sort of testing that realm of the algorithm where there is numerical breakdown, i.e., it's analyzing the limitation of the algorithm and how it falls apart for a larger, overparameterized system. Just exactly how the algorithm starts to fail is highly specific to the numerical methodsvery chaotic in the sense that one tiny perturbation will produce some completely different way it fails. And to illustrate, let's change the allowed iterations on my system
The number of iterations is almost like a random stopping criteria.
Wouldn't it be more useful to have a test that is acting in a more robust region of the algorithm?

@Dan: on my system, d(1) is
so, it is not possible to assert (d(1), 1). If we introduce a tolerance, then other bug reports will pop up sayng that the tolerance is not met. I would stay with if (isreal (d)) ...

@Marco: Good suggestion. I reached a similar conclusion because I think there isn't anything wrong with setting up the routines. It is simply that the library function is behaving slightly different in terms of convergence:
ip(4) is 1, not 3 in my case. It could be that convergence of the first one is solid, but the other two are tenuous and dependent on the compiler. Looking at what Mike listed:
if only that first eigenvalue is examined, then the followup code that checks for the imaginary component is not going to set "have_cplx_eig = true;". It's only the second and third cases where have_cplx_eig is set true in Mike's case.
In addition to the isreal() check, how about testing that d(1) == 1.0 or (1.0,0.0) in addition to the check on d(10)? It seems that's the reliable eigenvalue.

There is already a test in eigs of type
Look for bug #49523, starting from comment #41. It was a blasdependent result. So, the easy fix is to introduce ifelseendif in this test, too.
@Dmitri: can you check with different blas versions (ref, openblas, atlas)?

I can't reproduce on Ubuntu 18.04. Building with OpenBLAS, all tests pass, and I get

Looking at the code for ARPACK, I think the 'rvec' parameters are inputonly, so the changes you have made should be correct. But if they don't make any difference in your case, then needs more investigation.
I don't have a Ubuntu 18.04 to test on at the moment, but I'm sure I will at some point, so I can at least try to reproduce this.

Sorry, I didn't see the attachment until now. Yes, you have modified all of the functions as I was suggesting. I'm not sure if rvec is an in/out parameter or if it's inputonly. Either way, this was only a guess, since I have seen other problems with arrays of C++ bools being passed in to Fortran routines.

Dan  no, I meant it is following the right path. It's a very large ifelse so it was hard to associate your debug prints with line numbers. But it's supposed to be the condition where A is not sparse, not complex, A is not a function handle, sigma is zero, have_sigma is true, symmetric is false. I think it's the case you labeled " NOT SYMMETRIC".
I think you changed the wrong routine in eigsbase.cc, it should be 'EigsRealNonSymmetricMatrixShift'.
Since it works fine in Debian and in Fedora with the same version of ARPACK, this could be due to some difference like default compiler options used in Ubuntu (which have definitely made a difference in other situations).

FWIW  on Fedora 28 arpck 3.5.0  no problem
Dmitri.


You meant "not be doing the right thing"? I don't see the
std::cerr << "PROMOTED REAL TO COMPLEX\n";
printed out.
Anyway, I see what you mentioned, in the eigsbase.cc routine. Attached is the changes I've currently made. Are those F77_LOGICALs what you had in mind? I.e., basically do the following in all locations:
Those mods don't change the outcome.
(file #43976)

This case should have fallen into the "promoted real to complex" case, it seems to be doing the right thing for you.
I just noticed that the 'rvec' parameter is passed in as a C++ 'bool' type. But it is not properly converted to a 'F77_LOGICAL' type, it is just thrown on the call stack for each of the ARPACK routines DSEUPD, DNEUPD, ZNEUPD. Fixing that may resolve this for you.

OK, yeah 3.5.0+real2 here as well (64bit build).
Here's settings for what could be relevant:
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:
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:
The NOT SYMMETRIC case has
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 noncomplex 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.

I don't think you need to post on a bunch of other bugs, this is a relatively minor and wellcontained issue.
What you can do is post whether your system is 32bit or 64bit, whether you are using BLAS or ATLAS or OpenBLAS with Octave, and whether you have configured Octave for 64bit indexing or not.
Then if you want to debug this, you can step into the '__eigs__' function with gdb and see where the results you are getting diverge from another system where this test does not fail.

The library you installed is (I think) just a wrapper around the real ARPACK library, not a reimplementation, so you should have both installed. The package you really need is libarpack2dev, the '++' is not needed and not used by Octave.

Hmmmm, :
What library are we looking at? When I installed all the packages needed to make Octave build, I did a search for ARPACK and picked whatever appeared in the Muon (is it?) package manager. I see nothing there about 3.x.x series. Instead, what I've picked is
Have I picked the wrong packages?

With ARPACK 3.5.0 on Debian (Debian package version 3.5.0+real2, same as Ubuntu dev), I get no failures:

I have one test failing on a very uptodate Ubuntu 18.04 with KDE (not that the latter matters, I suppose), and it seems along the lines of this bug report. If it isn't the same, I'll create a separate bug report.
I will isolate just the failing test:
and a more detailed breakdown
So, it appears that there is nonconvergence but this is probably expected behavior because the test is disabling the printout of warning messages about nonconvergence. It may be the only issue is that ARPACK is generating "real" (float) values where another library is generating "complex" values. Perhaps there is a change in recent versions of ARPACK, don't know. If you want me to investigate further, I can.
