Call Array valued Fortran function from C 
Author Message
 Call Array valued Fortran function from C

Dear all,

I have a problem that can be reduced to the following situation:

- mainf.f90 (see below)
       a fortran function that returns an array containing
       n-copy's of a given double precision number.
- mainc.c (see below)
       a C-file that would like to use this function.

The problem is that it doesn't work (Bus Error on assigment of
array(I) in the Fortran code). One way or another, the output argument
is not available in the Fortran code.

I know one solution would be to change the function into a subroutine,
adding the result as an extra argument, but this would spoil the whole
design of the project. So I'd like to know if anyone has an
alternative way to solve this.

Best Regards

ps: for the sake of completeness: I'm using the latest stable
gFortran, gcc 4.3 and I'm working on a Intel Mac.

[1]: http://www.*-*-*.com/

---8<---- MAINF.F90 ---8<-------8<-------8<----
function repeatdouble( N, d ) result( array )
    integer, intent( in ) :: N
    double precision, dimension( N ) :: array
    double precision, intent( in ) :: d

    print *, "repeat ", N, " times ..."
    print *, "double ", d, "."
    print *, "array :", array
    do I = 1,N
        print *, "making copy nb", i, " of ", N
        array(I) = d
    end do
end function repeatdouble
--->8------->8------->8------->8------->8------->8----

---8<---- MAINC.C ---8<-------8<-------8<----
#include <stdio.h>
extern void repeatdouble_( double*, int*, double* );
int main(){

        double buf[9] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 };
        double* ptr2buf = buf;
        double d = 3.141592654;
        int  n = 4;
        int i;

        // make n copies of ch in sbf
        repeatdouble_( ptr2buf, &n, &d );

        for( i=0; i<9; i++){
                printf("buf[%i]=%d\n", i, buf[i]);
        }

Quote:
}

--->8------->8------->8------->8------->8------->8----


Thu, 07 Oct 2010 01:33:39 GMT  
 Call Array valued Fortran function from C

Quote:

> I have a problem that can be reduced to the following situation:
> - mainf.f90 (see below)
>        a Fortran function that returns an array containing
>        n-copy's of a given double precision number.
> - mainc.c (see below)
>        a C-file that would like to use this function.
> The problem is that it doesn't work (Bus Error on assigment of
> array(I) in the Fortran code). One way or another, the output argument
> is not available in the Fortran code.

Since C doesn't have array valued functions, this isn't going
to be easy.  You can have functions returning a struct containing
an array, though the array must be fixed size.

Quote:
> I know one solution would be to change the function into a subroutine,
> adding the result as an extra argument, but this would spoil the whole
> design of the project. So I'd like to know if anyone has an
> alternative way to solve this.

(snip)

Quote:
> ---8<---- MAINF.F90 ---8<-------8<-------8<----
> function repeatdouble( N, d ) result( array )
>     integer, intent( in ) :: N
>     double precision, dimension( N ) :: array
>     double precision, intent( in ) :: d

>     print *, "repeat ", N, " times ..."
>     print *, "double ", d, "."
>     print *, "array :", array
>     do I = 1,N
>         print *, "making copy nb", i, " of ", N
>         array(I) = d
>     end do
> end function repeatdouble
> ---8<---- MAINC.C ---8<-------8<-------8<----
> #include <stdio.h>
> extern void repeatdouble_( double*, int*, double* );

This is for a three argument subroutine, not a two argument
function.  It is possible that the calling sequence is the same,
but that is not standard.

Quote:
> int main(){
>    double buf[9] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 };
>    double* ptr2buf = buf;
>    double d = 3.141592654;
>    int  n = 4;
>    int i;
>    // make n copies of ch in sbf
>    repeatdouble_( ptr2buf, &n, &d );

This agrees with a call to a Fortran subroutine, not a function.

- Show quoted text -

Quote:
>    for( i=0; i<9; i++){
>            printf("buf[%i]=%d\n", i, buf[i]);
>    }
> }



