title: CUDA: Filling a column-major matrix
body: <p>I am fairly new to CUDA, and I am trying to offload to the GPU some cumbersome computations I am doing for a performance-critical project. On my computer I have two NVS 510 Graphic cards, but I am currently experimenting with one only.</p>
  
  <p>I have some big column-major matrix (1000-5000 rows x 1-5 M columns) to be filled. I was so far able to write the code  to fill the matrix like it were an array, and it works well for matrices of relatively small size.</p>
  
  <pre><code>__global__ void interp_kernel(fl_type * d_matrix, fl_type* weights, [other params], 
  int n_rows, int num_cols) {
     int index = blockIdx.x * blockDim.x + threadIdx.x;
     int column = index / n_rows;
     int row = index % n_rows;
     if (row &gt; n_sim || column &gt; num_cols) return;
     d_matrix[index] = …something(row, column,[other params]);
  }
  </code></pre>
  
  <p>The kernel is called:</p>
  
  <pre><code>fl_type *res;
  cudaMalloc((void**)&amp;res, n_columns*n_rows*fl_size);
  int block_size = 1024;
  int num_blocks = (n_rows* n_columns + block_size - 1) / block_size;
  std::cout &lt;&lt; "num_blocks:" &lt;&lt; num_blocks &lt;&lt; std::endl;
  interp_kernel &lt;&lt; &lt; num_blocks, block_size &gt;&gt; &gt; (res,[other params], n_rows,n_columns);
  </code></pre>
  
  <p>and everything works just fine.
  If I change the kernel to work with 2D threads:</p>
  
  <pre><code>__global__ void interp_kernel2D(fl_type * d_matrix, fl_type* weights, [other params], 
  int n_rows, int num_cols) {
  int column = blockIdx.x * blockDim.x + threadIdx.x;
  int row = blockIdx.y * blockDim.y + threadIdx.y;
  int index = column* n_rows + row;
  if (row &gt; n_rows || column &gt; num_cols) return;
     d_matrix[index] = …something(row, column,[other params]);
  }
  </code></pre>
  
  <p>and I invoke it</p>
  
  <pre><code>int block_size2 = 32; //each block will have block_size2*block_size2 threads
  dim3 num_blocks2(block_size2, block_size2);
  int x_grid = (n_columns + block_size2 - 1) / block_size2;
  int y_grid = (n_rows + block_size2 - 1) / block_size2;
  dim3 grid_size2(x_grid, y_grid);
  interp_kernel2D &lt;&lt;&lt; grid_size2, num_blocks2 &gt;&gt;&gt; (res,[other params], n_rows,n_columns);
  </code></pre>
  
  <p>the results are all zero and CUDA returns unknown error. What am I missing? the actual code, which compiles without error with VS2015 and CUDA 8.0, can be found here: <a href="https://pastebin.com/XBCVC7VV" rel="nofollow noreferrer">https://pastebin.com/XBCVC7VV</a> </p>
  
  <p>Here is the code from the pastebin link:</p>
  
  <pre><code>#include "cuda_runtime.h"
  #include "device_launch_parameters.h"
  #include &lt;stdio.h&gt;
  #include &lt;assert.h&gt;
  #include &lt;iostream&gt;
  #include &lt;random&gt;
  #include &lt;chrono&gt;
  typedef float fl_type;
  typedef int pos_type;
  typedef std::chrono::milliseconds ms;
  //declaration of the cuda function
  void cuda_interpolation_function(fl_type* interp_value_back, int result_size, fl_type * grid_values, int grid_values_size, fl_type* weights, pos_type* node_map, int  total_action_number, int  interp_dim, int n_sim);
  
  fl_type iterp_cpu(fl_type* weights, pos_type* node_map, fl_type* grid_values, int&amp; row, int&amp; column, int&amp; interp_dim, int&amp; n_sim) {
      int w_p = column*interp_dim;
      fl_type res = weights[w_p] * grid_values[row + node_map[w_p] * n_sim];
      for (int inter_point = 1; inter_point &lt; interp_dim; inter_point++) {
          res += weights[w_p + inter_point] * grid_values[node_map[w_p + inter_point] * n_sim + row];
      }
      return res;
  }
  
  
  __global__ void interp_kernel(fl_type * d_matrix, fl_type* weights, pos_type* node_map, fl_type* grid_values, int interp_dim, int n_sim, int num_cols) {
      int index = blockIdx.x * blockDim.x + threadIdx.x;
      int column = index / n_sim;
      int row = index % n_sim;
      int w_p = column*interp_dim;
      if (row &gt; n_sim || column &gt; num_cols) return;
      fl_type res = weights[w_p] * grid_values[row + node_map[w_p] * n_sim];
      for (int inter_point = 1; inter_point &lt; interp_dim; inter_point++) {
          res += weights[w_p + inter_point] * grid_values[row + node_map[w_p + inter_point] * n_sim];
      }
      d_matrix[index] = res;
  }
  
  __global__ void interp_kernel2D(fl_type * d_matrix, fl_type* weights, pos_type* node_map, fl_type* grid_values, int interp_dim, int n_sim, int num_cols) {
      int column = blockIdx.x * blockDim.x + threadIdx.x;
      int row = blockIdx.y * blockDim.y + threadIdx.y;
      int index = column*n_sim + row;
      int w_p = column*interp_dim;
      if (row &gt; n_sim || column &gt; num_cols) return;
      fl_type res = weights[w_p] * grid_values[row + node_map[w_p] * n_sim];
      for (int inter_point = 1; inter_point &lt; interp_dim; inter_point++) {
          res += weights[w_p + inter_point] * grid_values[row + node_map[w_p + inter_point] * n_sim];
      }
      d_matrix[index] = res;
  }
  
  void verify(fl_type *host, fl_type *device, int size) {
      int count = 0;
      int count_zero = 0;
      for (int i = 0; i &lt; size; i++) {
          if (host[i] != device[i]) {
              count++;
              //std::cout &lt;&lt;"pos: " &lt;&lt;i&lt;&lt; " CPU:" &lt;&lt;h[i] &lt;&lt; ",        GPU: " &lt;&lt; d[i] &lt;&lt;std::endl;
              assert(host[i] == device[i]);
              if (device[i] == 0.0)
                  count_zero++;
          }
      }
      if (count) {
          std::cout &lt;&lt; "Non matching: " &lt;&lt; count &lt;&lt; "out of " &lt;&lt; size &lt;&lt; "(" &lt;&lt; (float(count) / size * 100) &lt;&lt; "%)" &lt;&lt; std::endl;
          std::cout &lt;&lt; "Zeros returned from the device: " &lt;&lt; count_zero &lt;&lt;"(" &lt;&lt; (float(count_zero) / size * 100) &lt;&lt; "%)" &lt;&lt; std::endl;
      }
      else
          std::cout &lt;&lt; "Perfect match!" &lt;&lt; std::endl;
  }
  
  int main() {
      int fl_size = sizeof(fl_type);
      int pos_size = sizeof(pos_type);
      int dim = 5;             // range: 2-5
      int number_nodes = 5500; // range: 10.000-500.000
      int max_actions = 12;    // range: 6-200
      int n_sim = 1000;        // range: 1.000-10.000
      int interp_dim = std::pow(2, dim);
      int grid_values_size = n_sim*number_nodes;
      std::default_random_engine generator;
      std::normal_distribution&lt;fl_type&gt; normal_dist(0.0, 1);
      std::uniform_int_distribution&lt;&gt; uniform_dist(0, number_nodes - 1);
  
      double bit_allocated = 0;
      fl_type * grid_values;  //flattened 2d array, containing the value of the grid (n_sims x number_nodes)
      grid_values = (fl_type *)malloc(grid_values_size * fl_size);
      bit_allocated += grid_values_size * fl_size;
      for (int i = 0; i &lt; grid_values_size; i++)
          grid_values[i] = normal_dist(generator);
  
      pos_type * map_node2values_start; //vector that maps each node to the first column of the result matrix regarding that done
      pos_type * map_node2values_how_many; //vector that stores how many action we have per node  
      map_node2values_start = (pos_type *)malloc(number_nodes * pos_size);
      map_node2values_how_many = (pos_type *)malloc(number_nodes * pos_size);
  
  
      bit_allocated += 2 * (number_nodes * pos_size);
      for (int i = 0; i &lt; number_nodes; i++) {
          //each node as simply max_actions
          map_node2values_start[i] = max_actions*i;
          map_node2values_how_many[i] = max_actions;
      }
  
      //total number of actions, which is amount of column of the results
      int total_action_number = map_node2values_start[number_nodes - 1] + map_node2values_how_many[number_nodes - 1];
  
      //vector that keep tracks of the columnt to grab, and their weight in the interpolation
      fl_type* weights;
      pos_type * node_map;
      weights = (fl_type *)malloc(total_action_number*interp_dim * pos_size);
      bit_allocated += total_action_number * fl_size;
      node_map = (pos_type *)malloc(total_action_number*interp_dim * pos_size);
      bit_allocated += total_action_number * pos_size;
  
      //filling with random numbers
      for (int i = 0; i &lt; total_action_number*interp_dim; i++) {
          node_map[i] = uniform_dist(generator);      // picking random column
          weights[i] = 1.0 / interp_dim;              // uniform weights
      }
      std::cout &lt;&lt; "done filling!" &lt;&lt; std::endl;
      std::cout &lt;&lt; bit_allocated / 8 / 1024 / 1024 &lt;&lt; "MB allocated" &lt;&lt; std::endl;
  
      int result_size = n_sim*total_action_number;
      fl_type *interp_value_cpu;
      bit_allocated += result_size* fl_size;
  
  
  
      interp_value_cpu = (fl_type *)malloc(result_size* fl_size);
  
      auto start = std::chrono::steady_clock::now();
      for (int row = 0; row &lt; n_sim; row++) {
          for (int column = 0; column &lt; total_action_number; column++) {
              auto zz = iterp_cpu(weights, node_map, grid_values, row, column, interp_dim, n_sim);
              interp_value_cpu[column*n_sim + row] = zz;
          }
      }
      auto elapsed_cpu = std::chrono::steady_clock::now() - start;
      std::cout &lt;&lt; "Crunching values on the CPU (serial): " &lt;&lt; std::chrono::duration_cast&lt;ms&gt;(elapsed_cpu).count() / 1000.0 &lt;&lt; "s" &lt;&lt; std::endl;
      int * pp;
      cudaMalloc((void**)&amp;pp, sizeof(int)); //initializing the device, to not affect the benchmark
      fl_type *interp_value_gpu;
      interp_value_gpu = (fl_type *)malloc(result_size* fl_size);
      start = std::chrono::steady_clock::now();
      cuda_interpolation_function(interp_value_gpu, result_size, grid_values, grid_values_size, weights, node_map, total_action_number, interp_dim, n_sim);
      auto elapsed_gpu = std::chrono::steady_clock::now() - start;
      std::cout &lt;&lt; "Crunching values on the GPU: " &lt;&lt; std::chrono::duration_cast&lt;ms&gt;(elapsed_gpu).count() / 1000.0 &lt;&lt; "s" &lt;&lt; std::endl;
      float ms_cpu = std::chrono::duration_cast&lt;ms&gt;(elapsed_cpu).count();
      float ms_gpu = std::chrono::duration_cast&lt;ms&gt;(elapsed_gpu).count();
      int n_proc = 4;
      std::cout &lt;&lt; "Performance: " &lt;&lt; (ms_gpu- ms_cpu / n_proc) / (ms_cpu / n_proc) * 100 &lt;&lt; " % less time than parallel CPU!" &lt;&lt; std::endl;
      verify(interp_value_cpu, interp_value_gpu, result_size);
  
      free(interp_value_cpu);
      free(interp_value_gpu);
      free(grid_values);
      free(node_map);
      free(weights);
  }
  
  void cuda_interpolation_function(fl_type* interp_value_gpu, int result_size, fl_type * grid_values, int grid_values_size, fl_type* weights, pos_type* node_map, int total_action_number, int interp_dim, int n_sim) {
      int fl_size = sizeof(fl_type);
      int pos_size = sizeof(pos_type);
      auto start = std::chrono::steady_clock::now();
      //device versions of the inputs
      fl_type * grid_values_device;
      fl_type* weights_device;
      pos_type * node_map_device;
      fl_type *interp_value_device;
      int lenght_node_map = interp_dim*total_action_number;
      std::cout &lt;&lt; "size grid_values: " &lt;&lt; grid_values_size &lt;&lt;std::endl;
      std::cout &lt;&lt; "size weights: " &lt;&lt; lenght_node_map &lt;&lt; std::endl;
      std::cout &lt;&lt; "size interp_value: " &lt;&lt; result_size &lt;&lt; std::endl;
  
      //allocating and moving to the GPU the inputs
      auto error_code=cudaMalloc((void**)&amp;grid_values_device, grid_values_size*fl_size);
      if (error_code != cudaSuccess) {
          std::cout &lt;&lt; "Error during cudaMalloc of the grid_values" &lt;&lt; std::endl;
      }
      error_code=cudaMemcpy(grid_values_device, grid_values, grid_values_size*fl_size, cudaMemcpyHostToDevice);
      if (error_code != cudaSuccess) {
          std::cout &lt;&lt; "Error during cudaMemcpy of the grid_values" &lt;&lt; std::endl;
      }
      error_code=cudaMalloc((void**)&amp;weights_device, lenght_node_map*fl_size);
      if (error_code != cudaSuccess) {
          std::cout &lt;&lt; "Error during cudaMalloc of the weights" &lt;&lt; std::endl;
      }
      error_code=cudaMemcpy(weights_device, weights, lenght_node_map*fl_size, cudaMemcpyHostToDevice);
      if (error_code != cudaSuccess) {
          std::cout &lt;&lt; "Error during cudaMemcpy of the weights" &lt;&lt; std::endl;
      }
      error_code=cudaMalloc((void**)&amp;node_map_device, lenght_node_map*pos_size);
      if (error_code != cudaSuccess) {
          std::cout &lt;&lt; "Error during cudaMalloc of node_map" &lt;&lt; std::endl;
      }
      error_code=cudaMemcpy(node_map_device, node_map, lenght_node_map*pos_size, cudaMemcpyHostToDevice);
      if (error_code != cudaSuccess) {
          std::cout &lt;&lt; "Error during cudaMemcpy of node_map" &lt;&lt; std::endl;
      }
      error_code=cudaMalloc((void**)&amp;interp_value_device, result_size*fl_size);
      if (error_code != cudaSuccess) {
          std::cout &lt;&lt; "Error during cudaMalloc of interp_value_device " &lt;&lt; std::endl;
      }
      auto elapsed_moving = std::chrono::steady_clock::now() - start;
      float ms_moving = std::chrono::duration_cast&lt;ms&gt;(elapsed_moving).count();
      cudaDeviceSynchronize();
      //1d
      int block_size = 1024;
      int num_blocks = (result_size + block_size - 1) / block_size;
      std::cout &lt;&lt; "num_blocks:" &lt;&lt; num_blocks &lt;&lt; std::endl;
      interp_kernel &lt;&lt; &lt; num_blocks, block_size &gt;&gt; &gt; (interp_value_device, weights_device, node_map_device, grid_values_device, interp_dim, n_sim, total_action_number);
  
  
      //2d
      //int block_size2 = 32; //each block will have block_size2*block_size2 threads
      //dim3 num_blocks2(block_size2, block_size2);
      //int x_grid = (total_action_number + block_size2 - 1) / block_size2;
      //int y_grid = (n_sim + block_size2 - 1) / block_size2;
      //dim3 grid_size2(x_grid, y_grid);
      //std::cout &lt;&lt;"grid:"&lt;&lt; x_grid&lt;&lt;" x "&lt;&lt; y_grid&lt;&lt;std::endl;
      //interp_kernel2D &lt;&lt;&lt; grid_size2, num_blocks2 &gt;&gt;&gt; (interp_value_device, weights_device, node_map_device, grid_values_device, interp_dim, n_sim, total_action_number);
  
  
      cudaDeviceSynchronize();
      cudaError err = cudaGetLastError();
      if (cudaSuccess != err)
      {
          std::cout &lt;&lt; "Cuda kernel failed! " &lt;&lt; cudaGetErrorString(err) &lt;&lt;std::endl;
      }
      start = std::chrono::steady_clock::now();
      cudaMemcpy(interp_value_gpu, interp_value_device, result_size*fl_size, cudaMemcpyDeviceToHost);
      auto elapsed_moving_back = std::chrono::steady_clock::now() - start;
      float ms_moving_back = std::chrono::duration_cast&lt;ms&gt;(elapsed_moving_back).count();
  
      std::cout &lt;&lt; "Time spent moving the data to the GPU:" &lt;&lt; ms_moving &lt;&lt; " ms"&lt;&lt;std::endl;
      std::cout &lt;&lt; "Time spent moving the results back to the host: " &lt;&lt; ms_moving_back &lt;&lt; " ms" &lt;&lt; std::endl;
  
      cudaFree(interp_value_device);
      cudaFree(weights_device);
      cudaFree(node_map_device);
      cudaFree(grid_values_device);
  }
  </code></pre>
  
  <p>Moreover, I would extremely grateful for any direction on how to improve the performance of the code.</p>