Updating variable in a Fortran do while loop

Question

I'm trying to fix a problem with a code written in Fortran 90 thant seeks the solution of a non-linear equation by means of the secant method. As usual, it iterates until the difference between two subsequent iterations (called 'xdelta') is less than a fixed tolerance ('tolerr=1.0e-2'). But the DO WHILE loop does not seem to update the value of delta, and anytime the loop restarts delta takes its initial value (that I conventionally fixed to 2*tolerr):

real(kind=8) function secant(xa,xb,func) !This program find the solution of a generic equation f(x)=0 written in the external function "f_x":

real(kind=8) function secant(xa,xb,func)
implicit none
integer :: k,iter
real(kind=8) :: zerok,xa,xb,xdelta,tolerr,xkp1,xk
real(kind=8), external :: func
tolerr=1.0d-2
xdelta=2.0d0*tolerr
xk=xa
do while (xdelta>tolerr)
   write(*,*) '(1)',xdelta,abs(xkp1-xk),tolerr
   xkp1=zerok(xk,xb,func)
   write(*,*) '(2)',xdelta,abs(xkp1-xk),tolerr
   xdelta=abs(xkp1-xk)
   write(*,*) '(3)',xdelta,abs(xkp1-xk),tolerr
   xk=xkp1
   write(*,*) '(4)',xdelta,abs(xkp1-xk),tolerr
end do
secant=xk
return
end function secant

I obtain the following output:

(1) 2.0000000000000000E-002 749999.99998890515 1.0000000000000000E-002
(2) 2.0000000000000000E-002 1.5429133782163262E-003 1.0000000000000000E-002
(3) 1.5429133782163262E-003 1.5429133782163262E-003 1.0000000000000000E-002
(4) 1.5429133782163262E-003 0.0000000000000000 1.0000000000000000E-002
(1) 2.0000000000000000E-002 749999.99998890515 1.0000000000000000E-002
(2) 2.0000000000000000E-002 3.2309407833963633E-003 1.0000000000000000E-002
(3) 3.2309407833963633E-003 3.2309407833963633E-003 1.0000000000000000E-002
(4) 3.2309407833963633E-003 0.0000000000000000 1.0000000000000000E-002
(1) 2.0000000000000000E-002 749999.99998890515 1.0000000000000000E-002
(2) 2.0000000000000000E-002 5.0784338964149356E-003 1.0000000000000000E-002
(3) 5.0784338964149356E-003 5.0784338964149356E-003 1.0000000000000000E-002
(4) 5.0784338964149356E-003 0.0000000000000000 1.0000000000000000E-002
(1) 2.0000000000000000E-002 749999.99998890515 1.0000000000000000E-002
(2) 2.0000000000000000E-002 7.1007488295435905E-003 1.0000000000000000E-002
(3) 7.1007488295435905E-003 7.1007488295435905E-003 1.0000000000000000E-002
(4) 7.1007488295435905E-003 0.0000000000000000 1.0000000000000000E-002
(1) 2.0000000000000000E-002 749999.99998890515 1.0000000000000000E-002
(2) 2.0000000000000000E-002 9.3147760489955544E-003 1.0000000000000000E-002
(3) 9.3147760489955544E-003 9.3147760489955544E-003 1.0000000000000000E-002
(4) 9.3147760489955544E-003 0.0000000000000000 1.0000000000000000E-002
(1) 2.0000000000000000E-002 749999.99998890515 1.0000000000000000E-002
(2) 2.0000000000000000E-002 1.1739104287698865E-002 1.0000000000000000E-002
(3) 1.1739104287698865E-002 1.1739104287698865E-002 1.0000000000000000E-002
(4) 1.1739104287698865E-002 0.0000000000000000 1.0000000000000000E-002
(1) 1.1739104287698865E-002 0.0000000000000000 1.0000000000000000E-002
(2) 1.1739104287698865E-002 2813405822.4623795 1.0000000000000000E-002
(3) 2813405822.4623795 2813405822.4623795 1.0000000000000000E-002
(4) 2813405822.4623795 0.0000000000000000 1.0000000000000000E-002
(1) 2813405822.4623795 0.0000000000000000 1.0000000000000000E-002

then the program stops since it gives a segmentation fault, but for other reasons I know. As you can see, for each (1) and (2) the value of delta is restored to 2*tolerr=2e-2.

I think that the problem must be in the rest of the program, since the function 'secant' regularly works with other functions. The name 'xdelta' is not used in other parts of the program. Could you please have an idea of what may cause the problem?

Welcome, be sure to take the tour. Use tag fortran for all Fortran questions and add a version tag to version specific questions (this one is not version specific). See also please How to Ask and minimal reproducible example, because we need to see your code probably in a more complete example, not just one function. Do you use modules? Be aware that real(8) is ugly and not portable. — Vladimir F Героям слава
– Vladimir F Героям слава, Commented Sep 28, 2017 at 15:22

Pierre de Buyl · Accepted Answer · 2017-09-29 08:06:27Z

1

As you can see, for each (1) and (2) the value of delta is restored to 2*tolerr=2e-2.

For the first five calls to your function, xdelta is effectively lower than tolerr. So, the while loop exits and what you see with the next "(1)" is the next call to the function, hence the reset of xdelta to 0.02.

When the criterion is not met at the sixth iteration, the value of xdelta is retained at the loop "looping" as you would expect. At that point you have other issues that are related to the algorithm. See either MathWorld or WikiPedia for an explanation of the algorithm, for instance.

Also, the line with zerok(xk,xb,func) suggests that zerok is a function but it is only defined implicitly. If you want to ease future debugging, Fortran has more explicit ways to define interfaces. Look for "Fortran modules" and "interface".

EDIT: change comment on definition of zerok

edited Sep 29, 2017 at 8:06

answered Sep 29, 2017 at 7:26

Pierre de Buyl

7,3232 gold badges18 silver badges23 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

Pierre de Buyl Over a year ago

I meant implicitly. I always use explicit interfaces so it trips me up to see implicit ones. I edited the answer accordingly.

agentp Over a year ago

good catch on immediate loop exit. lesson to OP, a couple more debug write statements (eg before the loop) would have revealed that easily.

Collectives™ on Stack Overflow

Updating variable in a Fortran do while loop

1 Answer 1

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related