Создание CUSP coo_matrix из прошедших массивов FORTRAN

Question

Я работаю над интеграцией решателей CUSP в существующий код FORTRAN. В качестве первого шага я просто пытаюсь передать пару целых массивов и float (real * 4 в FORTRAN) из FORTRAN, которые будут использоваться для построения, а затем распечатать матрицу CUSP формата COO.

До сих пор я был в состоянии следовать этой теме и есть все, чтобы собрать и ссылку: Unresolved references using IFORT with nvcc and CUSP

К сожалению, программа по-видимому, посылая мусор в матрицу касп и заканчивается грохот со следующей ошибкой:

$./fort_cusp_test 
testing 1 2 3 
sparse matrix <1339222572, 1339222572> with 1339222568 entries 
libc++abi.dylib: terminating with uncaught exception of type thrust::system::system_error: invalid argument 

Program received signal SIGABRT: Process abort signal. 

Backtrace for this error: 
#0 0x10ff86ff6 
#1 0x10ff86593 
#2 0x7fff8593ff19 
Abort trap: 6 

код для Фортрана и CUDA источников являются:

cusp_runner.cu

#include <stdio.h> 
#include <cusp/coo_matrix.h> 
#include <iostream> 
#include <cusp/krylov/cg.h> 
#include <cusp/print.h> 

#if defined(__cplusplus) 
extern "C" { 
#endif 

void test_coo_mat_print_(int * row_i, int * col_j, float * val_v, int n, int nnz) { 

    //wrap raw input pointers with thrust::device_ptr 
    thrust::device_ptr<int> wrapped_device_I(row_i); 
    thrust::device_ptr<int> wrapped_device_J(col_j); 
    thrust::device_ptr<float> wrapped_device_V(val_v); 

    //use array1d_view to wrap individual arrays 
    typedef typename cusp::array1d_view< thrust::device_ptr<int> > DeviceIndexArrayView; 
    typedef typename cusp::array1d_view< thrust::device_ptr<float> > DeviceValueArrayView; 

    DeviceIndexArrayView row_indices(wrapped_device_I, wrapped_device_I + n); 
    DeviceIndexArrayView column_indices(wrapped_device_J, wrapped_device_J + nnz); 
    DeviceValueArrayView values(wrapped_device_V, wrapped_device_V + nnz); 

    //combine array1d_views into coo_matrix_view 
    typedef cusp::coo_matrix_view<DeviceIndexArrayView,DeviceIndexArrayView,DeviceValueArrayView> DeviceView; 

    //construct coo_matrix_view from array1d_views 
    DeviceView A(n,n,nnz,row_indices,column_indices,values); 

    cusp::print(A); 
} 
#if defined(__cplusplus) 
} 
#endif 

fort_cusp_test.f90

program fort_cuda_test 

    implicit none 

interface 
    subroutine test_coo_mat_print_(row_i,col_j,val_v,n,nnz) bind(C) 
     use, intrinsic :: ISO_C_BINDING, ONLY: C_INT,C_FLOAT 
     implicit none 
     integer(C_INT) :: n, nnz, row_i(:), col_j(:) 
     real(C_FLOAT) :: val_v(:) 
    end subroutine test_coo_mat_print_ 
end interface 

    integer*4 n 
    integer*4 nnz 

    integer*4, target :: rowI(9),colJ(9) 
    real*4, target :: valV(9) 

    integer*4, pointer :: row_i(:) 
    integer*4, pointer :: col_j(:) 
    real*4, pointer :: val_v(:) 

    n  = 3 
    nnz = 9 
    rowI = (/ 1, 1, 1, 2, 2, 2, 3, 3, 3/) 
    colJ = (/ 1, 2, 3, 1, 2, 3, 1, 2, 3/) 
    valV = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9/) 

    row_i => rowI 
    col_j => colJ 
    val_v => valV 

    write(*,*) "testing 1 2 3" 

    call test_coo_mat_print_(row_i,col_j,val_v,n,nnz) 

end program fort_cuda_test 

Если вы хотите попробовать сами, вот мой (довольно безвкусный) Makefile:

Test: 
    nvcc -Xcompiler="-fPIC" -shared cusp_runner.cu -o cusp_runner.so -I/Developer/NVIDIA/CUDA-6.5/include/cusp 
    gfortran -c fort_cusp_test.f90 
    gfortran fort_cusp_test.o cusp_runner.so -L/Developer/NVIDIA/CUDA-6.5/lib -lcudart -o fort_cusp_test 

