Wednesday, 12 March 2008

Another Fortran idiosyncracy bites me on the bum

I thought I'd have this program done by the end of the afternoon. So I start to write a function to calculate a nice simple Gaussian function. No NaNs and +Infs to contend with here, I thought. How wrong I was.
This time, I was finding that the intrinsic function exp(x) returned zero for x less than around -0.1E3. After some bashing, my compiler was sweet enough to tell me:
Result of EXP underflows its kind
Apparently, Fortran is not content with just "exp". It also has "dexp" and "cexp" for double-precision and complex floating point numbers respectively. So, dexp(x) it is then. But be careful how you define x, or you will get a shouty like this:
Type of argument 'x' in call to 'dexp' at (1) should be REAL(8), not REAL(4)
So one should try dexp(-0.1D3) (note the letter D!) rather than exp(-0.1E3), which will return zero, possibly grumbling in the process.

Monday, 10 March 2008

NaN, +Inf, -Inf, horrid things...

This is the answer.
Based on such, I have a sample function to detect a +Inf:

function ispinf(x)
! function to detect if x is +Inf
integer :: ispinf
double precision :: x

!local variables
integer :: IPInf
real :: PInf
data IPInf/B'01111111100000000000000000000000'/ ! +Infinity
PInf = transfer(IPinf,Pinf)

if (x.eq.PInf) then
ispinf = 1
ispinf = 0

end function ispinf