LAPACK Archives

[Lapack] Doubtful programing style leads to numerical inconsistency

Hi colleagues,

Working with LAPACK texts we've found some examples of programing style 
that is considered doubtful from the very beginning of hardware-assisted 
numerical computations.

One of examples is sger.f module (www.netlib.org/blas/sger.f); affected 
portion of code is given below.

IF (Y( JY ).NE.ZERO) THEN
  TEMP = ALPHA*Y( JY )
  DO 10, I = 1, M
    A( I, J ) = A( I, J ) + X( I )*TEMP
  10 CONTINUE
END IF

We see checking STRICT EQUALITY of two real numbers. It is not good. For 
example, if intermediate numerical computations are performed with 
redundant accuracy (as it takes place on Intel x86 and compatible 
processors using 80-bit floating point registers for numerical 
computations within non-SSE command subset), logic of execution changes 
if, for example,

A( I, J ) = A( I, J ) + X( I )*TEMP

is replaced by

TMP=X( I )*TEMP
A( I, J ) = A( I, J )+TMP

Second portion of code behaves much more likely to "correct" computers 
(that use 4-byte hardware FP arithmetics): a situation takes place in 
our tests, for example, when on i386 computer with the second variant we 
obtain strict zero for Y(JY), and for the first variant Y(JY) is about 
1E-18 that leads the process to another branch, and output results of 
the whole procedure are SIGNIFICANTLY DIFFERENT, being absolutely 
incorrect for unsplitted code line.

The question is not - to decompose arithmetic expressions or not, the 
question is - to compare reals without any tolerance level or not. If 
yes, the versatility of sources would be considerably affected. So, 
there is a choice - to strictly follow LAPACK sources that are de facto 
stantdard for numerical compurarions and  to have different numerical 
results for  different hardware platrforms - or to stop following LAPACK 
sources and to obtain uniform numerical results.

Ay, there is one more way out - to implement changes in official LAPACK 
sources if issues are found.

We have found a number of issues of such kind (sger.f mentioned above, 
?ggqrf and some others).

Is it possible to introduce changes into current LAPACK sources to fix 
these issues? It does not deem to be a hreat number of them. If yes, we 
can prepare proposals on how to do it and perform comprehensive 
numerical testing.

You are welcome to use my email address to communicate with us.

Sincerely yours,

Igor A. Kondrashkin,

Perflib QA Group,

Sun Microsystems, Inc.

<Prev in Thread] Current Thread [Next in Thread>


For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or