Thu, 07 Oct 2010 04:11:20 GMT  
 Call Array valued Fortran function from C

Quote:
> I have a problem that can be reduced to the following situation:
> - mainf.f90 (see below)
>       a Fortran function that returns an array containing
>       n-copy's of a given double precision number.
> - mainc.c (see below)
>       a C-file that would like to use this function.
> The problem is that it doesn't work (Bus Error on assigment of
> array(I) in the Fortran code). One way or another, the output argument
> is not available in the Fortran code.
> I know one solution would be to change the function into a subroutine,
> adding the result as an extra argument, but this would spoil the whole
> design of the project. So I'd like to know if anyone has an
> alternative way to solve this.
> ps: for the sake of completeness: I'm using the latest stable
> gFortran, gcc 4.3 and I'm working on a Intel Mac.
> [1]: http://docs.sun.com/source/806-3593/11_cfort.html
> ---8<---- MAINF.F90 ---8<-------8<-------8<----
> function repeatdouble( N, d ) result( array )
>    integer, intent( in ) :: N
>    double precision, dimension( N ) :: array
>    double precision, intent( in ) :: d
>    print *, "repeat ", N, " times ..."
>    print *, "double ", d, "."
>    print *, "array :", array
>    do I = 1,N
>        print *, "making copy nb", i, " of ", N
>        array(I) = d
>    end do
> end function repeatdouble
> --->8------->8------->8------->8------->8------->8----
> ---8<---- MAINC.C ---8<-------8<-------8<----
> #include <stdio.h>
> extern void repeatdouble_( double*, int*, double* );
> int main(){
> double buf[9] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 };
> double* ptr2buf = buf;
> double d = 3.141592654;
> int  n = 4;
> int i;
> // make n copies of ch in sbf
> repeatdouble_( ptr2buf, &n, &d );
> for( i=0; i<9; i++){
> printf("buf[%i]=%d\n", i, buf[i]);
> }
> }
> --->8------->8------->8------->8------->8------->8----

The first argument to repeatdouble should not be the address
of an array, rather that of the array descriptor.
Disassembling a gfortran call to a similar function, we get
something like:

type, bind(C) :: descr
   integer(C_INTPTR_T) address
   integer(C_INTPTR_T) unknown1
   integer(C_INTPTR_T) unknown2
   integer(C_INTPTR_T) stride
   integer(C_INTPTR_T) unknown3
   integer(C_INTPTR_T) last_element
end type descr
descr DD
double precision, target :: buf(9)
double precision d
integer n

buf = [(n,n=1,9)]
d = 3.141592654
n = 4

DD = descr(C_LOC(buf), 0, 537, 1, 0, n-1)
call repeatdouble_(descr, n, d)

But you should ask the gfortran folks if they can round up some
documentation about gfortran array descriptors for you.  They
are different in every compiler.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end



Thu, 07 Oct 2010 05:54:39 GMT  
 Call Array valued Fortran function from C

Quote:
> But you should ask the gfortran folks if they can round up some
> documentation about gfortran array descriptors for you.  They are
> different in every compiler.

Yep, that's on the TODO list. Here's the short version (in C syntax).

The descriptor for a given type (named here TYPE) is then:

struct {
  TYPE *data;
  size_t offset;
  ssize_t dtype;
  descriptor_dimension dim[7];

Quote:
}

where descriptor_dimension is defined as

typedef struct descriptor_dimension
{
  ssize_t stride;
  ssize_t lbound;
  ssize_t ubound;

Quote:
}

descriptor_dimension;

All fields have the meaning implied by their name, and dtype is such that:

  * (dtype & 0x07) is the rank
  * (dtype >> 6) is the size of an array element

The data component points to the first element in the array.  The offset
field is the position of the element (0, 0, ...). An element is accessed
by data[offset + index0*stride0 + index1*stride1 + ...]

