title: How to get "sum" of parallel arrays in cuda?
body: <p>my problem is about getting &quot;sum&quot; for some same length arrays. For example,I have a M*N(100 * 2000) length float array in all. I would like to get M(100) sum values of every N(2000) float numbers. I found two ways to do this job. One is with Cublas function in a for loop for M ,like <code>cublasSasum</code>. The other is self-written kernel function, adding numbers in loop. My problem is the speed of these two ways and how to choose between them.</p>
  <p>For Cublas method, no matter how big is N(4000~2E6), the time consuming is depending mainly on M, the loop number.</p>
  <p>For self-written kennel function, the speed varied much with N. In detail, if N is small, below 5000, it runs much faster than the Cublas way. Then the time consumption is increasing with N's increasing.</p>
  <p>N = 4000 |10000 | 40000 | 80000 | 1E6    | 2E6</p>
  <p>t = 254ms| 422ms | 1365ms| 4361ms| 5399ms | 10635ms</p>
  <p>If N is big enough, it runs much slower than Cublas way. My problem is how could I make a predition with M or N to decide which way I should use? My code might be used on different GPU device. Must I compare the speed in a parameter swept and then &quot;guess&quot; to make a choice in every GPU device, or I could inference from GPU device information?</p>
  <p>Besides, for the kernel function way，I also have problem in deciding the <code>blockSize</code> and <code>gridSize</code>. I must note here that what I concern more is speed not efficiency. Because the memory is limited. For example, if I got 8G memory. My dataformat is float in 4 bytes. N is 1E5. Then M is at most 2E4, which is smaller than the <code>MaxGridSize</code>. So I got two ways as below. I found have a bigger gridSize is always better, I don't know the reason. Is it about the usage of register number per thread? But I don't think it needs many registers per thread in this kernel function.</p>
  <p>Any suggestion or information would be appreciated. Thank you.</p>
  <p>Cublas way</p>
  <pre><code>for (int j = 0;j&lt;M;j++)
      cublasStatus = cublasSasum(cublasHandle,N,d_in+N*j,1,d_out+j);
  </code></pre>
  <p>self-written kernel way</p>
  <pre><code>__global__ void getSum(int M, int N, float* in, float * out)
  {
      int i = threadIdx.x + blockIdx.x * blockDim.x;
      if(i&lt;M){
          float tmp = 0;
          for(int ii = 0; ii&lt;N; ii++){
              tmp += *(in+N*i+ii);
          }
          out[i] = tmp;
      }
  }
  </code></pre>
  <p>Bigger gridSize is faster. I don't know the reason.</p>
  <pre><code>getSum&lt;&lt;&lt;M,1&gt;&gt;&gt;(M, N, d_in, d_out); //faster
  getSum&lt;&lt;&lt;1,M&gt;&gt;&gt;(M, N, d_in, d_out); 
  
  </code></pre>
  <p>This is a blockSize-time parameter swept result. M = 1E4.N = 1E5.</p>
  <pre><code>cudaEventRecord(start, 0);
  //blockSize = 1:1024;
  int gridSize = (M + blockSize - 1) / blockSize;
  getSum&lt;&lt;&lt;gridSize1,blockSize1&gt;&gt;&gt;...
  cudaEventRecord(stop, 0);
  cudaEventSynchronize(stop);
  cudaEventElapsedTime(&amp;time, start, stop);
  </code></pre>
  <p>It seems I should choose a relative small <code>blockSize</code>, like 10~200. I just would like to know why the full occupancy(blockSize 1024) is slower. I just post here for some possible reasons, registers number?latency?
  <a href="https://i.stack.imgur.com/pnNGq.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/pnNGq.png" alt="enter image description here" /></a></p>