C 中的Fortran多维阵列

Fortran multidimensional array in C++

本文关键字：阵列 Fortran 中的      更新时间：2023-10-16

我试图在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;
    }
};

，输出将为数组（1）= 1，并且数组（2）= 3。C 。我知道，诸如A（2,2）之类的2维数组将很容易找出C 中的扁平阵列。但是我认为，在3或4个维数阵列中找到C 中的元素可能会更具挑战性。

在C 中构造多维数组的正确方法是什么，以便我可以在fortran中的数组A（K，J，i）中的元素作为C 中的[I] [J] [K]？

预先感谢！

将指向标量的指针（在下面的示例中）铸造为指向（n-1）-dim阵列的指针可能有用，尽管编写数组视图类应该更灵活...

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（）例程。具体而言，我将上述（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
}
相关文章：