All this information and more can be found in the big comment in the
source at gcc/fortran/trans-types.c, line 949 in the current sources.

--
FX



Thu, 07 Oct 2010 06:23:33 GMT  
 Call Array valued Fortran function from C

Quote:
>> But you should ask the gfortran folks if they can round up some
>> documentation about gfortran array descriptors for you.  They are
>> different in every compiler.
> Yep, that's on the TODO list. Here's the short version (in C syntax).

Thanks for the documentation, FX!  Since I can never seem to get all
my ducks in a row on the first attempt, I tried a pure Fortran test
of my initial recommendation and found some syntax errors.  After
fixing them up I got:

C:\gfortran\clf\repeatdouble>type repeatdouble.f90
function repeatdouble( N, d ) result( array )
use iso_c_binding
   implicit none
   integer, intent( in ) :: N
   double precision, dimension( N ) :: array
   double precision, intent( in ) :: d
   integer I
   print *, "repeat ", N, " times ..."
   print *, "double ", d, "."
   print *, "array :", array
   do I = 1,N
      print *, "making copy nb", i, " of ", N
      array(I) = d
   end do
end function repeatdouble

C:\gfortran\clf\repeatdouble>type test3.f90
program main
   use ISO_C_BINDING
   implicit none
   type, bind(C) :: descr
      type(C_PTR) address
      integer(C_INTPTR_T) unknown1
      integer(C_INTPTR_T) unknown2
      integer(C_INTPTR_T) stride
      integer(C_INTPTR_T) unknown3
      integer(C_INTPTR_T) last_element
   end type descr
   interface
      subroutine repeatdouble(DD,n,d) bind(C,name='repeatdouble_')
         use ISO_C_BINDING
         import descr
         implicit none
         type(descr) DD
         integer(C_INT) n
         real(C_DOUBLE) d
      end subroutine repeatdouble
   end interface
   type(descr) DD
   double precision, target :: buf(9)
   double precision d
   integer n

   buf = [(n,n=1,9)]
   d = 3.141592654
   n = 4

   DD = descr(C_LOC(buf), 0, 537, 1, 0, n-1)
   call repeatdouble(DD, n, d)
   write(*,'(9(f4.1))') buf
end program main

C:\gfortran\clf\repeatdouble>c:\gcc_equation\bin\x86_64-pc-mingw32-gfortran  
-c r
epeatdouble.f90

C:\gfortran\clf\repeatdouble>c:\gcc_equation\bin\x86_64-pc-mingw32-gfortran
test
3.f90 repeatdouble.o -otest3

C:\gfortran\clf\repeatdouble>test3
 repeat            4  times ...
 double    3.1415927410125732      .
 array :
 making copy nb           1  of            4

So the test failed.  I found the reason in test3.f90:

C:\gfortran\clf\repeatdouble>type test3.f90
program main
   use ISO_C_BINDING
   implicit none
   type, bind(C) :: descr
      type(C_PTR) address
      integer(C_INTPTR_T) unknown1
      integer(C_INTPTR_T) unknown2
      integer(C_INTPTR_T) stride
      integer(C_INTPTR_T) unknown3
      integer(C_INTPTR_T) last_element
   end type descr
   interface
      subroutine repeatdouble(DD,n,d) bind(C,name='repeatdouble_')
         use ISO_C_BINDING
         import descr
         implicit none
         type(descr) DD
         integer(C_INT) n
         real(C_DOUBLE) d
      end subroutine repeatdouble
   end interface
   type(descr) DD
   double precision, target :: buf(9)
   double precision d
   integer n

   buf = [(n,n=1,9)]
   d = 3.141592654
   n = 4

   DD = descr(C_LOC(buf), 0, 537, 1, 0, n-1)
! Shows the error: start
write(*,*) DD
write(*,*) C_LOC(buf)
DD%address = C_LOC(buf)
write(*,*) DD
! Shows the error: end
   call repeatdouble(DD, n, d)
   write(*,'(9(f4.1))') buf
