NaN problem

[John W. Eaton]

[The following bug was reported to the bug-octave mailing list.  The
problem is really in the BLAS subroutine DGEMM.  I'm also posting this
to the help-octave mailing list because there has been some discussion
there about using vendor-specific or other optimized versions of the
BLAS with Octave.  I can fix the problem described below by modifying
the Fortran BLAS that is distrubuted with Octave, but I can't fix it
for other versions of the BLAS.  --jwe]


On  6-Feb-2000, Mirek Kwasniak <address@hidden> wrote:

| octave:72> [0;0]*[ inf nan ]
| ans =
| 
|            NaN           NaN
|            NaN           NaN
|
| octave:70> [ inf; nan ]*[0 0]
| ans =
| 
|   0  0
|   0  0
| 
| :((((((((((((((((((

In case it is not obvious, the first result is correct; the second is
not.

Octave uses the level-3 blas routine DGEMM for this calculation.
DGEMM uses the following algorithm:

*           Form  C := alpha*A*B + beta*C.
*
            DO 90, J = 1, N
               IF( BETA.EQ.ZERO )THEN
                  DO 50, I = 1, M
                     C( I, J ) = ZERO
   50             CONTINUE
               ELSE IF( BETA.NE.ONE )THEN
                  DO 60, I = 1, M
                     C( I, J ) = BETA*C( I, J )
   60             CONTINUE
               END IF
               DO 80, L = 1, K
                  IF( B( L, J ).NE.ZERO )THEN
                     TEMP = ALPHA*B( L, J )
                     DO 70, I = 1, M
                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
   70                CONTINUE
                  END IF
   80          CONTINUE
   90       CONTINUE

So, if B contains zeros, no mulitplication is done.  This optimization
is not really appropriate if there are also any Inf or NaN values in
the L-th column of A.

What do the vendor-specific versions of the BLAS do for this code?
What does ATLAS do?

jwe



-----------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.che.wisc.edu/octave/octave.html
How to fund new projects:  http://www.che.wisc.edu/octave/funding.html
Subscription information:  http://www.che.wisc.edu/octave/archive.html
-----------------------------------------------------------------------