FORTRAN calls C

[jayfortranuser]

Hi all,

The error of array bounds exceeded occurred when I use Compaq Visual Fortran 6.6 and Microsoft Visual C++ 6.0 to compile a very complex model developed by Microsoft Fortran PowerStation v4.0 and Visual C++ 4.0. Turning array bounds checking off is not an option because the results will be totally wrong.

I found out that after a C subroutine is called by FORTRAN, some variables not involved with the C subroutine is changed (for example, the layer was changed from 1 to 1148846080 and the time was changed from 1997 to 0). So, I think the problem must be related with the interaction of the two languages. The following is an example of FORTRAN calls C (csched.f calls flow.c). I wish you could find the error. Thank you. Jay

FORTRAN subroutine:

subroutine csched(cflow,protop,probot,

&froma1,tob1,

&froma2,tob2,

&frac,accum)

implicit none

include 'const.inc'

include 'zztim.inc'

real cflow, protop, probot, froma1, tob1, froma2, tob2,

&frac,accum(2)

INTERFACE

SUBROUTINE flow(from, to, when, howmuch)

!MS$ATTRIBUTES ALIAS:'_flow' :: flow

REAL from

REAL to

REAL when

REAL howmuch

END SUBROUTINE flow

END INTERFACE

real cflow1, cflow2

cflow2 = cfl ow * (protop / probot) * frac

cflow1 = cflow - cflow2

call flow(froma1,tob1,time,cflow1)

call flow(froma2,tob2,time,cflow2)

accum(UNLABL) = accum(UNLABL) + cflow1

accum(LABELD) = accum(LABELD) + cflow2

return

end

C subroutine:

#include

#include

#include "flow.h"

void flow(float *from, float *to, float *when, float *howmuch)

{

/* Increment the number of flows stored in the stack */

nflows += 1;

if (nflows > LENFST) {

/* Stack Overflow */

flow_err(1, *when);

exit(1);

} else {

/* Store the arguments in the stack */

flowstack[nflows].from = from;

flowstack[nflows].to = to;

flowstack[nflows].when = *when;

flowstack[nflows].amt = *howmuch;

}

return;

}

[Jugoslav_Dujic]

INTERFACE

SUBROUTINE flow(from, to, when, howmuch)

!MS$ATTRIBUTES ALIAS:'_flow' :: flow

Uh-oh. Don't do this. You're tweaking the interface to avoid linker error, but you're basically just lying to the compiler. The default calling convention for both FPS4 and CVF is __stdcall, while the default calling convention for C is __cdecl. They manage stack differently (and, also, IIRC, the order of arguments is reverse) so things can get really nasty in the run-time.

Instead, you have to tell the compiler to use C calling convention for flow:

INTERFACE
SUBROUTINE flow(from, to, when, howmuch)
!MS$ATTRIBUTES C, REFERENCE, ALIAS:'_flow' :: flow

I didn't look carefully at the rest of the code, but this should make things behaving much better. (Actually, ALIAS: '_flow' is now spurious, as compiler will expect '_flow' juston thebasis of C attribute).

Jugoslav

[jayfortranuser]

Thank you Jugoslav,

Should I add "C, REFERENCE" to all ITERFACE subroutines?

Should I change !MS to !DEC$ ?

Jay

[Jugoslav_Dujic]

C attribute sets _cdecl calling convention and implies that arguments are passed by value. REFERENCE attribute for the routine itself implies that arguments are passed by reference. Thus, if all your INTERFACEs are for C routines and all their arguments are pointers (references), yes, you need them everywhere. If you have mixed-case of value/reference somewhere, you should specify VALUE/REFERENCE attribute for individual arguments.

!MS$ vs. !DEC$ doesn't really matter -- the compiler recognizes both.

