From help-octave-request at bevo dot che dot wisc dot edu Fri Nov 28 14:57:30 2003
Subject: Re: strange problem
From: Mike Miller <mbmiller at taxa dot epi dot umn dot edu>
To: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
cc: Gerald Ebberink <ebberink at kvi dot nl>, help-octave@bevo.che.wisc.edu
Date: Fri, 28 Nov 2003 14:57:16 -0600 (CST)

I just want to say that John's message was absolutely amazing.  What a
numerical tour de force!  Every student in computer science or numerical
analysis should read it.  The seemingly simple computation has a vast
complexity behind it.  I'm very impressed.

Mike


On Fri, 28 Nov 2003, John W. Eaton wrote:

[snip]
> There has been a lot of guessing in this thread, which has probably
> not helped too much.  Since the source code of Octave is available,
> there is no need to guess -- we can look and see precisely what is
> happening.
[snip]
> | I've used
> | GNU Octave, version 2.1.49 (i686-pc-linux-gnu).
>                                    ^^
> I think this detail turns out to be important, because it is possible
> for floating point arithmetic to be done in 80-bit registers on these
> systems, and the quick answer is that this and smart compilers is what
> is causing the problem you see.
[snip]
> Octave uses dgemm from the blas to compute matrix multiplications.
[snip]
> I can verify the result that you received with Octave using the
> following small Fortran program and the Fortran version of dgemm:
[snip]
> so the problem (if there is one) is really in dgemm, not Octave.
[snip]
> Some hardware like the 8087 and numerical processors that followed
> have 80-bit extended precision registers and can do calculations using
> this 80-bit format.  The extended precision is usually helpful, but
> sometimes it can cause trouble.
[snip]
> Octave calls dgemm from the BLAS, which does all the computations in
> Fortran (or perhaps C or assembly if you are using ATLAS or some other
> vendor-supplied version of the BLAS) and which has the opportunity to
> use extended precision.
>
> So, in the calculation that resulted in
>
>  c( 1, 1) =   1. +  -0.01 *   100.
>     =  -2.08166817E-17
>
> the value -0.01 is only displayed this way; it's value is really not
> precisely -0.01, and when it is converted to extended precision and
> multiplied by 100, the result is not precisely -1.  So adding it to 1
> (again, in extended precision mode) produces a result that is not zero,
> and it happens that the difference shows up not only in the extra
> bits, but also in one of the mantissa bits that remain after
> converting back to the 64-bit representation that is eventually
> returned to Octave.
>
> Sometimes a compiler can be smart and eliminate temporaries, and get
> extended precision even if you explicitly write a series of
> calculations like the above.
>
> GCC (including g77 and g++) have an option -ffloat-store to force
> intermediate values to be stored in memory rather than extended
> precision registers, but this only affects variables that are
> explicitly stored in temporaries.  So if you write
>
>   r = a + b*c
>
> the intermediate result of b*c may still be stored in an extended
> precision register.
>
> Unfortunately, -ffloat-store can significantly reduce the speed of
> your computations because it forces data to be moved from on-chip
> registers to memory.  (It would be nice if you could simply request
> that the hardware not use the extended precision, but as far as I
> know, that is not possible on most if not all 8087-like processors).
>
> Compiling the Fortran code above on an x86 system with -O
> -ffloat-store doesn't change things.  But rearranging the innermost
> loop slightly to be
>
>             DO 70, I = 1, M
>                temp = B( L, J )*A( I, L )
>                C( I, J ) = C( I, J ) + temp
>    70       CONTINUE
>
> *and* compiling with -O -ffloat-store (or no optimization at all, but
> that is probably not what you really want to do), produces the results
> you expect:
>
>    0.  100.
>   -0.01  1.
>
> But we probably can't get everyone to agree that this is a good
> solution, and even if you do this, you will probably still find some
> other calculations that produce results that don't match precisely
> what you would get with exact decimal arithmetic.
>
> jwe



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

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