Wrong results from Fortran function called from Julia (0.6.2) on Windows

Question

I found a small helper subroutine from a larger legacy codebase that computes some statistics for a vector and have compiled it alone and am trying to call it from Julia (0.6.2) for practice.

To begin, I have found some forum posts and discourse threads, but I wanted to double check since the language is developing quickly and some may be using outdated syntax like this one. Additionally, I am using Windows and many of the forum posts I came across reference .so files, which to my understanding is for UNIX systems, not Windows?

The Fortran code I am referencing (called meansd.f) in its entirety is:

SUBROUTINE MEANSD (A,N,ABAR,VAR,STD)
IMPLICIT NONE
C     COMPUTE THE MEAN AND STANDARD DEVIATION OF A REAL VECTOR.
C     PARAMETERS:
C     A: REAL VECTOR TO COMPUTE MEAN AND STANDARD DEVIATION OVER.
C     N: LENGTH OF A.
C     ABAR: COMPUTED MEAN OF A.
C     VAR: VARIANCE OF A.
C     STD: COMPUTED STANDARD DEVIATION OF A.
C
REAL A(N)
INTEGER I,N
DOUBLE PRECISION APX, ABAR, VAR, STD, SUM, SUMSQ, DN, Z
C
DN=N
APX=0D0
DO 30 I=1,N
APX=APX+A(I)
30 CONTINUE
APX=APX/DN
SUM=0D0
SUMSQ=0D0
DO 40 I=1,N
Z=A(I)-APX
SUM=SUM+Z
SUMSQ=SUMSQ+Z*Z
40 CONTINUE
ABAR=SUM/DN
VAR=(SUMSQ-(SUM*ABAR))/(DN-1)
STD=DSQRT(VAR)
ABAR=ABAR+APX
RETURN
END

I have compiled this using gcc

gcc (x86_64-posix-seh-rev1, Built by MinGW-W64 project) 7.2.0

largely following this question, which completes with no errors:

gfortran -c meansd.f
gfortran -shared -fPIC  meansd.o -o meansd.dll

when I test with Libdl, everything seems to return fine and I get returned a pointer

path=joinpath(pwd(),"meansd.dll")
dl=Libdl.dlopen(path)

and then I try to build my ccall function as:

a_ref=Ref{Array{Float64,1}}([5.0,5.0,5.0,5.0,5.0])
i_ref=Ref{Int64}(1)
n_ref=Ref{Int64}(length([5.0,5.0,5.0,5.0,5.0]))
apx_ref=Ref{Float64}(0)
abar_ref=Ref{Float64}(0)
var_ref=Ref{Float64}(0)
std_ref=Ref{Float64}(0)
sum_ref=Ref{Float64}(0)
sumsq_ref=Ref{Float64}(0)
dn_ref=Ref{Float64}(0)
z_ref=Ref{Float64}(0)

ccall((:meansd_,"meansd.dll"),Void,(
Ref{Array{Float64,1}},#A
Ref{Int64}, #I
Ref{Int64}, #N
Ref{Float64}, #APX #I think I have to pass these too?
Ref{Float64}, #ABAR
Ref{Float64}, #VAR
Ref{Float64},  #STD
Ref{Float64}, #SUM
Ref{Float64}, #SUMSQ
Ref{Float64}, #DN
Ref{Float64}, #Z
),a_ref,i_ref,n_ref,apx_ref,
abar_ref,var_ref,std_ref,sum_ref,sumsq_ref,dn_ref,z_ref)

This completes which is as far as I have gotten, but abar_ref is returned as NaN. I then made a Float/Int 32 version as well, which completes but abar_ref is returned as 0.0. I thought I should be using 64 bit, since I compiled with 64 bit gfortran.

I think I am messing up on how to pass the array, but I have tried as many different approached I could think of with no success.

I am also confused as to when to use Julia types vs the ccall types (Float32 vs Cint?). I have seen examples of both scattered among the web, but if I am reading correctly, there are underlying convert statements built in behind the scenes at this point which handle that?

Hi, first issue to check: in fortran, the default "REAL" and "INTEGER" types depend on the compiler. On typical 64-bit systems, the common types are 32-bit for integers and reals (but again, this is "processor dependent"). In your case, the likely types are 32-bit real, 32-bit integer, then three 64-bit reals. As the calls rely on the c interface, you may want to look up the "iso c binding" feature of Fortran. — Pierre de Buyl, Jan 30 '18 at 08:30
Thanks, that is interesting to know. I had a few questions about the points you brought up. 1. When calling external libraries is it necessary to declare the intermediate variables (`APX` for example) or just the listed inputs/outputs? 2. How well do the C ISO standards work in terms of F77 code? [some quick searching makes it seem like this is represented more in the later formats](https://stackoverflow.com/questions/2902186/pass-fortran-77-function-to-c-c) — Casey, Jan 31 '18 at 02:20
I wouldn't worry too much about restrictions on F77 code: your example code isn't quite F77. gcc 7.2 has quite good support for C interoperability. (As it's meant to be fixed form though, perhaps you could look into indenting it?) — francescalus, Jan 31 '18 at 06:15
This may be a funny error bu I will admit I didn't know the code was supposed to be indented. I copied and pasted the code from the [source repo](https://sourceforge.net/p/open-fvs/code/HEAD/tree/trunk/base/src/meansd.f) and many of the files throughout the project look like that so I assumed that's how the fortran files were supposed to look.. I am new to compiled languages, I do not really understand the full capabilities of compilers vs the languages themselves. — Casey, Feb 02 '18 at 01:27

Casey · Accepted Answer · 2018-02-02T01:19:45.320

Edit: expanded notes

Following Pierre's comments, I went back and re-read through the documentation and managed to get a working version:

I first removed instantiating any variables other than those named in the function call:

SUBROUTINE MEANSD (A,N,ABAR,VAR,STD).

Next, I re-read through the documentation section on calling FORTRAN, in particular this warning section.

Warning

For string arguments (char*) the Julia type should be Cstring (if NUL- terminated data is expected) or either Ptr{Cchar} or Ptr{UInt8} otherwise (these two pointer types have the same effect), as described above, not String. Similarly, for array arguments (T[] or T*), the Julia type should again be Ptr{T}, not Vector{T}.

The resulting code is:

#declare variables

a=Float32[1.0,2.0,3.0,4.0,5.0]  #changed values from OP to have non-zero var/sd for demo
n=length(a) #works as Int32(length(a)) as well?
abar_ref=Ref{Float64}(0)
var_ref=Ref{Float64}(0)
std_ref=Ref{Float64}(0)

#call function

ccall((:meansd_,"meansd.dll"),Void,(
Ptr{Float32},#A
Ref{Int32}, #N
Ref{Float64}, #ABAR
Ref{Float64}, #VAR
Ref{Float64},  #STD
),a,n,abar_ref,var_ref,std_ref)

#check results
abar_ref[] #returns 3.00

var_ref[] #returns 2.50

std_ref[] #returns 1.58

added comments in answer above referring to pertinent points in the docs that helped me. I also found another [smaller helper function](https://sourceforge.net/p/open-fvs/code/HEAD/tree/trunk/base/src/tvalue.f) in the larger project and have I tried to get that working following what I learned here but have not been able to do so successfully yet, so my grasp on this is probably tenuous at best. — Casey, Feb 02 '18 at 01:31

Wrong results from Fortran function called from Julia (0.6.2) on Windows

1 Answers1