Fortran multidimensional array in C++

Question

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!

As I can see easiest way is to pass it as one dimesianal array (as in your example). You can access elements like array[j * n + i], where i and j are indices, n - number of rows. (It is not array[i * m + j], where m is number of columns, because in Fortran multidimensional arrays stored in column major order) You can also generalize it to 3 or more dimensions. — Nikolay K
– Nikolay K, Commented May 9, 2015 at 17:38
Which Fortran compiler are you using? E.g., with the current gfortran you may (depending on the Fortran version of your source) find that arrays are represented as something more complicated than you're expecting. See the stuff around the GFC_ARRAY_DESCRIPTOR macro in libgfortran.h in the gcc sources for details. — bg2b
– bg2b, Commented May 9, 2015 at 19:59
I was actually thinking about something more elegant, rather than just simplest. I think it would be simple to just use node major order to find elements in multidimensional arrays, but it isn't always the most appeasing structure. If it comes that I may not be able to find a more pleasing structure, I will probably stick with what works the best. — J. Brad Maeng
– J. Brad Maeng, Commented May 10, 2015 at 3:50
From what I understand, Fortran arrays, no matter what dimension you have, is referred to by a pointer to the first element. So that would mean when I pass it to C++, it would only require the pointer and the size of array. bg2b, are you suggesting that what I understand is not correct? Is there more to arrays than what I had just described? — J. Brad Maeng
– J. Brad Maeng, Commented May 10, 2015 at 3:53
@J.BradMaeng that used to be true in Fortran 77, from Fortran 90 onwords (it's year 2015! even 90 is now old!) it is not true. However, when you use external procedures which the compiler does not know about it will make a copy of the array that can be referred to with a single pointer, if necessary. — Vladimir F Героям слава
– Vladimir F Героям слава, Commented May 10, 2015 at 11:04

roygvib · Accepted Answer · 2015-05-12 04:02:36Z

8

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
}

edited May 12, 2015 at 4:02

answered May 10, 2015 at 10:45

roygvib

7,3952 gold badges23 silver badges37 bronze badges

Sign up to request clarification or add additional context in comments.

5 Comments

J. Brad Maeng Over a year ago

Hi roygbvib, thank you so much for such organized and complete answer to my question. I have tried your suggestions to higher multidimensional array test problems on my own and it seems to be working flawlessly. I do appreciate your extra answers on allocatable arrays, as I will be extensively using allocatable arrays in Fortran. I will make sure to put your suggestions into good uses!. As you might have guessed, I am not a good C++ programmer yet, but I strive to learn. Thanks again a bunch!

roygvib Over a year ago

I'm glad to be of some help :D Actually, I once made an effort to make an efficient array class for my use, but ironically, the above casting turned out to be most efficient. I would appreciate if you let us know any potential issues in practical use (because my experience is very limited...)

J. Brad Maeng Over a year ago

I am an Aerospace engineering PhD student, writing a Euler solver. Euler solver is a computer program to solve numerically a fundamental system of governing equations (called Euler equations) that describe ideal fluid flow where any diffusion terms are neglected. It is a subset of greater governing equation called Navier-Stoke's. Anyhow, what I am trying is relevant because I wrote all my code in Fortran, but a fellow researcher's code is in C++. So it is essential that we can coordinate the notations so that no error occurs while transferring data between codes.

roygvib Over a year ago

Then I guess we might need to be careful when the libraries from other colleagues define arrays of pointers for 2,3,...,N-1 dimensions, as typically done for multidimensional arrays in C/C++. If the data passing is always done via 1-dim array, there should be no problem; but if the libraries recieve double*** etc, we may need to be very careful... Hope things will go successfully with your case!

roygvib Over a year ago

@J.BradMaeng I have also made an example of more traditional approach for making a multi-dimensional array (in the EDIT2 above). I believe that there are a lot of good web articles about this so please check it out :) (because my example may include bugs...)

Collectives™ on Stack Overflow

Fortran multidimensional array in C++

1 Answer 1

5 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

5 Comments

Your Answer

Sign up or log in

Post as a guest

Related