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