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.
|