(Luckily, we didn't get !CPQ$, !HP$ and !INT$ in the meantime :-) ).

Jugoslav

[Steven_L_Intel1]

Don't think that we didn't get management suggesting that we change the prefix! Note that in Intel Visual Fortran, the meaning of DEC in the directive prefix is "Directive Enhanced Compilation". It's just a coincidence that it used to mean something else. :-)

[jayfortranuser]

Thank you Jugoslav.

Most of the routines in this model are FORTRAN calls C, but there is one occation that C calls FORTRAN. Can you have a look to seeif it is approaperate for the new FORTRAN/C versions? Thanks. Hong

C:

#include
#include
#include
extern float tcalc(float stemp, float teff[4]);

void calcdefac(int *texture, float *stemp, float *tfunc, float *bgwfunc,
float *agdefac, float *bgdefac, float *avgwfps,
float *teff, float *rprpet, int *idef, float *ppt)
{
.......

*tfunc = tcalc(*stemp, teff);

......
return;
}

FORTRAN:

real function tcalc(stemp, teff)

implicit none

real stemp, teff(4)

!MS$ATTRIBUTES C :: tcalc

real catanf
external catanf

real normalizer

normalizer = catanf(30.0, teff(1), teff(2), teff(3), teff(4))

tcalc = max(0.01,
& catanf(stemp, teff(1), teff(2), teff(3), teff(4)) /
& normalizer)

return
end

[Jugoslav_Dujic]

It looks fine -- ATTRIBUTES Cmeans _cdecl & call by value, and the C prototype extern float tcalc(float stemp, float teff[4]) matches that (arrays like teff are always passed by reference no matter what attributes you specify, so that's fine as well).

Only, take care that you cannot call tcalc from the rest of the Fortran code if you need to (unless you specify the interface block for it).

Jugoslav

[jayfortranuser]

Thank you Jugoslav.

The model ran smoothly now. However, results areslightly different like an annual accumulation is 0.105292 compared with 0.105231, which is ignorable for my research. I am just wondering why it occurred. Is that possible that different software versions (FORTRAN/C++) or different project setting could cause slight difference in results? Both default real and integer kinds are 4.

If not, maybe something still wrong with the compiling.

Thank you again.

[Jugoslav_Dujic]

Could be pretty much anything... what did you compare results with?

Theory of floating-point arithmetics is not really my field. If I'd have to bet, I'd say that it could be because x86 processors can (under certain levels of optimization) perform intermediate arithmetics in 80-bit FP registers, and only after doing full computation of the expression cast the result down to 64-bit/32-bit. There's an option "Enable Floating-point consistency" on "Fortran/Floating point" project settings, which, um, affects that, so you could play with it.

I'm sure Steve or someone else could tell you more about the subject.

Jugoslav

[Steven_L_Intel1]

There are many possible reasons for differences in floating point results, including different rounding, different order of operations, and more. There are numerous books and articles on the general topic of floating point - you could start with the one our Dave Eklund wrote called "The Perils of Real Numbers", published in the Visual Fortran Newsletter issues 8-10.

[jayfortranuser]

Thank you Jugoslav and Steve.

I am using Compaq Visual Fortran 6.6 and Microsoft Visual C++ 6.0 to compile a very complex model developed by Microsoft Fortran PowerStation v4.0 and Visual C++ 4.0.

The only thing I did was to add REFERENCE C in the MS$ATTRIBUTES statements, according to Jugoslavs advice. Differences in results (many outputs) between calculated using old executable and calculatedusing new executable were very small (about 0.05%), which were ignorable for my research. I worried that if I have different inputs, results might be wrong. If I cant solve the problem I would run both versions each time, but soon I will change the model. I have tried all different floating point settings, but got the same results. Are there any way we can find the problem, if there is any, without checking every subroutine?

Thanks

Jay

[Steven_L_Intel1]

Analyzing such differences is very difficult and time-consuming. You don't know which result is closer to "correct", and you also have to consider the precision and accuracy (two different concepts) of the input data, plus what is done to the data. Numerical analysis is a complex process, and if the difference really doesn't matter to you, it may not be worth going down this road.

[jayfortranuser]

Thank you.

Jay