Fortran multidimensional array in C ++
I am trying to pass a multidimensional Fortran array to a C ++ program in a C ++ Fortran interop program. I have a basic understanding of how arrays are passed from Fortran to C ++; you are passing the location of the array from Fortran to C ++. Then C ++ takes a flattened array and you need to do some algebraic calculation to find an element in a given multidimensional array.
I was able to successfully test this idea on a scalar array. It is not that hard to compute the index of an element in C ++ because it is mapped linearly from the Fortran index in C ++ with an offset of -1. Sample codes for Fortran and C ++:
! Fortran main program
program fprogram
integer :: i
real*8 :: array(2)
array(1) = 1.0
array(2) = 2.0
! call cpp function
call cppfuncarray(array, 2)
write(*,*) array
end program
//cpp file
extern "C" {
void cppfuncarray_(double *array, int *imax);}
void cppfuncarray_(double *array, int *imax) {
int iMax = *imax;
for ( int i = 0; i < iMax; i++ ) {
array[i] = array + 1.0*i;
}
};
and the output will be array (1) = 1 and array (2) = 3. Now I am trying to pass multidimensional arrays like A (2,2) or (2,3,2) from Fortran to C ++. I know that a 2 dimensional array like A (2,2) would be easy to define a flattened array in C ++. But I believe it would be a little more difficult to find elements in C ++ for 3 or 4 dimensional arrays.
What is the correct way to construct a multidimensional array in C ++ so that I can refer to the elements in the array A (k, j, i) in Fortran as A [i] [j] [k] in C ++?
Thanks in advance!
source to share
Casting a pointer to a scalar (int * in the example below) into a pointer to an (N-1) -dim-array can be useful, although notation of the array representation class should be more flexible ...
fortsub.f90:
subroutine fortsub()
implicit none
integer a(4,3,2) !! this line will be changed in the EDIT below
integer ndims(3), i
ndims(:) = [ ( size( a, i ), i = 1, 3 ) ]
call cppfunc( a, ndims )
print *, "a(1,1,1) = ", a(1,1,1) !! gives 10101
print *, "a(2,2,1) = ", a(2,2,1) !! gives 20201
print *, "a(4,3,2) = ", a(4,3,2) !! gives 40302
end subroutine
cppfunc.cpp:
extern "C" {
void fortsub_( void );
void cppfunc_( int *A, int *n )
{
typedef int (*A3d_t)[ n[1] ][ n[0] ];
A3d_t A3d = (A3d_t) A; // get a pointer to 2-dim subarray
// 3-dim access
for( int k = 0; k < n[2]; k++ )
for( int j = 0; j < n[1]; j++ )
for( int i = 0; i < n[0]; i++ ) {
A3d[ k ][ j ][ i ] = (i+1)*10000 + (j+1)*100 + (k+1); // set test data
}
}
}; // extern "C"
int main()
{
fortsub_();
return 0;
}
Compile:
$ g++ fortsub.f90 cppfunc.cpp -lgfortran # tested with gcc >=4.4.7
EDIT: While it may be off-topic, I've also tried passing a dedicated or massive pointer to the same cppfunc () routine. Specifically, I changed the declaration (4,3,2) above to one of the following:
!! case 1
integer, allocatable :: a(:,:,:)
allocate( a(4,3,2) )
!! case 2
integer, pointer :: a(:,:,:)
allocate( a(4,3,2) )
!! case 3 : an array view for contiguous memory
integer, target :: b(4,3,2,5)
integer, pointer :: a(:,:,:)
a => b( :, :, :, 5 )
!! case 4 : an array view for non-contiguous memory
integer, target :: c(5,4,3,2)
integer, pointer :: a(:,:,:)
a => c( 5, :, :, : )
When compiled with
$ g++ fortsub.f90 cppfunc.cpp -lgfortran -fcheck-array-temporaries
all cases give the correct result. Only in case 4, the compiler creates a temporary array, passes the address of its first element to cppfunc (), and copies the resulting data back into the actual argument. Otherwise, the compiler passes the address of the first element of the actual array directly to cppfunc (), as in Fortran77.
EDIT 2: Some routines might want to get an N-dimensional array as an array of pointers. In this case, a more traditional approach would be something like this:
getptr3d.hpp:
template <typename T>
T*** get_ptr3d( T* A, int* n )
{
typedef T (*A3d_t)[ n[1] ][ n[0] ];
A3d_t A3d = (A3d_t) A;
T*** p = new T** [ n[2] ];
for( int k = 0; k < n[2]; k++ ) {
p[ k ] = new T* [ n[1] ];
for ( int j = 0; j < n[1]; j++ ) {
p[ k ][ j ] = (T*) &( A3d[ k ][ j ][ 0 ] );
}
}
return p;
}
template <typename T>
void free_ptr3d( T*** p, int*n )
{
for( int k = 0; k < n[2]; k++ ) { delete[] p[ k ]; }
delete[] p;
}
Changed cppfunc.cpp:
#include "getptr3d.hpp"
...
void cppfunc_( int* A, int* n )
{
int*** A3d = get_ptr3d( A, n ); // can also be used for double***
... // use A3d[ k ][ j ][ i ]
// or pass A3d to other functions like func( int*** B, ... )
free_ptr3d( A3d, n ); // when calculation finished
}
source to share