So here is the fortran90 with C problem

This is the example I concocted to help demonstrate the problem and help me learn how to solve it.

The problem is, I solved this … without any fancy interface blocks or anything strange like that.

Go figure.


The makefile


CC	= icc
FC	= ifort
LD	= ifort
CFLAGS	= -g
FFLAGS	= -g
LFLAGS	= -g


all:	test.exe

test.exe:	test_f.o test_c.o
	$(LD) $(LFLAGS) test_f.o test_c.o -o test.exe

test_f.o:	test.f90
	$(FC) $(FFLAGS) -c $< -o test_f.o

test_c.o:	test.c
	$(CC) $(FFLAGS) -c $< -o test_c.o


clean:
	rm -f *.o *.exe

the test.f90 code:


	program test
	
	integer, allocatable, dimension(:,:)	:: a
	integer			  		:: b
	integer			  		:: n,m,i,j
	real, allocatable, dimension(:,:)	:: c
	
	! -- allocate a and c as 10x10 arrays, and fill with sequential integers
	allocate(a(10,10))
	allocate(c(10,10))
	
	do i=1,10
	 do j=1,10
	  a(i,j)=10*(i-1)+j
	  c(i,j)=dble(10*(i-1)+j)
	 enddo
	enddo
	
	n=10
	m=10
	
	print *,'In test.f90: printing out 1 column of a and then c'
	
	print *,'a:'
	do j=1,10
	 print *,'a(1,',j,') = ',a(1,j)
	enddo
	
	print *, ' '

	print *,'c:'
	do j=1,10
	 print *,'c(1,',j,') = ',c(1,j)
	enddo 
	call test_c(n,m,a,c)
	
	end program

the test.c code:


#include 
void test_c_(
	      int *n,
	      int *m,
	      int *a,
	      float *c
	    )
{
  int i,j;
  printf("in test.c: n=%i, m=%i, a=%x, c=%x\n",*n,*m,a,c);
  for(i=0;i<*m;i++)
   {
    printf("a(%i,0)=%i\n",i,*(a+ (*m)*i+0));
   }
  for(i=0;i<*m;i++)
   {
    printf("c(%i,0)=%f\n",i,*(c+ (*m)*i+0));
   }

  return;
}

What is of interest is that Fortran includes metadata in its internal symbol tables that C doesn't have access to. So you can't write the C arrays as a[i][j], as this symbol metadata never traversed the language barrier. You have to do pointer arithmetic. Notice how it is done.


     a[i][j] -> *(a + j * DIMENSION_ALONG_J + i)

which is equivalent to a(j-1,i-1) in Fortran90.

Nice little Rosetta stone we have here 🙂

When you run the code, this is what you get

~/testcase$ ./test.exe 
 In test.f90: printing out 1 column of a and then c
 a:
 a(1,           1 ) =            1
 a(1,           2 ) =            2
 a(1,           3 ) =            3
 a(1,           4 ) =            4
 a(1,           5 ) =            5
 a(1,           6 ) =            6
 a(1,           7 ) =            7
 a(1,           8 ) =            8
 a(1,           9 ) =            9
 a(1,          10 ) =           10
  
 c:
 c(1,           1 ) =    1.000000    
 c(1,           2 ) =    2.000000    
 c(1,           3 ) =    3.000000    
 c(1,           4 ) =    4.000000    
 c(1,           5 ) =    5.000000    
 c(1,           6 ) =    6.000000    
 c(1,           7 ) =    7.000000    
 c(1,           8 ) =    8.000000    
 c(1,           9 ) =    9.000000    
 c(1,          10 ) =    10.00000    
in test.c: n=10, m=10, a=129a030, c=129a1d0
a(0,0)=1
a(1,0)=2
a(2,0)=3
a(3,0)=4
a(4,0)=5
a(5,0)=6
a(6,0)=7
a(7,0)=8
a(8,0)=9
a(9,0)=10
c(0,0)=1.000000
c(1,0)=2.000000
c(2,0)=3.000000
c(3,0)=4.000000
c(4,0)=5.000000
c(5,0)=6.000000
c(6,0)=7.000000
c(7,0)=8.000000
c(8,0)=9.000000
c(9,0)=10.000000

So it works. Without an interface block. And indexing is basically correct. The hex address pointers appear to be valid.

Now how to translate this to the code we are having trouble with ...

[updated ... fixed my indexing errors...]

Viewed 7997 times by 1662 viewers

4 thoughts on “So here is the fortran90 with C problem

  1. You didn’t already know that? Seriously? You never checked any code that crosses Fortran and C/C++, like GNU Octave, Cactuscode, etc.? Or any of the BLAS implementations in C that are called from LAPACK (and hence can be called with an ALLOCATABLE array from one level up)?

    A[i][j] in C means find the ith pointer in the array addressed by A and then the jth element in the array of that pointer. It’s not a multidimensional array. Nothing to do with symbol tables. If you declare int A[2][3];, the compiler acts differently, yes. That’s an inherited legacy that does not play well with the rest of the language.

    If you follow a convention of X and ldX for naming the array and the leading dimension, try #define AREF(A,i,j) (A)[i + j * LDNAME(A)] and #define LDNAME(A) ld##A. You need the second expansion (LDNAME) to force name expansion if you pass a macro as AREF’s A argument. Then just AREF(A, i, j) throughout, and &AREF(A, i, j) works as expected.

  2. @anon

    No I knew about the differences between C and Fortran indexing arrays … its a simple polynomial evaluation, and I just finished coding it into some inline-able functions to make the debugging easier. I also know quite well that a[2][3] is not the same as *(a + 2*lda +3) in a general sense, if a has been built as a pointer array to pointers.

    Fortran doesn’t do this. There they try to allocate one big contiguous memory region for the array. Tends to make optimizations easier at the compiler level.

    Been using LINPACK/EISPACK since the late 80s, and LAPACK more “recently” (early 90s). Did this as functions rather than macros for the moment. Easier to debug the functions …

    Everything is “working” though we are now getting a huge array that looks like it was allocated on the heap rather than the stack, and hitting a segv as soon as we try to index it … which is very strange, as I print values from this array previously in this same routine.

    Likely a pilot/editor error.

  3. Well, it kind of makes sense. Some time ago I was trying hard to call lapack functions from C and then I realized these nifty little things.

    BTW: is the dynamic memory allocation in fortran supposed to happen in heap or stack? Heap, I’d have imagined.

  4. @Rohit

    It looks like it uses a simple rule. If the allocation is below some (fixed) size, it happens on the stack, otherwise it happens in heap. gcc/gfortran and the intel compilers let you play with this a little.

    What I have found more recently is a need to use the mcmodel=medium options in compilers rather than the default (small). This is a data size issue, if your data is larger than 4GB of ram, you need this.

Comments are closed.