Fortran - error in simple loop !!

Author

error in simple loop !!

nakisa

2007-04-05, 7:11 pm

Hi everybody
I have encountered with a strange error !
Here is a code that I wrote :

program amplification
implicit none

real f
real u1, tE
real b1, theta , t0 ,t
real p
real x ,y
real xi,xf,stepx
integer i

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

pi=4.0*atan(1.0)
theta=pi/4

t0=0
tE=10
b1=0.05

xi=-6.
xf=0
stepx=0.001

i=0
t=-0.85
u1=b1*sin(theta)+((t-t0)/tE)*cos(theta)

open(1,file='test.txt')
y=0
do x=xi,xf,stepx
f=abs(u1-(x-x/(x**2+y**2)))

if (f<0.005) then
i=i+1
write (1,*),i,x,y,f

end if ! if of fnew
end do ! do of x

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

close (1)

end program amplification

Error take places in selecting initial value of loop (i.e xi) , when I
select xi=-7. , the answers are completely different with when I
select xi=-6 !!!!
I can't understand logic of program ! answers are near -1.0 , why xi
affect them ??
Can anybody say me how fortran thinks about this code ??
Best,nakisa

Beliavsky

2007-04-05, 7:11 pm

On Apr 5, 2:08 pm, "nakisa" <nakisa.noor...@gmail.com> wrote:
> Hi everybody
> I have encountered with a strange error !
> Here is a code that I wrote :

<snip>

> real x ,y

<snip>

> do x=xi,xf,stepx
> f=abs(u1-(x-x/(x**2+y**2)))
>
> if (f<0.005) then
> i=i+1
> write (1,*),i,x,y,f
>
> end if ! if of fnew
> end do ! do of x
>
> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

> close (1)
>
> end program amplification

You are using a DO loop with a non-integer loop variable, which is not
standard Fortran 95. The committee removed it from Fortran because it
is an unreliable construct.

I suggest rewriting your code using an integer loop variable,
something like

n = (xf-xi)/xi ! not sure if nint should be used here
do ix=1,n
x = (ix-1)*stepx + xi
....

> do x=xi,xf,stepx
> f=abs(u1-(x-x/(x**2+y**2)))

Beliavsky

2007-04-05, 7:11 pm

On Apr 5, 2:14 pm, "Beliavsky" <beliav...@aol.com> wrote:

<snip>

> I suggest rewriting your code using an integer loop variable,
> something like
>
> n = (xf-xi)/xi ! not sure if nint should be used here
> do ix=1,n
> x = (ix-1)*stepx + xi

I should have written that one DEFINES the grid in terms of xi, n, and
stepx. Within the loop, one could treat ix=n as a special case to set
the last grid value "exactly".

Dick Hendrickson

2007-04-05, 7:11 pm

nakisa wrote:
> Hi everybody
> I have encountered with a strange error !
> Here is a code that I wrote :
>
> program amplification
> implicit none
>
>
> real f
> real u1, tE
> real b1, theta , t0 ,t
> real p
> real x ,y
> real xi,xf,stepx
> integer i
>
> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

> pi=4.0*atan(1.0)

Are you sure this is the program you actually tried? Because of the
implicit none statement, this line should have generated an error when
you compiled the program. There is no declaration for "pi". If you
are looking at this program and running a different one it will be
very hard to understand the results.

> theta=pi/4
>
> t0=0
> tE=10
> b1=0.05
>
> xi=-6.
> xf=0
> stepx=0.001
>
> i=0
> t=-0.85
> u1=b1*sin(theta)+((t-t0)/tE)*cos(theta)
>
> open(1,file='test.txt')
> y=0
> do x=xi,xf,stepx

It isn't usually a good idea to use real variables for DO loop
indexes. Values like .001 are not exactly represented in
a computers binary arithmetic. When you ask for -6.0
-6.0 + .001, -6.0 + .001 + .001, etc. you can get a large
amount of round off error, just because the .001 is inaccurate.
If you start with xi = -7.0 one one series, it will not have the
value -6.0 after 1000 iterations. It will have some value near
-6.0, but it won;t be exact. Small differences in x might cause
large differences in f.

Dick Hendrickson
> f=abs(u1-(x-x/(x**2+y**2)))
>
> if (f<0.005) then
> i=i+1
> write (1,*),i,x,y,f
>
> end if ! if of fnew
> end do ! do of x
>
> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

