real compare should work but doesn't

Hi - new to silverfrost and also learning FORTRAN. I have a short Fortran 95 program which hits an overflow because 2 variables get the same value at one point. Should be easy to manage that, right? But all my comparisons don't work, see code below for the latest try....

The statement that should pick it up is: if (abs(tanh(SthisX)-tanh(SlastX)).EQ.zero) ...... but whenever tanh(SthisX) becomes the same value as tanh(SlastX), (both =1.000E0) the program crashes with an overflow - AFTER the if condition. IE the If condition doesn't pick it up, but the calculation after that does and divides by 0. I have also tried testing the difference against my 'zero' parameter, same result.

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

It appears that your program is written to solve tanh(x) = 0, which has just one root, x = 0. If you really meant to solve that trivial equation instead of the more likely equation tanh(x) = C, given a value of C not equal to zero, there was no need to write a program. Your program is working reasonably, but since the estimates of the root approach zero, the secant formula will involve calculation of 0 / 0 when x is too small to represent in a machine representation.

In any case, you do not need to write a program to obtain the root, since we know that 2 arctanh(C) = ln [ (1 + C) / (1 - C)].

A few ideas ...

It would be a good idea to try to indent the code to make it easier to see its structure.

Do you really want to test for the absolute difference in the two function values to be equal to your 'zero' parameter? Should you replace .EQ. with .LE. ?

Your constants are default real kind and [u:64e7e6ff7a]not[/u:64e7e6ff7a] real(kind=dp). You need to use the following to get the precision you want:

real(kind=dp), parameter :: zero=TINY(0.0_dp), infinity=HUGE(0.0_dp), root=0.0_dp

Print these out or use the debugger to see the difference between, for example, tiny(0.0) and tiny(0.0_dp).

Generally, using tiny(0.0_dp) for your tolerance value is too small. Just because you can represent a value doesn't mean you can divide by it. So consider using a larger value.

Thanks for replies, all very useful. . My code is indented in silverfrost, but it gets lost in the copy&paste to the forum - how do I keep that formatting please? . It is a trivial problem - designed only to demonstrate the method using a simple function, for the 'single root'. I'm learning Fortran as part of a degree by correspondence, so the extra support is really helpful. . Testing with .LT. instead of .EQ. did the trick , I will take comparing reals much more seriously now 😃 . Thanks for the new knowledge of '_dp', the course I am doing kind of assumes we can teach ourselves fortran so all pure fortran advice is much appreciated.

Hi again, as an aside - I was thinking about dp, does anyone know why we declare it like this? integer, parameter :: dp = selected_real_kind(15, 307) Why use 'integer'? I think because I don't understand it well enough, it looks a little clumsy?

Why is dp declared integer? Because kind numbers are integers that serve to define the size of the type. The particular selection of numbers depends on the compiler. Some compilers use 4 and 8 for the kind number of 4-byte and 8-byte reals. FTN95 uses kind=1 and 2 for these types. The name 'dp' is user selected, as it is suggestive of 'double precision'.

To retain the indents, use the 'code' button before and after the pasted info. You will soon find that there is a limit to the length of a code section and you may have to post it as several individual posts. Ian

Quoted from IanLambley You will soon find that there is a limit to the length of a code section and you may have to post it as several individual posts.

Indeed it is so. I believe that phpBB (the engine behind this forum -- bottom of page) BBCode could be updated so the 'code' window would include a vertical scrollbar for long listings. I think it is just a question of changing the html code that [code] links too. (Note as users we can't use HTML directly on this board)

A belated thanks to all, I'm making progress 😉