Генерация CUSP coo_matrix из переданных массивов FORTRAN
Я работаю над интеграцией решателей CUSP в существующий код FORTRAN. В качестве первого шага я просто пытаюсь передать пару целочисленных массивов и число с плавающей запятой (действительное * 4 в FORTRAN) из FORTRAN, которое будет использоваться для построения и последующей печати матрицы CUSP формата COO.
До сих пор я был в состоянии следовать этой теме и получил все для компиляции и ссылки: Неразрешенные ссылки с использованием IFORT с nvcc и CUSP
К сожалению, программа, очевидно, отправляет мусор в матрицу 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 и fortran следующий:
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
Если вы хотите попробовать это сами, вот мой (довольно не элегантный) make-файл:
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?
со снятым блоком интерфейса и добавлением оператора print в начале функции C я вижу, что массивы передаются правильно, но n и 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
Решение
Существует два метода передачи аргументов из процедур Fortran в C: первый — использовать блок интерфейса (новый подход в современном Fortran), а второй — не использовать блок интерфейса (старый подход, действительный даже для Fortran77).
Во-первых, следующее о первом способе использования интерфейсного блока. Поскольку подпрограмма C ожидает получения указателей C (row_i, col_j и val_v), нам нужно передать адрес для этих переменных со стороны Фортрана. Для этого мы должны использовать звездочку (*), а не двоеточие (:) в блоке интерфейса, как показано ниже. (Если мы используем двоеточие, то это говорит компилятору 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 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 содержимое row_i и т. Д. Печатается, чтобы убедиться, что аргументы переданы правильно. Если напечатанные значения верны, то я думаю, что проблема будет более вероятной при вызове подпрограмм тяги (включая способ передачи информации о размере, такой как n и nnz).
[1] Указатели Fortran, объявленные как «real, pointer :: a (:)», фактически представляют что-то вроде класса представления массива (на жаргоне C ++), который отличается от фактических указанных данных. Здесь необходимо отправить адрес фактических данных, а не адрес этого объекта просмотра массива. Кроме того, звездочка в блоке интерфейса (a (*)) представляет массив предполагаемого размера, который является старым методом передачи массива в Fortran. В этом случае адрес первого элемента массива передается, как и ожидалось. [2] Если n и nnz объявлены как указатели в подпрограмме C (как во втором подходе), то этот атрибут VALUE должен не быть присоединенным, потому что подпрограмма C хочет адрес фактических аргументов, а не их значений.Другие решения