> close (1)
>
> end program amplification
>
> Error take places in selecting initial value of loop (i.e xi) , when I
> select xi=-7. , the answers are completely different with when I
> select xi=-6 !!!!
> I can't understand logic of program ! answers are near -1.0 , why xi
> affect them ??
> Can anybody say me how fortran thinks about this code ??
> Best,nakisa
>

glen herrmannsfeldt

2007-04-06, 8:04 am

Dick Hendrickson wrote:

> nakisa wrote:
(snip)

> It isn't usually a good idea to use real variables for DO loop
> indexes. Values like .001 are not exactly represented in
> a computers binary arithmetic. When you ask for -6.0
> -6.0 + .001, -6.0 + .001 + .001, etc. you can get a large
> amount of round off error, just because the .001 is inaccurate.
> If you start with xi = -7.0 one one series, it will not have the
> value -6.0 after 1000 iterations. It will have some value near
> -6.0, but it won;t be exact. Small differences in x might cause
> large differences in f.

While there are many cases that should not be done with REAL
loop variables, I am not convinced that this is one.
It should especially not be done when one depends on the number
of times the loop is executed. In this case, though, yes
there will be accumulated rounding error, and yes the results
might be different because of that, but in that case the program
is too sensitive to the exact value of the variable, and will be
wrong in any case.
[color=darkred]

for non-zero values of ul or y it is not obvious at all
when the WRITE statement should be executed, and there
could be a large dependence of f on the exact value of x.

-- glen

Herman D. Knoble

2007-04-06, 8:04 am

Nakisa: Try the following slight modification of your code

module SetPrecision
implicit none
integer, public, parameter :: DP = selected_real_kind(15)
end module SetPrecision

program amplification
use SetPrecision
real(kind=DP) :: f, pi, u1, tE, b1, theta, t0, t, x, y
real(kind=DP) :: xi, xf, stepx
integer i, j, nsteps

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

pi=4.0*atan(1.0_DP)
theta=pi/4.0_DP

t0=0
tE=10
b1=0.05_DP

xi=-6.
! xi=-7.
xf=0
stepx=0.001_DP

i=0
t=-0.85_DP
u1=b1*sin(theta)+((t-t0)/tE)*cos(theta)

! open(1,file='test.txt')
y=0
! do x=xi,xf,stepx
nsteps=(xf-xi)/stepx + 0.5_DP
do j=1,nsteps
x=xi + (j-1)*stepx
f=abs(u1-(x-x/(x**2+y**2)))

if (f<0.005_DP) then
i=i+1
write (*,*) i,x,y,f

end if ! if of fnew
end do ! do of x

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! close (1)

end program amplification

Skiip Knoble

On 5 Apr 2007 11:08:02 -0700, "nakisa" <nakisa.nooraee@gmail.com> wrote:

-|Hi everybody
-|I have encountered with a strange error !
-|Here is a code that I wrote :
-|
-| program amplification
-| implicit none
-|
-|
-| real f
-| real u1, tE
-| real b1, theta , t0 ,t
-| real p
-| real x ,y
-| real xi,xf,stepx
-| integer i
-|
-| !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

-| pi=4.0*atan(1.0)
-| theta=pi/4
-|
-| t0=0
-| tE=10
-| b1=0.05
-|
-| xi=-6.
-| xf=0
-| stepx=0.001
-|
-| i=0
-| t=-0.85
-| u1=b1*sin(theta)+((t-t0)/tE)*cos(theta)
-|
-| open(1,file='test.txt')
-| y=0
-| do x=xi,xf,stepx
-| f=abs(u1-(x-x/(x**2+y**2)))
-|
-| if (f<0.005) then
-| i=i+1
-| write (1,*),i,x,y,f
-|
-| end if ! if of fnew
-| end do ! do of x
-|
-| !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

-| close (1)
-|
-| end program amplification
-|
-|Error take places in selecting initial value of loop (i.e xi) , when I
-|select xi=-7. , the answers are completely different with when I
-|select xi=-6 !!!!
-|I can't understand logic of program ! answers are near -1.0 , why xi
-|affect them ??
-|Can anybody say me how fortran thinks about this code ??
-|Best,nakisa

Dr Ivan D. Reid

2007-04-06, 7:07 pm

On Fri, 06 Apr 2007 08:56:54 -0400,
Herman D Knoble <SkipKnobleLESS@SPAMpsu.DOT.edu>
wrote in <qogc13943pl3ss081i9oo280u1o9pi6t1o@4ax.com>:
> Nakisa: Try the following slight modification of your code

