Call Array valued Fortran function from C
koentj...@gmail.com
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://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----
glen herrmannsfeldt
> 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.
> 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)
> ---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.
> 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.
James Van Buskirk
news:80a7c031-33c1-4463...@l64g2000hse.googlegroups.com...
> 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
> [1]: http://docs.sun.com/source/806-3593/11_cfort.html
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
FX
> 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];
}
where descriptor_dimension is defined as
typedef struct descriptor_dimension
{
ssize_t stride;
ssize_t lbound;
ssize_t ubound;
}
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
James Van Buskirk
news:fudrd5$9ts$1...@nef.ens.fr...
>> 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.
Sebastian Hanigk
> 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
FX
> 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.
> 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
James Van Buskirk
>> 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.
>> 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]
> $ 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.
nf.e...@gmail.com
> 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
FX
> 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
James Van Buskirk
news:yLGdnReybdwp_ZbV...@comcast.com...
> 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 debugger, 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://msdn2.microsoft.com/en-us/library/x4ea06t0(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://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html#index-fomit_002dframe_002dpointer-592
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?