[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: strange problem
From: |
Mike Miller |
Subject: |
Re: strange problem |
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
-------------------------------------------------------------
- Re: strange problem, (continued)
Re: strange problem, Laurent Jacques, 2003/11/28
Message not available
strange problem, John W. Eaton, 2003/11/28
- Re: strange problem,
Mike Miller <=
RE: strange problem, THOMAS Paul Richard, 2003/11/28
RE: strange problem, THOMAS Paul Richard, 2003/11/28
Re: strange problem, pkienzle, 2003/11/29