Problems with CSYSVXX / ZSYSVXX / SSYSVXX / DSYSVXX tests

[mickpont]

I built the latest version of LAPACK, 3.4.2, and told it to use the extended precision XBLAS by editing make.inc and setting these lines:

USEXBLAS = Yes
XBLASLIB = libxblas.a

where libxblas.a is a compiled copy of XBLAS version 1.0.248 that I also got from netlib.

I used gfortran for the build (a development branch build of gfortran but I don't think that's relevant).

If I build LAPACK without the use of XBLAS, then I see no problems at all - since that's the default I imagine that's what most people do.

With an XBLAS-using LAPACK build, most tests still seem to work fine, but I found some trouble when running the test program xlintstc using ctest.in as input, getting error messages like these:

CSYSVXX, FACT='F', UPLO='U', N = 1, type 1, test 5, ratio = 0.83886E+07
CSYSVXX, FACT='N', UPLO='U', N = 1, type 1, test 5, ratio = 0.83886E+07
CSYSVXX, FACT='F', UPLO='L', N = 1, type 1, test 5, ratio = 0.83886E+07
CSYSVXX, FACT='N', UPLO='L', N = 1, type 1, test 5, ratio = 0.83886E+07

and about 150 more messages (different values of N and test number).

I tracked some of the problems down to the test routine TESTING/LIN/cdrvsyx.f. After line 603 in that file:

* --- Test CSYSVXX ---

the code goes like this:

...
CALL CSYSVXX( FACT, UPLO, N, NRHS, A, LDA, AFAC,
...
* Compute residual of the computed solution.
*
CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
CALL CSYT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
$ LDA, RWORK( 2*NRHS+1 ), RESULT( 2 ) )
RESULT( 2 ) = 0.0
*
* Check solution from generated exact solution.
*
CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
$ RESULT( 3 ) )
*
* Check the error bounds from iterative refinement.
*
CALL CPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
$ XACT, LDA, RWORK, RWORK( NRHS+1 ),
$ RESULT( 4 ) )
ELSE
K1 = 6
END IF
*
* Compare RCOND from CSYSVXX with the computed value
* in RCONDC.
*
RESULT( 6 ) = SGET06( RCOND, RCONDC )

There are a couple of problems here:

(o) First of all, notice that after the call of CSYT02, which computes RESULT( 2 ), the next line sets RESULT( 2 ) = 0.0 - this looks like a mistake because test 2 can then never fail.

(o) Then, notice that the call of CPOT05 passes array RWORK as the 12th argument. This is argument FERR of CPOT05, and is supposed to contain forward error estimates. The trouble is that forward error estimates were never placed into RWORK by the test program (actually, RWORK was used as a genuine workspace argument in the call of CSYSVXX).

Comparing this code in cdrvsyx.f with similar code in the equivalent REAL code sdrvsyx.f, I think I can see what has happened here. In sdrvsyx.f, uring the test of SSYSVX, RWORK is used as the FERR argument of SSYSVX and so does contain the forward error estimates. Test routine sdrvsyx.f then goes on to test SSYSVXX, but avoids overwriting the contents of RWORK(1:N). It then again tests that RWORK contains good forward error estimates, this time assuming that they came from SSYSVXX, when in fact they are still hanging around from the earlier call of SSYSVX. So, I think sdrvsyx.f is incorrect in this regard.

In the COMPLEX case, the code to test the forward errors returned by CSYSVXX also passes RWORK as the forward error estimates, but unlike in the REAL case, RWORK doesn't contain any estimates because the earlier call of CSYSVX didn't use RWORK to hold them. So, this explains some of the reported test failures for CSYSVXX. I think that something like this may be needed immediately before the call of CPOT05 to copy error estimates into RWORK:

do k = 1, nrhs
rwork(k) = errbnds_n(k,2)
end do

However, I'm not completely sure of the difference between errbnds_n and errbnds_c (whenever I examine the contents of the two arrays they are identical, though I see that they are defined differently in CSYSVXX). Furthermore, adding this extra code does get rid of a few of the reported errors, but most still remain, and in these cases it seems to me that the contents of errbnds_n are too optimistic and so the test program is correct to report errors. For example, in a case I examined where CSYSVX said that the forward error was around 0.3, CSYSVXX said that it was about 1.0e-7 even though the results returned in X didn't look very good and were not much different to those returned by CSYSVX. So, I think there may be something wrong with CSYSVXX itself, but I've not tried to diagnose it.

Of course, all the comments above relating to CSYSVXX and SSYSVXX apply equally to the ZSYSVXX and DSYSVXX test programs.