1

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);*/

}

これが不可能な場合は、別のアプローチに移行できますが、正確性について確信が持てないため、決定できません。どんな助けでも大歓迎です。

4

0 に答える 0