Reported by Edouard Canot, Oct 24, 2013
At least for one SPD matrix (A_7x7_SPD.dat, attached), the dsaupd/dseupd routines called with the 'BE' option (Both-End eigenvalues) fails. Actually, it returns the 'LM' (Largest Magnitude eigenvalues). After some investigation, it appears that the problem is in the dsaup2 routine where, at the end of the Lanczos iteration, the swap of the two parts (both end) for the requested ritz values is done. In this swap, a new value of 'np' must be used, based on 'nev0' and not on 'nev' (because the rule kplusp=nev+np must always hold, if 'nev0' is used, 'np' must be changed accordingly). I can't believe that such a bug remains in Arpack for a so long time! The 'dsaup2.f.patch' is a proposed patch. It has been created by the command: $ diff -u dsaup2.f dsaup2_new.f > dsaup2.f.patch The bug can also be reproduced directly with Scilab-5.4.1, using the same matrix and the commands stored in the file 'bug_dsaup2.sci'. Curiously enough, Octave-3.6.4 doesn't have this problem and returns, as Matlab-8.2, the correct results. ------ Steps to reproduce the problem: 1. Read the attached A_7x7_SPD.dat matrix in your program or in Scilab-5.4.1; 2. Ask for some eigenvalues (2, 3 or 4) with the 'BE' option Expected result: 0.25222 0.25222 1.37395 12.75510 Actual result: 0.41465 1.37395 1.37395 12.75510
Comment 1 by Edouard Canot, Oct 24, 2013
If the bug is confirmed, it affects also the ssaup2 routine (single precision), because the code is the same.
Comment 2 by Edouard Canot, Oct 24, 2013
As a side-effect, the modification of np (fourth line below) eliminates also a run-time error, which may occur when (old value of) np is null. Let's consider the lines 549-557 of the patched dsaup2.f source file, reproduced hereafter: nevd2 = nev0 / 2 nevm2 = nev0 - nevd2 if ( nev .gt. 1 ) then np = kplusp - nev0 call dswap ( min(nevd2,np), ritz(nevm2+1), 1, & ritz( max(kplusp-nevd2+1,kplusp-np+1) ), 1) call dswap ( min(nevd2,np), bounds(nevm2+1), 1, & bounds( max(kplusp-nevd2+1,kplusp-np+1)), 1) end if When np=0, kplusp-np+1 is strictly greater then kplusp, which is the size of the ritz vector. Programs compiled with array bounds checking (e.g. '-fbounds-check' when using GNU-gfortran) will then stop, claiming that the index of ritz (6th line above) is out of range. With the current fix, np cannot be zero, because nev0 < n = kplusp, so np > 0 (if I'm not wrong kplusp=n, n being the matrix size; and the constraint NEV < N is clearly stated in the doc of dsaupd)
Comment 3 by Edouard Canot, Oct 24, 2013
Why the three attached files in the main report disappeared? I'm sure that there was there one hour ago, at the creation of the report. A bug in forge.scilab.org? Find them here again
Comment 4 by Edouard Canot, Oct 24, 2013
And pretty sure that the current bug affects also the pdsaup2 routine in PARPACK, which is based of the same algorihtm...
Comment 5 by Sylvestre Ledru, Nov 6, 2013
For the internal arpack-ng test system, do you think you could write a fortran test case ?
Comment 6 by Sylvestre Ledru, Nov 6, 2013
Fixed by commit 30ee9d7c3d09f5da67ad2bb435cc8ad4695b66c9 Leaving it open until we have a unitary test
Comment 7 by vedant yadav, Aug 7, 2019
Hello guys now come here and visit this site for the work because this https://abc2xyz.co/sdfghjkl site help you for the gaming here in game.
Comment 8 by suman desai, Aug 19, 2019