help-octave
[Top][All Lists]

## NaN problem

 From: John W. Eaton Subject: NaN problem Date: Sun, 6 Feb 2000 23:02:54 -0600 (CST)

```[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
-----------------------------------------------------------------------

```