CODEWITHSUNDEEP
c++arraysmultidimensional-arrayfortranfortran90
J. Brad Maeng
J. Brad Maeng

Reputation: 33

Fortran multidimensional array in C++

I am trying to be able to pass a multidimensional Fortran array to a C++ program, in a C++ Fortran interoperating program. I have a basic idea of how passing the arrays from Fortran to C++ works; you pass a location of the array from Fortran to C++. Then C++ takes the flattened array and you have to do some algebraic calculation to find the element in a given multidimensional array.

I was able to successfully test this idea on a scalar array. It is not so hard to figure out the index of element in C++, because it is linearly mapped from Fortran index to C++ with offset of -1. Sample codes for Fortran and C++ are:

! 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 such as A(2,2) or A(2,3,2) from Fortran to C++. I know that 2 dimensional array such as A(2,2) will be easy to figure out the flattened array in C++. But I reckon it might be a bit more challenging to locate 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 elements in array A(k,j,i) in Fortran as A[i][j][k] in C++?

Thanks in advance!

Upvotes: 3

Views: 4023

Answers (1)

roygvib
roygvib

Reputation: 7385

Casting a pointer to scalar (int* in the example below) to a pointer to (N-1)-dim array may be useful, although writing an array view 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: Although this may be off-topic, I also tried passing an allocatable array or array pointer to the same cppfunc() routine. Specifically, I changed the declaration of a(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 compiling with

$ g++ fortsub.f90 cppfunc.cpp -lgfortran -fcheck-array-temporaries

all the cases give the correct result. Only in case 4 the compiler creates an array temporary, pass the address of its first element to cppfunc(), and copies the obtained data back to the actual argument. Otherwise, the compiler passes the address of the first element of an actual array directly to cppfunc(), as in Fortran77.

EDIT 2: Some routines may want to receive an N-dim array as an array of pointers. In this case a more traditional approach will 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;
}

Modified 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
}

Upvotes: 7

Related Questions