[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 bugoctave mailing list. The
problem is really in the BLAS subroutine DGEMM. I'm also posting this
to the helpoctave mailing list because there has been some discussion
there about using vendorspecific 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 6Feb2000, 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 level3 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 Lth column of A.
What do the vendorspecific 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

 NaN problem,
John W. Eaton <=