- 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

Actually, in many cases it will still work (althrough not entirely correct), but for a certain matrix we've got:

- Code: Select all
`werr[19]=2.002852E-15 w[19]=4.223298E-15 // not converged`

werr[20]=4.485400E-16 w[20]=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[20] < w[19]). 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[19]=2.002852E-15 w[19]=4.223298E-15`

werr[20]=3.388356E-13 w[20]=...

Now, w[19] < werr[20], and DLARRB finally converges.

/Victor