end program main

C:\gfortran\clf\repeatdouble>c:\gcc_equation\bin\x86_64-pc-mingw32-gfortran  
-c r
epeatdouble.f90

C:\gfortran\clf\repeatdouble>c:\gcc_equation\bin\x86_64-pc-mingw32-gfortran
test
3.f90 repeatdouble.o -otest3

C:\gfortran\clf\repeatdouble>test3
                    0                    0                  537
   1                    0                    3
              2293200
              2293200                    0                  537
   1                    0                    3
 repeat            4  times ...
 double    3.1415927410125732      .
 array :  1.00000000000000000        2.0000000000000000
3.000000000000000
0        4.0000000000000000
 making copy nb           1  of            4
 making copy nb           2  of            4
 making copy nb           3  of            4
 making copy nb           4  of            4
 3.1 3.1 3.1 3.1 5.0 6.0 7.0 8.0 9.0

So it seems that the structure constructor for type(descr) isn't
assigning that first member, address.  If you do it by hand
everything works.  This is a pretty recent version of gfortran,
so I think that the above example shows a bug in the compiler.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end



Thu, 07 Oct 2010 07:14:50 GMT  
 Call Array valued Fortran function from C

Quote:
> But you should ask the gfortran folks if they can round up some
> documentation about gfortran array descriptors for you.  They
> are different in every compiler.

