Fortran から呼び出した CUSP ソルバーを使用して、CG/BicGSTAB を繰り返し解きたいと考えています。転送を避けるために、Fortran データを CUSP に直接渡しています。コードはコンパイルされますが、実行時にフラグが立てられます。
terminate called after throwing an instance of 'thrust::system::system_error'
what(): invalid argument
terminate called recursively
Aborted (core dumped)
コードのコアは言うまでもなく、印刷ストリームも発生していません。もちろんコードは準備段階ですが、どこが悪いのか気になります。
extern "C" void bicgstab_(int *device_I, int *device_J, float *device_V, float *device_x, float *device_b, int *n, int *nnz){
int N = *n;
int NNZ = *nnz;
std::cout << N << " " << NNZ << " " << *device_I << std::endl;
for(int i=0; i<N;i++)std::cout << device_I[i] << " "; std::cout << std::endl;
for(int i=0; i<NNZ;i++)std::cout << device_J[i] << " "; std::cout << std::endl;
for(int i=0; i<NNZ;i++)std::cout << device_V[i] << " "; std::cout << std::endl;
for(int i=0; i<N;i++)std::cout << device_x[i] << " "; std::cout << std::endl;
for(int i=0; i<N;i++)std::cout << device_b[i] << " "; std::cout << std::endl;
// *NOTE* raw pointers must be wrapped with thrust::device_ptr!
thrust::device_ptr<int> wrapped_device_I(device_I);
thrust::device_ptr<int> wrapped_device_J(device_J);
thrust::device_ptr<float> wrapped_device_V(device_V);
thrust::device_ptr<float> wrapped_device_x(device_x);
thrust::device_ptr<float> wrapped_device_b(device_b);
// use array1d_view to wrap the individual arrays
typedef typename cusp::array1d_view< thrust::device_ptr<int> > DeviceIndexArrayView;
typedef typename cusp::array1d_view< thrust::device_ptr<float> > DeviceValueArrayView;
std::cout << wrapped_device_I[3];
/*
DeviceIndexArrayView row_indices (wrapped_device_I, wrapped_device_I + (N+1));
DeviceIndexArrayView column_indices(wrapped_device_J, wrapped_device_J + NNZ);
DeviceValueArrayView values (wrapped_device_V, wrapped_device_V + NNZ);
DeviceValueArrayView x (wrapped_device_x, wrapped_device_x + N);
DeviceValueArrayView b (wrapped_device_b, wrapped_device_b + N);
// std::cout << device_x[0] ;
// for(int i=0;i<NNZ;i++)std::cout << column_indices[i] << std::endl;
// combine the three array1d_views into a csr_matrix_view
typedef cusp::csr_matrix_view<DeviceIndexArrayView,
DeviceIndexArrayView,
DeviceValueArrayView> DeviceView;
// construct a csr_matrix_view from the array1d_views
DeviceView A(N, N, NNZ, row_indices, column_indices, values);
// set stopping criteria: // iteration_limit = 100 // relative_tolerance = 1e-5
cusp::verbose_monitor<float> monitor(b, 100, 1e-5);
// solve the linear system A * x = b with the Conjugate Gradient method
// cusp::krylov::bicgstab(A, x, b);*/
}
これが不可能な場合は、別のアプローチに移行できますが、正確性について確信が持てないため、決定できません。どんな助けでも大歓迎です。