clean: 
    rm *.o *.so 

Пути библиотеки будет необходимо изменить соответственно.

Может ли кто-нибудь указать мне в правильном направлении, как правильно передать необходимые массивы из кода fortran?

---

с интерфейсом блоком удален и добавлением заявления печати в начале функции C можно увидеть, что массивы передается правильно, но п и NNZ вызывает проблемы. я получаю следующий результат:

$ ./fort_cusp_test 
testing 1 2 3 
n: 1509677596, nnz: 1509677592 
    i, row_i, col_j,  val_v 
    0,  1,  1, 1.0000e+00 
    1,  1,  2, 2.0000e+00 
    2,  1,  3, 3.0000e+00 
    3,  2,  1, 4.0000e+00 
    4,  2,  2, 5.0000e+00 
    5,  2,  3, 6.0000e+00 
    6,  3,  1, 7.0000e+00 
    7,  3,  2, 8.0000e+00 
    8,  3,  3, 9.0000e+00 
    9,  0, 32727, 0.0000e+00 
    ... 
    etc 
    ... 
    Program received signal SIGSEGV: Segmentation fault - invalid memory reference. 

Backtrace for this error: 
#0 0x105ce7ff6 
#1 0x105ce7593 
#2 0x7fff8593ff19 
#3 0x105c780a2 
#4 0x105c42dbc 
#5 0x105c42df4 
Segmentation fault: 11 

fort_cusp_test

interface 
     subroutine test_coo_mat_print_(row_i,col_j,val_v,n,nnz) bind(C) 
      use, intrinsic :: ISO_C_BINDING, ONLY: C_INT,C_FLOAT 
      implicit none 
      integer(C_INT),value :: n, nnz 
      integer(C_INT) :: row_i(:), col_j(:) 
      real(C_FLOAT) :: val_v(:) 
     end subroutine test_coo_mat_print_ 
    end interface 

     integer*4 n 
     integer*4 nnz 

     integer*4, target :: rowI(9),colJ(9) 
     real*4, target :: valV(9) 

     integer*4, pointer :: row_i(:) 
     integer*4, pointer :: col_j(:) 
     real*4, pointer :: val_v(:) 

     n  = 3 
     nnz = 9 
     rowI = (/ 1, 1, 1, 2, 2, 2, 3, 3, 3/) 
     colJ = (/ 1, 2, 3, 1, 2, 3, 1, 2, 3/) 
     valV = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9/) 

     row_i => rowI 
     col_j => colJ 
     val_v => valV 

     write(*,*) "testing 1 2 3" 

     call test_coo_mat_print_(rowI,colJ,valV,n,nnz) 

    end program fort_cuda_test 

cusp_runner.cu

#include <stdio.h> 
    #include <cusp/coo_matrix.h> 
    #include <iostream> 
    // #include <cusp/krylov/cg.h> 
    #include <cusp/print.h> 

    #if defined(__cplusplus) 
    extern "C" { 
    #endif 

    void test_coo_mat_print_(int * row_i, int * col_j, float * val_v, int n, int nnz) { 

     printf("n: %d, nnz: %d\n",n,nnz); 

     printf("%6s, %6s, %6s, %12s \n","i","row_i","col_j","val_v"); 
     for(int i=0;i<n;i++) { 
      printf("%6d, %6d, %6d, %12.4e\n",i,row_i[i],col_j[i],val_v[i]); 
     } 
     if (false) { 
     //wrap raw input pointers with thrust::device_ptr 
     thrust::device_ptr<int> wrapped_device_I(row_i); 
     thrust::device_ptr<int> wrapped_device_J(col_j); 
     thrust::device_ptr<float> wrapped_device_V(val_v); 

     //use array1d_view to wrap individual arrays 
     typedef typename cusp::array1d_view< thrust::device_ptr<int> > DeviceIndexArrayView; 
     typedef typename cusp::array1d_view< thrust::device_ptr<float> > DeviceValueArrayView; 

     DeviceIndexArrayView row_indices(wrapped_device_I, wrapped_device_I + n); 
     DeviceIndexArrayView column_indices(wrapped_device_J, wrapped_device_J + nnz); 
     DeviceValueArrayView values(wrapped_device_V, wrapped_device_V + nnz); 

     //combine array1d_views into coo_matrix_view 
     typedef cusp::coo_matrix_view<DeviceIndexArrayView,DeviceIndexArrayView,DeviceValueArrayView> DeviceView; 

     //construct coo_matrix_view from array1d_views 
     DeviceView A(n,n,nnz,row_indices,column_indices,values); 

     cusp::print(A); } 
    } 
    #if defined(__cplusplus) 
    } 
    #endif 