> module SetPrecision
> implicit none
> integer, public, parameter :: DP = selected_real_kind(15)
> end module SetPrecision

> program amplification

Should there be an 'implicit none' here as well? I recall getting
bitten once because I'd thought that USEing a module with IMPLICIT NONE
carried over to the USEing routine (guess I was thinking alog the lines
of a C "#include"). Or maybe it was vice-versa...

> use SetPrecision
....

--
Ivan Reid, School of Engineering & Design, _____________ CMS Collaboration,
Brunel University. Ivan.Reid@[brunel.ac.uk|cern.ch] Room 40-1-B12, CERN
KotPT -- "for stupidity above and beyond the call of duty".

nakisa

2007-04-08, 7:03 pm

hi and thanks alot
your hints about using non-integer steps in loop make my mind clear
and
by changing form false "do " ,the error has been removed .

i have another problem here . i wrote k=0.0001 in program ( i define
k as real variable ) ,and then
say fortran to wrie k , it writes : 9.9999997E-06 .
in one part of my program numbers will compare with k , now i there is
diffrence between i think and what fortran
calculaute ! am i wrong ?? why fortran computes k with this form ??
thanks alot
nakisa

Richard Maine

2007-04-08, 7:03 pm

nakisa <nakisa.nooraee@gmail.com> wrote:

> i have another problem here . i wrote k=0.0001 in program ( i define
> k as real variable ) ,and then
> say fortran to wrie k , it writes : 9.9999997E-06 .

This is closely related to your DO loop problem and understanding it is
absolutely critical to doing anything much useful with scientific
programming.

Almost all numbers in computers today are stored as binary - not
decimal. There are exceptions, but they are rare and for your current
purposes, they can be ignored. (For example, there are no exceptions in
current Fortran compilers).

The number 0.0001 cannot be represented exactly in a finite binary
representation. I repeat the "cannot". Nothing you can do will change
that. Anything you code that depends on having this number exactly will
fail. The 0.0001 that you write is converted to the closest
approximation available in the particular binary representation being
used. When you write the number back out, the binary number is converted
to a decimal character representation for the printout. That conversion
will also be an approximation.

Let me repeat that understanding this is *CRITICAL* to scientific
programming. Floating point arithmetic involves approximations. You
cannot afford to sweep the matter aside. A large part of the study of
scientific computation involves analysis of those approximations and
their effects.

--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain

nakisa

2007-04-08, 10:04 pm

Richard , i am compeletly agree with you
that understanding is ciritcal in scientific programming.
in my program , i need values of function be compared with to 0.001 .
is there any way that i
can work with such precesion ?
best ,nakisa

glen herrmannsfeldt

2007-04-09, 7:04 pm

nakisa wrote:
> Richard , i am compeletly agree with you
> that understanding is ciritcal in scientific programming.
> in my program , i need values of function be compared with to 0.001 .
> is there any way that i
> can work with such precesion ?

You should at least use 0.001D0 to get double precision.
It still won't be exact, as no binary fraction will be, but it
will be much closer.

-- glen

John Harper

2007-04-10, 7:05 pm

In article <1hw8yeh.gwdomt1o6smz0N%nospam@see.signature>,
Richard Maine <nospam@see.signature> wrote:
>nakisa <nakisa.nooraee@gmail.com> wrote:
>
>
>This is closely related to your DO loop problem and understanding it is
>absolutely critical to doing anything much useful with scientific
>programming.

I agree with everything Richard said. The "executive summary" is the
quote I gave here in July 2002:

"Floating-point numbers are like piles of sand. Every time you do
anything with a pile of sand you lose a little sand and you pick
up a little dirt."

I asked then who said that. I still don't know.

-- John Harper, School of Mathematics, Statistics and Computer Science,
Victoria University, PO Box 600, Wellington 6140, New Zealand
e-mail john.harper@vuw.ac.nz phone (+64)(4)463 5341 fax (+64)(4)463 5045

Steven G. Kargl

2007-04-10, 7:05 pm

In article <1176244069.752464@bats.mcs.vuw.ac.nz>,
harper@mcs.vuw.ac.nz (John Harper) writes:
>
> "Floating-point numbers are like piles of sand. Every time you do
> anything with a pile of sand you lose a little sand and you pick
> up a little dirt."
>
> I asked then who said that. I still don't know.
>

google claims that it is in Kernighan and Plauger.

--
Steve
http://troutmask.apl.washington.edu/~kargl/

Walter Spector

2007-04-10, 7:05 pm

