## DSYEVR hangs in DLARRB on certain matrices

Open discussion regarding features, bugs, issues, vendors, etc.

### DSYEVR hangs in DLARRB on certain matrices

There seems to be a problem in DLARRB (and I can imagine, any xLARRB), which causes it to hang on certain matrices. We're doing matrix eigen-decomposition using DSYEVR on RedHat Enterprise Linux 4, Intel Xeon x 2 using gcc 3.3.3 and the latest CLAPACK (and there were no updates to DLARRB code in the latest Fortran version). DLARRB contains this loop for unconverged intervals:
Code: Select all
`**     While( NCNVRG.LT.NEIG )*      INITI1 = I1      INITI2 = I2   80 CONTINUE      IF( NCNVRG.LT.NEIG ) THEN<...skipped...>  100    CONTINUE         NINT = NINT - OLNINT + I2         GO TO 80      END IF`

Inside this loop an interval gap is calculated like this:
Code: Select all
`                  IF( I.EQ.IFIRST ) THEN                     GAP = WERR( I+1 ) - W( I )                  ELSE IF( I.EQ.ILAST ) THEN                     GAP = WERR( I ) - W( I-1 )                  ELSE                     GAP = MIN( WERR( I+1 )-W( I ), WERR( I )-W( I-1 ) )                  END IF`

I tried to reverse engineer logic behind this function (as I don't know if there is more documentation to it), and this code looks wrong to me. Correct me if I misunderstood anything:
* for eigenvalues that converged before the loop above (IWORK[i] == 1), W[i] contains eigenvalues, and WERR[i] contains half a length of uncertainty interval
* for eigenvalues that did not converge before the loop above (IWORK[i] == 0), W[i] and WERR[i] are reused, and now they mean right and left boundaries of uncertainty interval
Given this, gap calculations above are NOT correct if the gap is referring to (i-1)-th and/or (i+1)-th eigenvalue data, and such an eigenvalue happen to converge before the loop above! In this case, W[i-1] and/or WERR[i+1] have different meaning.
Actually, in many cases it will still work (althrough not entirely correct), but for a certain matrix we've got:
Code: Select all
`werr=2.002852E-15 w=4.223298E-15 // not convergedwerr=4.485400E-16 w=3.392842E-13 // pre-converged`

See the difference in semantics if WERR and W in these 2 cases. But when gap is calculated according to the formulas above, it gets negative (werr < w). Width is always non-negative, and even though it converges to 0, the eigenvalue convergence condition
Code: Select all
`IF( WIDTH.LT.THRESH*GAP ) THEN`

is never fulfilled because the gap in this case is negative, so the function was never able to converge this last eigenvalue, so it was just stuck in this loop. Our fix (sorry, it's in C, for CLAPACK) was to put this before the loop:
Code: Select all
`    initi2 = i2; //       INITI2 = I2    if (initi1 > *ifirst) {        --initi1;        tmp = werr[initi1];        werr[initi1] = w[initi1] - werr[initi1];        w[initi1] = w[initi1] + tmp;    }    if (initi2 < *ilast) {        ++initi2;        tmp = werr[initi2];        werr[initi2] = w[initi2] - werr[initi2];        w[initi2] = w[initi2] + tmp;    }L80: //   80 CONTINUE    if (ncnvrg < neig) { //       IF( NCNVRG.LT.NEIG ) THEN`

This converts W and WERR for eigenvalues adjacent to unconverged cluster to the same format, and it is converted back together with the rest of W and WERR for unconverged values in the last loop of the function. After this, everything works, and for our matrix corrected values are:
Code: Select all
`werr=2.002852E-15 w=4.223298E-15werr=3.388356E-13 w=...`

Now, w < werr, and DLARRB finally converges.

/Victor
Victor Mozgin

Posts: 1
Joined: Tue Jun 27, 2006 3:46 pm

It seems like a bug to me. Is it going to be patched in the future?
Michael Agranat

Posts: 1
Joined: Thu Aug 10, 2006 11:22 am 