You could try to use the Chasm library
(<http://chasm-interop.sourceforge.net/>) which provides a more or less
consistent C API to different compiler's dope vector structures; I've
used to write interface wrappers in C which accept Fortran-pointer
arguments (passing arrays between Fortran and C was the goal).

Sebastian



Thu, 07 Oct 2010 15:44:34 GMT  
 Call Array valued Fortran function from C

Quote:
>    type, bind(C) :: descr
>       type(C_PTR) address
>       integer(C_INTPTR_T) unknown1
>       integer(C_INTPTR_T) unknown2
>       integer(C_INTPTR_T) stride
>       integer(C_INTPTR_T) unknown3
>       integer(C_INTPTR_T) last_element
>    end type descr

No, these are not the types I indicated.

Quote:
> So it seems that the structure constructor for type(descr) isn't
> assigning that first member, address.

I got to the following reduced testcase:

$ cat a.f90
program main
   use ISO_C_BINDING
   implicit none
   type, bind(C) :: descr
      type(C_PTR) :: address
   end type descr
   type(descr) :: DD
   double precision, target :: buf(1)

   buf = (/ 0 /)
   DD = descr(c_loc(buf))
   print *, transfer(DD%address, 0_c_intptr_t), &
     transfer(c_loc (buf), 0_c_intptr_t)
end program main
$ gfortran a.f90 && ./a.out
           0 -1073744200

where both gfortran and g95 give weird results (Sun doesn't want to
compile it, for some reason I don't understand).

--
FX



Thu, 07 Oct 2010 18:49:38 GMT  
 Call Array valued Fortran function from C

Quote:
>>    type, bind(C) :: descr
>>       type(C_PTR) address
>>       integer(C_INTPTR_T) unknown1
>>       integer(C_INTPTR_T) unknown2
>>       integer(C_INTPTR_T) stride
>>       integer(C_INTPTR_T) unknown3
>>       integer(C_INTPTR_T) last_element
>>    end type descr
> No, these are not the types I indicated.

I know.  These are what I came up with via reverse engineering
gfortran assembly language output.  Given documentation instead I
would be more inclined towards:

use ISO_C_BINDING, only: C_PTR, C_SIZE_T
type, bind(C) :: descriptor_dimension
   integer(C_SIZE_T) stride
   integer(C_SIZE_T) lbound
   integer(C_SIZE_T) ubound
end type descriptor_dimension
type, bind(C) :: descriptor_1
   type(C_PTR) data
   integer(C_SIZE_T) offset
   integer(C_SIZE_T) dtype
   type(descriptor_dimension) dim(1)
end type descriptor_1

The first kind type that I noticed changed from 32 to 64 bits
on going from 32-bit Windows to 64-bit Windows was C_INTPTR_T,
quite analogous to DVF/CVF/ifort int_ptr_kind().  I guess
size_t does the same thing and ssize_t follows along for the
ride which is why all the descriptor fields ended up being
64 bits wide.  It's hard to tell how they got that way from
reverse engineering alone.

Quote:
>> So it seems that the structure constructor for type(descr) isn't
>> assigning that first member, address.
> I got to the following reduced testcase:

[snip]

Quote:
> $ gfortran a.f90 && ./a.out
>           0 -1073744200
> where both gfortran and g95 give weird results (Sun doesn't want to
> compile it, for some reason I don't understand).

Thanks for taking a look at this.  I try to keep you informed when
I find this kind of behavior.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end



Thu, 07 Oct 2010 23:26:43 GMT  
 Call Array valued Fortran function from C
FX,

Quote:
> where both gfortran and g95 give weird results (Sun doesn't want to
> compile it, for some reason I don't understand).

It's the constructor that is the problem:

Replace DD = descr(c_loc(buf))
by      DD%address = c_loc(buf)

and all is well.

I have to take a look at structure constructors for other reasons
tonight, so I'll look out for this too.

Cheers

Paul



Thu, 07 Oct 2010 23:56:10 GMT  
 Call Array valued Fortran function from C

Quote:
> I have to take a look at structure constructors for other reasons
> tonight, so I'll look out for this too.

By the way, I filed it as GCC PR #35983
(http://gcc.gnu.org/bugzilla/show_bug.cgi?id=35983).

--
FX



Fri, 08 Oct 2010 00:11:16 GMT  
 Call Array valued Fortran function from C


Quote:
>          These are what I came up with via reverse engineering
> gfortran assembly language output.

Oh yeah, another thing I noticed from the assembly language code
was that if you compiled with the -fomit-frame-pointer switch the
compiler generated code that accessed at negative offsets from rsp.

I'm sure the gcc guys have been around he block on this one several
times already, but is it really safe to do this?  In SWConventions.doc
that comes with the Windows SDK, it says:

"First, all memory beyond the current address of RSP is considered
volatile:  The OS, or a de{*filter*}, may overwrite this memory during
a user debug session, or an interrupt handler.  Thus, RSP must
always be set before attempting to read or write values to a stack
frame."

The same verbiage may be found at
http://www.*-*-*.com/ (VS.80).aspx .
This is not the case in Linux where there is a 128-byte red zone
below rsp, according to Agner Fog.  Now going to
http://www.*-*-*.com/ #index-fomit_0...
it is documented that -fomit-frame-pointer makes debugging
impossible on some machines, so this covers one rationale for
Microsoft's proscription on gcc's usage, and the other rationale
seems dubious as well because the current task's context should
be saved in its TSS or some other place than its ring-3 stack,
and it would be senseless for ring-0 code such as an interrupt
handler to assume that a ring-3 stack has any space left to
safely write to.

Even so, writing to negative offsets from rsp seems to be contrary
to Microsoft's documentation, so has someone at gcc taken this up
with Microsoft and made sure that it's OK?

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end



Fri, 08 Oct 2010 11:21:14 GMT  
 
 [ 11 post ] 

 Relevant Pages 

1. VC++ calling fortran function and fortran function calling a c++ function

2. A question regarding array-valued functions in fortran 90

3. FORTRAN dynamic arrays used as parameters of functions called from C++

4. Passing array valued functions as argument to function.

5. Copying one array value into subsequent array values

6. array valued functions

7. array valued functions

8. array valued functions

9. Memory leaks with array-valued functions

10. Array valued function and memory

11. rank-reduced array-valued function (like SUM)

12. Allocatable array as return value in function

 

 
Powered by phpBB® Forum Software