matlabfortranargumentsnag-fortran

Function returns different answers with same arguments

I'm transitioning from MATLAB to Fortran and encountering all sorts of weird behaviors I'd never expect from MATLAB. Here's one that's got me puzzled:

Program pruebanormal

double precision :: g01eaf, x, y

character :: T*1

integer :: Iffail

Iffail = 0

T = 'L'

x = 0.0

y = g01eaf('L',0.0,Iffail)

write(*,*) 'y = ',y

end program pruebanormal

I have this fairly simple program in which I'm trying to find the pdf at x=0 of a standard N(0,1) variable (should be 0.5). g01eaf() is the NAG library function that does this for me. I'm using gfortran to compile.

Leaving the rest of the program unchanged, depending on how I write the arguments in g01eaf(), I get different answers:

a) g01eaf(T,x,Iffail)

b) g01eaf(T,0.0,Iffail)

c) g01eaf(T,x,0)

Now, under MATLAB, I would get the same (correct) answer either way: y = 0.500000. Under Fortran, however, I get:

a) y = 0.500000

b) y = 1.000000

c) Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0xB766C163
#1  0xB766C800
#2  0xB77763FF
#3  0x804982E in g01eafn_
Violación de segmento (`core' generado)

I have no words for the answer in (b) and no clue what (c) even means.

Solution

  • The very quick answer for the "wrong result" is that in

    y = g01eaf('L',0.0,Iffail)

    you are passing a different type of real variable than in

    double precision x
x = 0.0   ! x is still double precision.
y = g01eaf('L',x,Iffail)

    The function g01eaf probably expects double precision: you should very carefully read NAG's documentation.

    y = g01eaf('L', 0d0, Iffail)

    Now to elaborate.

    You don't want these problems to come about too often. You want to ensure an interface is available for the function call to g01eaf. Your compiler would then complain about passing a real of default kind to the function.

    Assuming you have an up-to-date version of the library, you want to do something like

    use nag_library, only : g01eaf, nag_wp
implicit none

integer Iffail
real(kind=nag_wp) y
y = g01eaf('L', 0._nag_wp, Iffail)

end

    Again, see the documentation. Both for the library, and for meaning of modules, etc.

    For older versions, one should still have a module available, but it may be called something different, and nag_wp may not be defined (meaning you must carefully choose kinds).

    The module will also lead to a complaint that Iffail requires to be able to be set, so must be a variable, not 0. This explains (c).