Issue could be reproduced when GEMV and GEMM routines compute different values when order of multiplications is different because of accumulated errors.

The issue was reproduced with MKL BLAS.

Current algorithm could choose zero pivot in the case when non-zero pivot could be used.

Below is tracing of the case (by lines) when the issue could be reproduced (CLAHEF_ROOK, UPLO = 'L'):

In the brackets are lines from clahef_rook.f file (revision 1598).

[753] ABSAKK == 0

[761] COLMAX > 0; COLMAX = W(IMAX, K) (value A in the matrix below)

[793] ABSAKK < COLMAX

[815] Compute W(K, K + 1) == B by CGEMV call. In fact computation of W(K, K + 1) value could be skipped because it should be equal to CONJG(A) == CONJG(W(IMAX, K)) - it can be seen from picture below because a matrix is symmetric/hermitian. But it turns out to be zero value: W(K, K + 1) == 0. This behavior is caused by rounding error in GEMM or GEMV.

[828] ROWMAX == 0, W(IMAX, K+1) == 0

[847] Case(2): W(IMAX, K+1) >= ROWMAX condition is true because ROWMAX == 0

[853] Use zero pivot without testing Case(4). It caused getting NaNs(0 / 0) in line 992.

- Code: Select all
`*`

**

***

***D B

***a* b

***a**b

IMAX ***A**b

***a**b*

K

One of the options to fix this issue is to replace value "B" by already computed value "A" (or CONJG(A) for complex case). This fix allows to satisfy condition ROWMAX >= B > 0 when A > 0 (because B will be equal to A > 0). In this case Case(2) [847] will be skipped when W(IMAX, K+1) == 0 and pivot will be computed correctly.

This kind of fix was done for all 6 routines LAHEF_ROOK and LASYF_ROOK for Lower and Upper cases. Fix is attached.

It's a little more complicated than just assignment after GEMV call because it's necessary to update index value for Case(4) when pivot wasn't found because of following column swapping.