"Steven G. Kargl" wrote:
>
> In article <1176244069.752464@bats.mcs.vuw.ac.nz>,
> harper@mcs.vuw.ac.nz (John Harper) writes:
>
> google claims that it is in Kernighan and Plauger.

K&P mention the quote in their book "The Elements of Programming
Style". The actual quote (2nd edition, page 117) is:

As a wise programmer once said, "floating point numbers are like
sandpiles: every time you move one, you lose a little sand and you
pick up a little dirt."

However K&P do not mention who the "wise programmer" was...

W.

John Harper

2007-04-10, 10:03 pm

In article <evh3pn$dht$1@gnus01.u.washington.edu>,
Steven G. Kargl <kargl@troutmask.apl.washington.edu> wrote:
>In article <1176244069.752464@bats.mcs.vuw.ac.nz>,
> harper@mcs.vuw.ac.nz (John Harper) writes:
>
>google claims that it is in Kernighan and Plauger.

Thank you - so it is. But they don't claim to have said it first.
The exact quote from K&P p.92 is 'As a wise programmer once said,
"Floating point numbers are like sandpiles: every time you move one,
you lose a little sand and you pick up a little dirt."'

Dennis Ritchie said in comp.compilers on 29 Sep 1998
"Kernighan just assured me that it was indeed Vic Vyssotsky
from whom the quote was taken (my backup guess was Hamming)."

-- John Harper, School of Mathematics, Statistics and Computer Science,
Victoria University, PO Box 600, Wellington 6140, New Zealand
e-mail john.harper@vuw.ac.nz phone (+64)(4)463 5341 fax (+64)(4)463 5045

Steven G. Kargl

2007-04-10, 10:03 pm

In article <1176254785.693477@bats.mcs.vuw.ac.nz>,
harper@mcs.vuw.ac.nz (John Harper) writes:
> In article <evh3pn$dht$1@gnus01.u.washington.edu>,
> Steven G. Kargl <kargl@troutmask.apl.washington.edu> wrote:
>
> Thank you - so it is. But they don't claim to have said it first.
> The exact quote from K&P p.92 is 'As a wise programmer once said,
> "Floating point numbers are like sandpiles: every time you move one,
> you lose a little sand and you pick up a little dirt."'

Yeah, after posting I went down the hall to the APL library and
found the passage. I also scanned Knuth's tome, and it appearred
that he did not write the passage.

> Dennis Ritchie said in comp.compilers on 29 Sep 1998
> "Kernighan just assured me that it was indeed Vic Vyssotsky
> from whom the quote was taken (my backup guess was Hamming)."

Interesting. I would have guessed Cody or Kahan. Perhaps,
an email to Kernighan is in order.

--
Steve
http://troutmask.apl.washington.edu/~kargl/

John Harper

2007-04-22, 7:04 pm

In article <1176254785.693477@bats.mcs.vuw.ac.nz>,
John Harper <harper@mcs.vuw.ac.nz> wrote:

>"Floating point numbers are like sandpiles..." [: Vic Vyssotsky]

G95 offers a good example because it can do the same(?) thing in two
different ways, giving two very slightly different answers. Andy said
that in g95 compile-time math is done by fp emulation, but run-time
math is done by the cpu itself. I'm not sure if this is an appropriate
use of the transfer intrinsic, though.

Both answers are of course as good as you can expect from double
precision. The test program and its output, with either NetBSD or
Sun Solaris, are:

PROGRAM testdpconst ! initialization sometimes /= assignment ?
INTEGER,PARAMETER :: long = selected_int_kind(15)
DOUBLE PRECISION :: dx = 2d0
PRINT 666,'sqrt(dx)**2 ',sqrt(dx)**2 &
, 'in hex format',transfer(sqrt(dx)**2,1_long)
PRINT 666,'sqrt(2d0)**2 ',sqrt(2d0)**2 &
, 'in hex format',transfer(sqrt(2d0)**2,1_long)
666 FORMAT (1X,A," = ",F18.16,/,1X,A," = ",Z16)
END PROGRAM testdpconst

sqrt(dx)**2 = 2.0000000000000004
in hex format = 4000000000000001
sqrt(2d0)**2 = 1.9999999999999996
in hex format = 3FFFFFFFFFFFFFFE

-- John Harper, School of Mathematics, Statistics and Computer Science,
Victoria University, PO Box 600, Wellington 6140, New Zealand
e-mail john.harper@vuw.ac.nz phone (+64)(4)463 5341 fax (+64)(4)463 5045

Sponsored Links