LAPACK slaln2.f: more accuracy hurts if not consistent

[langou]

This is a bug report from David Hough

Subject: LAPACK slaln2.f: more accuracy hurts if not consistent
Date: Wed, 29 Oct 2014 15:29:20 -0700 (PDT)

Another example turned up by enabling fused FMA operations to replace
multiply/add pairs at the compiler's discretion.

slaln2.f solves 1x1 or 2x2 systems of equations. Couldn't be simpler, right?

http://www.netlib.org/lapack/explore-html/da/d2d/slaln2_8f.html

But its invocation from the TESTING programs schkee.f/schkec.f/sget31.f
and the test input data file sec.in
gets into trouble with too-high residuals computed in sget31.f. Since slaln2.f uses total pivoting, how can it possibly have problems
with matrices of dimension <= 2?

After all, slaln2.f computes the answers with code like:

298 * Real 1x1 system.
299 *
300 * C = ca A - w D
301 *
302 csr = caa( 1, 1 ) - wrd1
...
323 x( 1, 1 ) = ( b( 1, 1 )*scale ) / csr

and

327 * Complex 1x1 system (w is complex)
328 *
329 * C = ca A - w D
330 *
331 csr = caa( 1, 1 ) - wrd1
332 csi = -wid1
...
354 CALL sladiv( scaleb( 1, 1 ), scale*b( 1, 2 ), csr, csi,
355 $ x( 1, 1 ), x( 1, 2 ) )

and sget31 checks them with code like:

221 res = abs( ( caa( 1, 1 )-wrd1 )*
222 $ x( 1, 1 )-scale*b( 1, 1 ) )

279 res = abs( ( caa( 1, 1 )-wrd1 )*
280 $ x( 1, 1 )+( wid1 )x( 1, 2 )-
281 $ scaleb( 1, 1 ) )
282 res = res + abs( ( -wid1 )x( 1, 1 )+
283 $ ( caa( 1, 1 )-wr*d1 )x( 1, 2 )-
284 $ scaleb( 1, 2 ) )

A very straighforward residual computation: to check x = bs / c compute
cx - b*s

The problems arise when higher precision is used capriciously, as when
only one of sget31 or slaln2 replaces the ca-wd with an FMA, because they
are separately compiled, perhaps with different compiler options.

Then the computed residual is nowhere near where it should be if ca-wd
mostly cancels in single precision. And cchkec.f has this to say about
that:

Error in SLALN2: RMAX = 0.357E+07
LMAX = 60625 NINFO= 0 71553 KNT= 230400

My working solution - again not canonical LAPACK because it uses higher
precision - is to promote the candidates for fusing to higher precision:
In slaln2:

< CSR = CAA( 1, 1 ) - WRD1

CSR = dble(CA)*dble(A( 1, 1 )) - dble(WR)*dble(D1)

< CSR = CAA( 1, 1 ) - WRD1

CSR = dble(CA)*dble(A( 1, 1 )) - dble(WR)*dble(D1)

In sget31:

< RES = ABS( ( CAA( 1, 1 )-WRD1 )*
< $ X( 1, 1 )-SCALE*B( 1, 1 ) )

RES = ABS( ( dble(CA)dble(A( 1, 1 ))-dble(WR)dble(D1 ))
$ (X( 1, 1 ))-(SCALE)(B( 1, 1 )) )

< RES = ABS( ( CAA( 1, 1 )-WRD1 )*
< $ X( 1, 1 )+( WID1 )X( 1, 2 )-
< $ SCALEB( 1, 1 ) )
< RES = RES + ABS( ( -WID1 )X( 1, 1 )+
< $ ( CAA( 1, 1 )-WR*D1 )X( 1, 2 )-
< $ SCALEB( 1, 2 ) )

RES = ABS( x ( dble(CA)dble(A( 1, 1 ))-dble(WR)dble(D1) ) X( 1, 1 )+
x ( (WI)(D1) )(X( 1, 2 ))- SCALEB( 1, 1 ) x )
RES = RES + ABS( $ ( dble(CA)dble(A( 1, 1 ))-dble(WR)dble(D1) )X( 1, 2 )+
x (-(WI)(D1) )(X( 1, 1 ))- SCALEB( 1, 2 ) x )

That seems to do the trick. Interestingly, in my test environment,
corresponding bugs did not come up in the [cdz] versions of the same
functions.

It's not a question of whether higher precision is used. It's a question
of whether higher precision is used consistently. In a world where fusing multiple-add is mostly advantageous, or at least harmless, and is
either on or off globally, sensitive code needs to be protected.

The right solution is to have ways of indicating locally within the program text
where higher precision or other value-changing optimizations should be
avoided. Those places are relatively few, but hard to find after the
fact. However the right solution awaits a more cooperative attitude
from language designers and implementers.