Home > Archive > Fortran > July 2004 > Precision issue?
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
|
|
| Richard Felkins 2004-07-16, 3:58 pm |
| <rich@compt> 4% cat h.f
real*4 y,yy
yy=1.10
y=0.0
do j=1,1000000
y=y+yy
enddo
write(*,*)' single precision y = ',y
stop
end
<rich@comp> 5% ./h
single precision y = 1110920.
<rich@comp> 6%
The answer should be 1110000
Here is a simple computation provided to me to help in
isolating the cause in this precision drift. The answer
using 'real' drifts, after a w , this variance is off the
scale. Previous f77 and f90 Fortran distribution produced
the expected results. The latest upgrade Dec Alpha and the
newest Fortran from HP produces these incorrect results.
I might also note the Sun Forte7 compiler has issues
along these lines also. I know this is a form of under-
flow in a sense, but I really need a work around suggestion.
There are many, many lines of code involved here.
Thanks,
Rich.
| |
| Richard Maine 2004-07-16, 3:58 pm |
| Richard Felkins <richardf@san.rr.com> writes:
[code adds single precision 1.1 to itself a million times.]
> single precision y = 1110920.
....
> The answer should be 1110000
No, that's not what the answer "should" be. You appear to
be expecting exact answers from floating point airthmetic;
you won't get it. If you happen to get exact answers in
some special cases, those are special cases instead of
general ones. Floating point arithmetic on all current
machines in all languages is done in binary. The number
1.1 does not have an exact representation in binary in the
first place, even before you start talking about arithmetic
with it.
The resolution of single precision floating point on most
compilers is about 6 to 7 decimal digits. Your real*4
is explicit about specifying that resolution (in a nonstandard,
but widely used way). When you start adding a million things
with finite resolution together, you'll loose some precision.
All in all, the results you are seeing look well within the
bounds of reasonable expectation for the resolution you are
using.
You could use double precision to improve your results. That
may be good enough for your problem. (But be sure to change
constants such as the 1.1 to double precision forms as well).
However, even double precision is finite. I'd say that you need to
get a little better understanding about floating point arithmetic. If
you don't have at least a basic understanding of it, you will
eventually get in trouble with any finite precision, no matter how
large.
There is a paper with a title something like "What Every
Programmer Should Know About Fooating Point Arithmetic" (that title
isn't necessarily exact, but it is something like that), which
is cited regularly in this newsgroup. It might be a good choice.
> Previous f77 and f90 Fortran distribution produced
> the expected results. The latest upgrade Dec Alpha and the
> newest Fortran from HP produces these incorrect results.
Other than the correction that the new results are not
incorrect, I'll note thatthere are many reasons why you could
get different results from different compiler versions for
code like this. Data might be temporarily stored and operated
on in higher-precision registers; computations can be re-ordered
(though you don't have much computation to re-order here).
Just because you got the results that you expected in some cases
does not mean that other results are wrong or that your method
of getting the results is robust.
> I really need a work around suggestion.
> There are many, many lines of code involved here.
My best guess is that you need double precision, but without studying
the code in detail (which might take a lot of time if it is many, many
lines of code) it is hard for me to say whether that will really solve
your problems or not.
--
Richard Maine
email: my last name at domain
domain: summertriangle dot net
| |
| Dr Ivan D. Reid 2004-07-16, 8:58 pm |
| On Fri, 16 Jul 2004 08:34:44 -0700, Richard Felkins <richardf@san.rr.com>
wrote in <40F7F594.7050106@san.rr.com>:
><rich@compt> 4% cat h.f
> real*4 y,yy
> yy=1.10
> y=0.0
> do j=1,1000000
> y=y+yy
> enddo
> write(*,*)' single precision y = ',y
> stop
> end
><rich@comp> 5% ./h
> single precision y = 1110920.
><rich@comp> 6%
> The answer should be 1110000
1100000, Shirley?
> Here is a simple computation provided to me to help in
> isolating the cause in this precision drift. The answer
> using 'real' drifts, after a w, this variance is off the
> scale. Previous f77 and f90 Fortran distribution produced
> the expected results. The latest upgrade Dec Alpha and the
> newest Fortran from HP produces these incorrect results.
Are you sure that your previous installations didn't declare default
switches to promote reals to double precision?
> I might also note the Sun Forte7 compiler has issues
> along these lines also. I know this is a form of under-
> flow in a sense, but I really need a work around suggestion.
> There are many, many lines of code involved here.
As Richard has already said, it's not so much underflow as the
iherent imprecision in *binary* floating-point representation of decimal
fractions.
As examples, take these extended versions of your programme (g77
under cygwin):
$ cat tsp.f
implicit none
integer i,j,k
real a,y(0:6)
a=1.1
do j=0,6
k=10**j
y(j)=0.
do i=1,k
y(j)=y(j)+a
enddo
enddo
write(*,*) y
stop
end
$ ./tsp
1.10000002 11.000001 109.999901 1099.98901 10999.5156 110101.477
1110920.5
$ cat tspd.f
implicit none
integer i,j,k
double precision a,y(0:6)
a=1.1d0
do j=0,6
k=10**j
y(j)=0.0d0
do i=1,k
y(j)=y(j)+a
enddo
enddo
write(*,*) y
write(*,'(f20.10)')y
stop
end
$ ./tspd
1.1 11. 110. 1100. 11000. 110000. 1100000.
1.1000000000
11.0000000000
110.0000000000
1100.0000000000
11000.0000000020
110000.0000001741
1099999.9999886872
....and, for completeness:
$ cat tspds.f
implicit none
integer i,j,k
double precision a,y(0:6)
a=1.1e0
do j=0,6
k=10**j
y(j)=0.0d0
do i=1,k
y(j)=y(j)+a
enddo
enddo
write(*,*) y
write(*,'(f20.10)')y
stop
end
$ ./tspds
1.10000002 11.0000002 110.000002 1100.00002 11000.0002 110000.002
1100000.02
1.1000000238
11.0000002384
110.0000023842
1100.0000238419
11000.0002384186
110000.0023841858
1100000.0238418579
[Note that g77 gives identical results here if a is defined REAL rather
than DOUBLE PRECISION, which implies that SP storage into DP is zero-filled.]
--
Ivan Reid, Electronic & Computer Engineering, ___ CMS Collaboration,
Brunel University. Ivan.Reid@brunel.ac.uk Room 40-1-B12, CERN
KotPT -- "for stupidity above and beyond the call of duty".
| |
| Eric K. 2004-07-16, 8:58 pm |
| Richard Felkins <richardf@san.rr.com> wrote in message news:<40F7F594.7050106@san.rr.com>...
> real*4 y,yy
> yy=1.10
> y=0.0
> do j=1,1000000
> y=y+yy
> enddo
> write(*,*)' single precision y = ',y
> stop
> end
> <rich@comp> 5% ./h
> single precision y = 1110920.
> <rich@comp> 6%
> I know this is a form of under-
> flow in a sense, but I really need a work around suggestion.
> There are many, many lines of code involved here.
One workaround that will mitigate, though not eliminate, the
roundoff issue is to use an expression like j*yy instead of the
intermediate sum y wherever it is needed. You can reduce (but
not eliminate) the roundoff error even more with "real*8 yy";
just be sure to initialize it as yy=1.10d0, NOT yy=1.10.
Another possibility is to scale the problem so that all intermediate
results are stored as sufficiently wide integers (integer*8, for instance).
Then convert to floating-point values only when needed. Even here,
however, some rounding error is likely in the conversion to floating point.
Quantities like 1.1 are not exactly representable in binary arithmetic,
so the last bits in the value of j*yy wobble a bit depending on the
value of j. You can't expect an "exact" result in binary floating-point
arithmetic.
| |
| Michael Prager 2004-07-19, 4:00 pm |
| Richard Maine <nospam@see.signature> wrote:
[...]
>
>There is a paper with a title something like "What Every
>Programmer Should Know About Fooating Point Arithmetic" (that title
>isn't necessarily exact, but it is something like that), which
>is cited regularly in this newsgroup. It might be a good choice.
The Perils of Floating Point, by Bruce Bush:
http://www.lahey.com/float.htm
It outlines the general issues (though it seems that's been done
by other posters as well.).
--
Mike Prager, NOAA, Beaufort, NC
Address spam-trapped; remove color to reply.
* Opinions expressed are personal and not represented otherwise.
* Any use of tradenames does not constitute a NOAA endorsement.
| |
| Paul Van Delst 2004-07-19, 4:00 pm |
| Richard Felkins wrote:
> <rich@compt> 4% cat h.f
>
> real*4 y,yy
> yy=1.10
> y=0.0
> do j=1,1000000
> y=y+yy
> enddo
> write(*,*)' single precision y = ',y
> stop
> end
> <rich@comp> 5% ./h
> single precision y = 1110920.
> <rich@comp> 6%
>
> The answer should be 1110000
Others have replied why the answer isn't what you think it should be. The way the code you
posted is accumulating is simply a bad way of adding numbers because it doesn't take the
imprecision of floating point calculations into account.
I realise your example was probably oversimplified for posting to clf, but can you do
something like:
y = real(i)*yy
and avoid the summation altogether?
If you are adding a set of numbers, e.g.
y = 0.0
do j = 1, 1000000
y = y + x(j)
end do
you may want to look into some sort of compensated summation algorithm. Kahan's method is
a commonly known one, but there are plenty of others (most of them with their own caveats.)
Apart from all that, double precision will reduce the "precision drift" of the actual
answer from the expected answer to a very small fraction of 1.0. (If you print the result
with the "right" format string, you will get the exactly the expected answer printed. :o)
cheers,
paulv
>
> Here is a simple computation provided to me to help in
> isolating the cause in this precision drift. The answer
> using 'real' drifts, after a w, this variance is off the
> scale. Previous f77 and f90 Fortran distribution produced
> the expected results. The latest upgrade Dec Alpha and the
> newest Fortran from HP produces these incorrect results.
>
> I might also note the Sun Forte7 compiler has issues
> along these lines also. I know this is a form of under-
> flow in a sense, but I really need a work around suggestion.
> There are many, many lines of code involved here.
>
> Thanks,
>
> Rich.
>
>
| |
| beliavsky@aol.com 2004-07-19, 4:00 pm |
|
Paul Van Delst <paul.vandelst@noaa.gov> wrote:
>Others have replied why the answer isn't what you think it should be. The
way the code
>you
>posted is accumulating is simply a bad way of adding numbers because it
doesn't take
>the
>imprecision of floating point calculations into account.
>
>I realise your example was probably oversimplified for posting to clf, but
can you do
>
>something like:
>
> y = real(i)*yy
>
>and avoid the summation altogether?
>
>If you are adding a set of numbers, e.g.
>
> y = 0.0
> do j = 1, 1000000
> y = y + x(j)
> end do
>
>you may want to look into some sort of compensated summation algorithm.
Kahan's method
>is
>a commonly known one, but there are plenty of others (most of them with
their own caveats.)
In Fortran 90/95 one can write y = SUM(x). I wonder how sophisticated the
implementation of SUM is in Fortran compilers. Do they just use the loop
above, perhaps with loop unrolling?
To reduce roundoff error, one should sum the numbers in order of increasing
size if they are all positive. Sometimes one has information about the magnitudes
and ordering of data. I wonder if it would make sense to add optional arguments
to the SUM function and other intrinsic functions, so that one could write
y = SUM(x,order="descending",sign="positive")
and have the compiler use the optimal summation algorithm.
Optional arguments could also help a compiler optimize. If one had a statement
if (sum(x,sign="positive") > y) do_something
The compiler could stop evaluating the recursive sum as soon as it exceeded
y.
Can someone suggest some easily generated sequences of floating point numbers
for which the difference between the true sum and that generated by naive
recursive summation is large? I want to see what compilers do.
----== Posted via Newsfeed.Com - Unlimited-Uncensored-Secure Usenet News==----
http://www.newsfeed.com The #1 Newsgroup Service in the World! >100,000 Newsgroups
---= 19 East/West-Coast Specialized Servers - Total Privacy via Encryption =---
| |
| Paul Van Delst 2004-07-19, 4:00 pm |
| beliavsky@aol.com wrote:
> Paul Van Delst <paul.vandelst@noaa.gov> wrote:
>
>
>
> way the code
>
>
> doesn't take
>
>
> can you do
>
>
> Kahan's method
>
>
> their own caveats.)
>
> In Fortran 90/95 one can write y = SUM(x). I wonder how sophisticated the
> implementation of SUM is in Fortran compilers. Do they just use the loop
> above, perhaps with loop unrolling?
If a loop is generated, unrolling it appropriately would be seemly I would think.
>
> To reduce roundoff error, one should sum the numbers in order of increasing
> size if they are all positive.
Yeah, the sorting of data to improve some summation algorithms was pretty much what I
meant by "caveat". The sorts of data I typically sum (1000000's of points of spectral
atmospheric transmittances) are always 0<x<1 so sorting doesn't provide a big enough
improvement to balance the actual sorting (in both time and code complexity). Have I
mentioned previously that I do pretty simple stuff? :o)
> Sometimes one has information about the magnitudes
> and ordering of data. I wonder if it would make sense to add optional arguments
> to the SUM function and other intrinsic functions, so that one could write
>
> y = SUM(x,order="descending",sign="positive")
>
> and have the compiler use the optimal summation algorithm.
Hmm. I don't know about that. I sort of like having to think about that sort of stuff
myself. Not because I'm a masochist, but it makes me appreciate and cogitate about
numerical precision/problems for subsequent problems.
cheers,
paulv
| |
| Richard Maine 2004-07-19, 4:00 pm |
| "beliavsky@aol.com" <beliavsky@127.0.0.1:7501> writes:
> To reduce roundoff error, one should sum the numbers in order of increasing
> size if.... [and other discussion]
Of course, all this assumes that the summation is done in isolation as
an end to itself. Hard to tell for sure from the simplified sample
shown, but often there are other things going on in the same loop
as the summation, with the summation itself a secondary issue. A
prime example is something like solving a differential equation; in
cases such as that, re-ordering the values is not likely to be
even sensible.
Paul's suggestion to substitute multiplication in place of addition
is often (but not always) a good approach there. Hard to say without
a bigger picture of the real application.
--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain
| |
| glen herrmannsfeldt 2004-07-19, 8:57 pm |
| Richard Felkins wrote:
> real*4 y,yy
> yy=1.10
> y=0.0
> do j=1,1000000
> y=y+yy
> enddo
> write(*,*)' single precision y = ',y
> stop
> end
> single precision y = 1110920.
> The answer should be 1110000
The exact answer should be 1100000 not 1110000, but other
than that, it seems to be a case where rounding errors
accumulate instead of cancel on average.
One solution might be probabilistic rounding, where rounding
is done with the help of a random number generator. I don't
know of any implementations, though. Normally when adding
different numbers rounding errors don't accumulate as they
do in this case.
-- glen
| |
| Richard Felkins 2004-07-19, 8:57 pm |
|
Thanks for all that have responded. You have provided
several nice alternatives to be evaluated, I really
do appreciate your contributions.
Richard.
Richard Maine wrote:
> "beliavsky@aol.com" <beliavsky@127.0.0.1:7501> writes:
>
>
>
>
> Of course, all this assumes that the summation is done in isolation as
> an end to itself. Hard to tell for sure from the simplified sample
> shown, but often there are other things going on in the same loop
> as the summation, with the summation itself a secondary issue. A
> prime example is something like solving a differential equation; in
> cases such as that, re-ordering the values is not likely to be
> even sensible.
>
> Paul's suggestion to substitute multiplication in place of addition
> is often (but not always) a good approach there. Hard to say without
> a bigger picture of the real application.
>
| |
| glen herrmannsfeldt 2004-07-19, 8:57 pm |
| Richard Felkins wrote:
>
> Thanks for all that have responded. You have provided
> several nice alternatives to be evaluated, I really
> do appreciate your contributions.
I tried a few more tests on this, to be sure that
I understood the problem.
Note that the problem is not due to rounding in the binary
representation of 1.1, as one might expect. The low bit is
only important for the first few additions. As y gets larger,
successively more significant bits of yy become important to
the rounding. 1.1 as a binary fraction to 24 bits is:
1.00011001100110011001101
Most of the effect comes in about the last half of the additions,
where only the most significant bits are important. If you try
1.00011, decimal 1.0975, the sum is only slightly smaller.
-- glen
| |
| Tim Prince 2004-07-20, 3:57 am |
|
"beliavsky@aol.com" <beliavsky@127.0.0.1:7501> wrote in message
news:40fc046a$1_1@127.0.0.1...
> In Fortran 90/95 one can write y = SUM(x). I wonder how sophisticated the
> implementation of SUM is in Fortran compilers. Do they just use the loop
> above, perhaps with loop unrolling?
>
8 parallel sums is a common scheme. This often produces better accuracy
than a single sequential sum without extra precision, but "maintain
precision" options may disable this parallel treatment.
> To reduce roundoff error, one should sum the numbers in order of
increasing
> size if they are all positive. Sometimes one has information about the
magnitudes
> and ordering of data. I wonder if it would make sense to add optional
arguments
> to the SUM function and other intrinsic functions, so that one could write
>
> y = SUM(x,order="descending",sign="positive")
>
> and have the compiler use the optimal summation algorithm.
Summing with extra precision is far more effective, when applicable.
The compiler is supposed to divine what you consider optimal?
I've been bewildered lately by the differences in results from various
compilers with the C++ STL partial_sum() operator. An equally bad, probably
slower, but self-explanatory, Fortran version:
forall(i=1:n)a(i)=sum(b(:i))
(Unless your figure of merit is the number of redundant operations
performed)
The accuracy problem could be solved by promoting sum() to the next bigger
KIND than the default (same KIND as b(:)). What's the most succinct and
portable way to express that?
| |
| Gerry Thomas 2004-07-20, 3:57 am |
|
"Tim Prince" <tprince@nospamcomputer.org> wrote in message
news:YQ%Kc.94130$3K5.67908@newssvr29.news.prodigy.com...
> The accuracy problem could be solved by promoting sum() to the next
bigger
> KIND than the default (same KIND as b(:)).
Unmitigated nonsense!
> What's the most succinct and
> portable way to express that?
Don't!
--
You're Welcome,
Gerry T.
______
"Don't play dumb! You're not as good at it as I am!" -- Col Sam Flagg,
ICORPS dropin to the 4077th M*A*S*H
| |
| Arjen Markus 2004-07-20, 3:57 am |
| "beliavsky@aol.com" wrote:
>
>
>
> Can someone suggest some easily generated sequences of floating point numbers
> for which the difference between the true sum and that generated by naive
> recursive summation is large? I want to see what compilers do.
>
>
I do not remember who posted this code, but here is a program that
appeared
in this group about a year ago.
Regards,
Arjen
------
! Program to analyse the accuracy of summing
!
! See: Y. He and C. Ding: Using accurate arithmetics to improve
! numerical reproducibility and stability in parallel applications
!
program acc_sum
implicit none
integer, parameter :: wp = kind(1.0d0)
real (kind = wp), dimension(1:100,1:127) :: fvalue
real (kind = wp) :: sum_i_first
real (kind = wp) :: sum_j_first
real (kind = wp) :: sum_via_sum
real (kind = wp) :: expected_sum
integer :: i
integer :: j
call fill_array( fvalue, expected_sum )
sum_i_first = 0.0_wp
sum_j_first = 0.0_wp
do j = 1,size(fvalue,2)
do i = 1,size(fvalue,1)
sum_i_first = sum_i_first + fvalue(i,j)
enddo
enddo
do i = 1,size(fvalue,1)
do j = 1,size(fvalue,2)
sum_j_first = sum_j_first + fvalue(i,j)
enddo
enddo
sum_via_sum = sum( fvalue )
write(*,*) 'Expected sum ', expected_sum
write(*,*) 'Sum: i first ', sum_i_first
write(*,*) 'Sum: j first ', sum_j_first
write(*,*) 'Sum: via sum ', sum_via_sum
stop
contains
subroutine fill_array( fvalue, expected_sum )
real (kind=wp), dimension(:,:) :: fvalue
real (kind=wp) :: alpha = 0.05_wp
real (kind=wp) :: beta = 0.01_wp
real (kind=wp) :: pi = 3.1415926_wp
real (kind = wp) :: expected_sum
integer :: i
integer :: j
integer :: m
m = size(fvalue,1)
expected_sum = beta * size(fvalue)
do j = 1,size(fvalue,2)
do i = 1,size(fvalue,1)
fvalue(i,j) = sin( pi * (2.0_wp*i-m-1)/2.0_wp/(m-1) ) * &
exp( alpha*j) + beta
enddo
enddo
end subroutine fill_array
end program acc_sum
| |
| Salvatore Filippone 2004-07-20, 8:58 am |
| Richard Felkins <richardf@san.rr.com> wrote in message news:<40F7F594.7050106@san.rr.com>...
> <rich@compt> 4% cat h.f
>
> real*4 y,yy
> yy=1.10
> y=0.0
> do j=1,1000000
> y=y+yy
> enddo
> write(*,*)' single precision y = ',y
> stop
> end
> <rich@comp> 5% ./h
> single precision y = 1110920.
> <rich@comp> 6%
>
> The answer should be 1110000
.................
> Thanks,
>
> Rich.
This is a well known issue, so well known that I use it for students
who need an introduction to floating-point computations.
To get all the info you need (and much more) you should consult the
book:
N. Higham
Accuracy and stability of numerical algorithms
SIAM
It has an entire chapter devoted to the problem of summation.
Cheers
Salvatore
|
|
|
|
|