Issue 1329: dsaupd and 'BE' option returns wrong eigenvalues for a SPD matrix

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 

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 
  $ 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 
2. Ask for some eigenvalues (2, 3 or 4) with the 'BE' option

Expected result:

Actual result:

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 

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 

A bug in

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 site help you for the gaming here 
in game.

Created: 6 years 10 months ago by Edouard Canot

Updated: 1 year 1 month ago

Status: New

Followed by: 3 persons