1

I want to repetitively solve the CG/BicGSTAB using CUSP solver, called from Fortran. To avoid transfers I am passing the Fortran data directly to CUSP. The code compiles but breaks at the run time flagging:

terminate called after throwing an instance of 'thrust::system::system_error'
  what():  invalid argument
terminate called recursively
Aborted (core dumped)

Let alone the core of the code, even the print stream is not happening. The code of course is in the preliminary stage, but I wonder what is wrong with it.

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

}

If this is not feasible, I can move over to another approach,but as I am not sure about the correctness, I am unable to decide. Any help is appreciated.

The Hiary
  • 109
  • 1
  • 13

0 Answers0