Force the residual vector to be in the range of OP. (#79)

[caliarim]

Pull request for issue (#79) (force the initial residual to be in the range of the operator).

[sylvestre]

Many thanks for your PR!

Seems that this is failing:
https://travis-ci.org/opencollab/arpack-ng/jobs/333185780
The following tests FAILED:
77 - bug_79_double_complex (Not Run)

[caliarim]

I forgot to add the test to add_dependencies in CMakeLists.txt?

[sylvestre]

@turboencabulator @thrasibule
any opinion on this?

[turboencabulator]

Build script changes look good to me.

[thrasibule]

Why is this needed? We only end up here if ido != 0.

[caliarim]

I'm away from my laptop, but aftet the very first OP application ido is -1, isn't it? So we need to copy the result of OP application back to resid.

[thrasibule]

This seems pretty low risk, and it fixes a test case, so looks good to me.

[caliarim]

Shall we change the documentation somewhere about this change? For instance, in X{s,n}aupd it is written

c RESID Complex*16 array of length N. (INPUT/OUTPUT)
c On INPUT:
c If INFO .EQ. 0, a random initial residual vector is used.
c If INFO .NE. 0, RESID contains the initial residual vector,
c possibly from a previous run.

No mention to the range of the operator. Basically the same in X{s,n}aup2. In Xgetv0 it is written

c Force the residual vector to be in the range of the operator OP.

In the book, page 72, it is written

In either case, it the (shift-invert) computational mode calls for it, the vector is also forced to be in the range of OP.

But with this change the vector is forced to be in the range for the regular mode, too. And it has always been forced to be in the range for the generalized problem.
And what happens if somebody really would like to test her/his own inital vector?

[sylvestre]

@caliarim is that ready to be merged? Thanks

[caliarim]

I think so. With respect to my previous comment, in Xgetv0 it is written that the initial vector is in the range, and now this is (always) true. And of course we cannot change the book...

[sylvestre]

thanks!

[ViralBShah]

Apologies for chiming in on a closed PR. For Julia's Arpack.jl, we believe that this has introduced a bug (it is our best guess):

The test case in JuliaLinearAlgebra/Arpack.jl#118 demonstrates the following behaviour:

1. Fails about 2% of the time with arpack-ng 3.5. If initialized with a random v0, it passes reliably.
2. Fails always on later versions of arpack-ng, including 3.8, even when initialized with a random v0 (demonstrated in Test ARPACK 3.8 JuliaLinearAlgebra/Arpack.jl#145)

[sylvestre]

@ViralBShah nice to see you there and no worries for this.

@caliarim do you want to propose a PR ? :)
thanks

[fghoussen]

@ViralBShah: could revert both d04bdf1 and 31854ca related commits and check things are getting better at your side? Making sure the reg you see comes from that

AFAIU, this would imply reopening #79 and #142!...