C++ 中的 Fortran 多维数组

Fortran multidimensional array in C++

我试图在 C++ Fortran 互操作程序中将多维 Fortran 数组传递给 C++ 程序。我对将数组从 Fortran 传递到 C++ 的工作方式有一个基本的了解;您将数组的位置从 Fortran 传递给 C++。然后 C++ 采用展平数组,您必须进行一些代数计算才能在给定的多维数组中找到元素。

我能够在标量数组上成功地测试这个想法。在 C++ 中找出元素的索引并不难,因为它是从 Fortran 索引到 C++ 的线性映射,偏移量为 -1。 Fortran 和 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;
    }
};

并且输出将是 array(1) = 1 和 array(2) = 3。现在我正在尝试从中传递多维数组,例如 A(2,2) 或 A(2,3,2) Fortran 到 C++。我知道像 A(2,2) 这样的二维数组很容易在 C++ 中找出扁平数组。但我认为在 C++ 中为 3 或 4 维数组定位元素可能更具挑战性。

在 C++ 中构建多维数组的正确方法是什么,以便我可以在 Fortran 中将数组 A(k,j,i) 中的元素引用为 C++ 中的 A[i][j][k]?

提前致谢!

将标量指针(下例中的 int*)转换为指向 (N-1)-dim 数组的指针可能会有用,尽管编写数组视图 class 应该更灵活。 .

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

编译:

$ g++ fortsub.f90 cppfunc.cpp -lgfortran  # tested with gcc >=4.4.7

编辑: 虽然这可能是题外话,但我也尝试将可分配数组或数组指针传递给同一个 cppfunc() 例程。具体来说,我将上面的 a(4,3,2) 的声明更改为以下之一:

!! 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, :, :, : )                                                           

编译时使用

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

所有情况都给出了正确的结果。仅在情况 4 中,编译器创建一个临时数组,将其第一个元素的地址传递给 cppfunc(),并将获得的数据复制回实际参数。否则,编译器会将实际数组的第一个元素的地址直接传递给 cppfunc(),就像在 Fortran77 中一样。

编辑 2:某些例程可能希望接收 N-dim 数组作为指针数组。在这种情况下,更传统的方法是这样的:

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

已修改 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
}