Ah, there we go (ARPACK, The Next Generation)...
https://github.com/opencollab/arpackng/blob/master/SRC/znaupd.f

@Dan: that is original (old) arpack. Look for arpackng.

znaupd doc is in the header of
http://www.caam.rice.edu/software/ARPACK/UG/node138.html

I found it. which has to be set to 'LM' in the fortran code. My fault. Plain arpack "converges" to wrong eigenvalues, too. This is maybe expected, since B is not HPD. I think we should use the idea in section 3.2.2 of the documentation. I will elaborate on it. Of course, after 4.4.
@Dan: some time ago I asked for help into trying to make a review and write a sort of documentation of eigs.m>__eigs__.cc>eigsbase.cc>ARPACK, something similar to https://arxiv.org/pdf/cs/0604006. If you are interested, please privately write me.
I close the bug report.

@Dan: which documentation are you referring to? I cannot find it.

OK, the invalid number of converged Ritz makes sense.
I'm wondering about the following comment in the code:
The documentation says that when IDO=1:
That comment above makes it sound as though Z = M X, but then why would the documentation state Y = OP Z? That would be Y = inv(A) M M X. I don't follow that part of things.

@Dan: when there is an error in znaupd/zneupd (like maximum number of iterations reached) the number of converged eigenvalues is not set in ARPACK. The final printing of my fortran file (taken by an existing driver) is meant for converged situations. So, no real problem in my fortran file or ARPACK library.

