[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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

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?


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

Octave's home on the web:
How to fund new projects:
Subscription information:

reply via email to

[Prev in Thread] Current Thread [Next in Thread]