Многомерный массив Fortran в C++

Question

Я пытаюсь передать многомерный массив Fortran в программу на C++ в программе взаимодействия C++ Fortran. У меня есть общее представление о том, как передаются массивы из 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. Теперь я пытаюсь передать многомерные массивы, такие как A (2 , 2) или A (2,3,2) от Fortran до C++. Я знаю, что 2-мерный массив, такой как A (2,2), будет легко определить сплющенный массив в C++. Но я считаю, что было бы немного сложнее найти элементы на C++ для 3 или 4-мерных массивов.

Каков правильный способ построения многомерного массива в C++, чтобы я мог ссылаться на элементы в массиве A (k, j, i) в Fortran как A [i] [j] [k] в C++?

Заранее благодарен!

Answer 1

Кастинг указатель на скаляр (интермедиат * в приведенном ниже примере) на указатель (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; 
} 

Compile:

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

EDIT: Хотя это может быть не по теме, я также пробовал передать выделенный или массивный указатель на ту же процедуру 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.

EDIT 2: Некоторые подпрограммы могут захотеть получить N-мерный массив в виде массива указателей. В этом случае более традиционный подход будет что-то вроде этого:

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 
}