@Dan: sorry for the zip file, which misses debug.h and contains subfolders. I uploaded a new one. I did not notice the problem with the number of converged ritz values: so, there is probably something in my fortran file or in ARPACK itself.
In this example sigma is zero. So inv(Asigma M) is inv(A) which is a diagonal matrix. Therefore, I divided x by the diagonal of A.
(file #44033)

@Marco: I've looked at the code just a little bit. I have no concrete answers, but just observations.
First, I agree with your thinking that the FORTRAN (nonsparse construct) and C (sparse construct) implementations should behave pretty much the same in terms of convergences. So, this is a good test.
For the FORTRAN version, I had to remove the debug header file
I then run the program and get some loop outputs for the residuals, then the final summary:
Everything looks reasonable there except the number of converged Ritz values reported is unrealistic. It should report 0. log2(538187800) is 29.004, so it's not like this could be 1 in disguise or anything.
The other thing that has me wondering is the IDO value reported keeps repeating with the pattern of 2,2,1,2,2,1,2,2,1,2,2,1...
The above does look to be about 17 repeated patterns (i.e., the number of computed OP*x). However, I'm wondering why in your code for the IDO=1 and IDO=1 cases there is the division by something resembling the eigenvalues or Ritz values (is it supposed to be (Asigma*M)^1?). Won't that throw the estimates out of sorts?
From documentation
http://www.caam.rice.edu/software/ARPACK/UG/node138.html
when IDO is 1, the routine wants
and what the code is computing is
Otherwise, all the FORTRAN vs. C indexing looks about right. The mx() routine looks to match what the eigs() example is set up as.

I run the original test (asking for 4 eigenvalues) with a plain fortran file, and arpack does not converge, which is better than silently returning the wrong eigenvalues as octave.
Arpack calls several times znaupd which returns a vector x and ask for an operation y = OP(x) (either y = B*x or y = A\x). I displayed the first element of x and y before and after each of the operations. They coincide, in plain fortran and octave, till the end of the first major iteration (about 30 calls to znaupd). Then, eigenvalues are computed. They coincide, even if they are in a different order. This is already strange, but maybe I missed a sorting operation when Arpack directly displays results. Then the first residual has to be computed. In order to do that, a B*x operation is needed. Well, at this point, x is different in plain fortran and octave. At this point, resid, which is an input/output argument for znaupd, is the same as input, but different as output. It seems that znaupd behaves differently when called by the fortran file or when called by octave. I have only one set of blas/lapack/arpack libraries on my system.
I include a zip file, with the modified eigsbase.cc which produces output, with zndrv4_53719.f (gfortran zndrv4_53719.f lblas llapack larpack), and with the output files, with an entry LOOK HERE to find the divergence point. I still have no idea.
(file #44025)

I verified that the relative errors in the original Test_GEVP are now all zero. Marking as fixed and closing report.

Marco: I imported your changes, made a couple of minor edits to the commit messages, then pushed on stable, and merged with default. All the tests for eigs pass on my system.

I decided for two patches. The first is about the segmentation fault. Moreover, it fixes some conditions in the subfunction select in eigs.m. It should certainly be applied (eigs1.diff).
The second changes the switch between using ARPACK or eig inside eigs.m. Unfortunately, it was not as easy as I thought in comment #3. In fact, we should compare p and rows (a). By default, p is 2 * k, but it can be overwritten by opts.p. I also changed a lot of BIST tests in order to really use ARPACK and not eig. Since the default number of eigenvalues is 6, and the default Krylov dimension is 2 * 6, I forced eigs to use eig when rows (a) < 13 (it was rows (a) < 9). This second patch (eigs2.diff) needs an accurate review (we should be sure that all tests pass again on different platforms) and does not solve the undetected illconditioning of the problem in this bug report.
When eig is used, the output argument flag is set to 0. What about setting it to 0 in order to know when eig was used instead of ARPACK?
(file #44003, file #44004)

For question 2:
2) the diff fixes in general the crash when B is not SPD, but maybe a bug in sparse_chol is now hidden (is it correct that fact.P() does not exist when info != 0?)
This is correct. Downstream code should not rely on any of the outputs from the factorization if the factorization has failed.
P is not, in fact a data member of the sparse_chol class, but actually a function. The code is
So P is not calculable if Lsparse hasn't been computed. There is an optional fourth argument that we could pass to the sparse_chol to force some kind of factorization. Instead of this,
call
but I doubt that is a good idea.
I will commit your patch whenever you think it is sufficiently tested.

The matrix B is not Hermitian positive definite, as required in the documentation, and this is the main problem here. Looking for small eigenvalues does not require explicit inversion of B by cholesky, and therefore this problem is not detected. Matlab succeeds when asking 5 eigenvalues, but returns an error with 4 (illconditioned problem). I believe that Matlab uses eig when asking 5 eigenvalues. This is reasonable: since the default choice is to double the number of eigenvalues for the Krylov space, when the Krylov space has the same dimension of the problem, we should use eig. This is possible by replacing
with
in eigs.m. It seems that ARPACK does not detect problems in convergence and therefore Octave silently returns wrong eigenvalues (when asking 4 eigenvalues). I will implement the example in fortran, to be 100% sure it is a notdetectedproblem in ARPACK.
For the crash, there is a problem in eigsbase:make_cholb. For the sparse versions, when the factorization fails, it seems that fact.P () is not available. The test should be as in the full case
fact.P() (which is the second output of chol) is probably set later, I don't know. I attach a diff (not yet a changeset): I had to change a lot of tests, otherwise call_eig was used. Morever, I fixed some conditions in eigs.m.
To summarize:
1) the diff fixes the problem by calling eig (which is reasonable for this example)
2) the diff fixes in general the crash when B is not SPD, but maybe a bug in sparse_chol is now hidden (is it correct that fact.P() does not exist when info != 0?)
3) the diff does not detect that the original problem is illconditioned (this is probably a defect in ARPACK, I will check).
(file #43992)

It seems a problem in _eigs_.cc. I will investigate.

Very timely bug report. We're looking at similar issues here:
https://savannah.gnu.org/bugs/index.php?53700
(but that case is one which follows the "not sparse" path).
Just noting that these three cases work:
[Evectors Evalues] = eigs(A, B, 10, 'SM');
[Evectors Evalues] = eigs(A, B, 9, 'SM');
[Evectors Evalues] = eigs(A, B, 8, 'SM');
but going to
[Evectors Evalues] = eigs(A, B, 7, 'SM');
throws the eigenvalues out of alignment. However, I think there is some requirement in the _eigs_.cc and/or eigsbase.cc code that the NEV should satisfy 2 <= NCVNEV, so perhaps internally those first three above are handled with a different algorithmdon't know.

Attached I send a very small octave script showing that eigs returns wrong results in case of a general eigenvalue problem with complex matrices.
In addition octave crashes if we search for the largest eigenvalues ('LM').
I tested it in windows as well in linux. Both with the same behaviour.
