[ognik]

I have a Fortran 95 program (using silverfrost compiler) which hits an overflow because 2 variables are the same. Should be easy to manage that, right? But all my comparisons don't work, see code below.

The statement that should pick it up is:
if (abs(tanh(SthisX)-tanh(SlastX)).EQ.zero) ......
but whenever tanh(SthisX) becomes the same as tanh(SlastX), the program crashes with an overflow - AFTER the if condition. IE the If doesn't pick it up.

Any ideas? Thanks

--------------------------------------------------
program secTanHroots
!Summary: find positive root of tanh(C) using Secant method
implicit none
integer, parameter :: dp = selected_real_kind(15, 307)
real(kind=dp) :: guess, SthisX, SnextX, SlastX, Sdiff
real(kind=dp), parameter :: zero=TINY(0.0), infinity=HUGE(0.0), root=0.0
real :: C=5
integer :: iter
integer, parameter :: n=58
character(len=n ) :: stars=REPEAT('*',n)

!----- header ----
print *, " "
print *, 'Program secTanHroots starts'
print 5, 'Enter a guess to estimate the positive root of Tanh(', C,')'
print *, ' (enter a negative number to stop)'
read (*,*) guess
print *, stars
print *, " "
print 10, 'Iter.','X(i) ', 'X(i+1)"root"', 'tanh(SthisX)','tanh(SlastX)', 'diff'
!----- end header ----

do while (guess.GE.0) ! allows trials of diff. values of guess
!----- Initialise ----
SthisX= guess
Sdiff=1
SlastX=1
iter=0
Do while (abs(Sdiff).GT.zero)
iter=iter+1
print *, tanh(SthisX)-tanh(SlastX)
if (abs(tanh(SthisX)-tanh(SlastX)).EQ.zero) then
print *, '++++++++++++++++ overflow +++++++++++++++'
go to 100
end if
SnextX= SthisX-tanh(SthisX)*((SthisX-SlastX)/(tanh(SthisX)-tanh(SlastX))) !Secant method
Sdiff = root-SnextX
print 20, iter, SthisX, SnextX, tanh(SthisX),tanh(SlastX), Sdiff

SlastX=SthisX
SthisX=SnextX
if (iter.GT.20) then ! to prevent overrun
print *, 'forced stop, cntr > 20'
stop
end if
end do

100 print *, stars
print *, " "
print *, 'Enter another guess (or negative number to stop)'
read *, guess
print *, " "
end do

print *, 'User selected end'
print *, stars
5 format (A, F3.1,A)
10 format (A, 4X,5(A, 6X))
20 format (I2, 5(3X, ES14.6))
end program secTanHroots

[opedrofunk]

If "SthisX == SlastX" then you are dividing by zero.

SnextX= SthisX-tanh(SthisX)*((SthisX-SlastX)/(tanh(SthisX)-tanh(SlastX))) !Secant method

Also, the correct way to compare floating point numbers is not with .eq., but with an error tolerance. For example:

if (abs(sthisx - slastx) <= epsilon(real) then ...

Hope this helps,
Peter

[ognik]

The .LE. instead of .EQ. sorted it thanks. I guess I was in denial about comparing reals in Fortran :-)