Answer 1

Есть два способа передачи аргументов из Fortran подпрограмм C: прежде всего, использовать интерфейсный блок (новый подход в современном Fortran) , а второй - не использовать интерфейсный блок (старый подход действителен даже для Fortran77).

Прежде всего, речь идет о первом методе использования блока интерфейса. Поскольку подпрограмма C ожидает получения указателей C (row_i, col_j и val_v), нам нужно передать адрес для этих переменных со стороны Fortran. Для этого нам нужно использовать звездочку (*), а не двоеточие (:) в интерфейсном блоке, как показано ниже. (Если мы используем двоеточие, то это говорит компилятору Fortran отправить адрес объектов указателя Fortran [1], что не является желаемым поведением.) Кроме того, поскольку n и nnz в подпрограмме C объявляются как значения (не указатели) , интерфейсный блок должен иметь атрибут VALUE для этих переменных, так что компилятор Fortran отправляет значения n и nnz, а не их адреса.Подводя итог, в первом подходе, подпрограммы C и Fortran выглядеть следующим образом:

Fortran routine: 
... 
interface 
    subroutine test_coo_mat_print_(row_i,col_j,val_v,n,nnz) bind(C) 
     use, intrinsic :: ISO_C_BINDING, ONLY: C_INT,C_FLOAT 
     implicit none 
     integer(C_INT) :: row_i(*), col_j(*) 
     real(C_FLOAT) :: val_v(*) 
     integer(C_INT), value :: n, nnz  !! see note [2] below also 
    end subroutine test_coo_mat_print_ 
end interface 
... 
call test_coo_mat_print_(rowI, colJ, valV, n, nnz) 

C routine: 
void test_coo_mat_print_ (int * row_i, int * col_j, float * val_v, int n, int nnz) 

---

Ниже о втором подходе без интерфейсного блока. При таком подходе сначала удалить блок интерфейса и массив указателей полностью, и изменить Fortran код следующим образом

Fortran routine: 

integer rowI(9), colJ(9), n, nnz  !! no TARGET attribute necessary 
real  valV(9) 

! ...set rowI etc as above... 

call test_coo_mat_print (rowI, colJ, valV, n, nnz) !! "_" is dropped 

и подпрограмму C следующим

void test_coo_mat_print_ (int* row_i, int* col_j, float* val_v, int* n_, int* nnz_) 
{ 
    int n = *n_, nnz = *nnz_; 

    printf("%d %d \n", n, nnz); 
    for(int k = 0; k < 9; k++) { 
     printf("%d %d %10.6f \n", row_i[ k ], col_j[ k ], val_v[ k ]); 
    } 

    // now go to thrust... 
} 

Обратите внимание, что N_ и nnz_ объявлены как указатели в процедуре C, потому что без блока интерфейса компилятор Fortran всегда отправляет адрес фактических аргументов в процедуру C. Также обратите внимание, что в приведенной выше процедуре C содержимое строки и т. Д. Печатается, чтобы убедиться, что аргументы переданы правильно. Если напечатанные значения верны, то я думаю, что проблема будет более вероятна в вызове подпрограмм (например, как передать информацию о размере, например n и nnz).

[1] Указатели Fortran, объявленные как «real, pointer :: a (:)», фактически представляют собой нечто вроде класса представления массива (на языке C++), который отличается от фактических данных. Здесь необходимо отправить адрес фактических данных, а не адрес этого объекта с видом массива. Кроме того, звездочка в блоке интерфейса (a (*)) представляет собой массив предполагаемого размера, который является старым методом передачи массива в Fortran. В этом случае адрес первого элемента массива передается, как и ожидалось.

[2] Если п и NNZ объявлены как указатели в подпрограмме C (как во втором подходе), то этот атрибут VALUE должен не быть прикреплен, потому что подпрограмма С хочет адрес фактических аргументов, а не их значения.