title: CUDA - Optimize mean of matrix rows calculation using shared memory
body: <p>I am trying to optimize the computation of the mean of each row in my <code>512w x 1024h</code> image, and then subtract the mean from the row from which it was computed. I wrote a piece of code which does it in <code>1.86 ms</code>, but I want to reduce the speed. This piece of code works fine, but does not use shared memory, and it utilizes for loops. I want to do away with them. </p>
  
  <pre><code>__global__ void subtractMean (const float *__restrict__ img, float *lineImg, int height, int width) {
  
    // height = 1024, width = 512
  
    int tidy = threadIdx.x + blockDim.x * blockIdx.x; 
  
    float sum = 0.0f; 
    float sumDiv = 0.0f; 
  
    if(tidy &lt; height) { 
  
        for(int c = 0; c &lt; width; c++) { 
  
            sum += img[tidy*width + c];
        }
        sumDiv = (sum/width)/2;
  
        //__syncthreads(); 
  
        for(int cc = 0; cc &lt; width; cc++) { 
  
            lineImg[tidy*width + cc] = img[tidy*width + cc] - sumDiv;
        }
  
    }
  
    __syncthreads(); 
  </code></pre>
  
  <p>I called the above kernel using: </p>
  
  <pre><code>subtractMean &lt;&lt;&lt; 2, 512 &gt;&gt;&gt; (originalImage, rowMajorImage, actualImHeight, actualImWidth);
  </code></pre>
  
  <p>However, the following code I wrote uses shared memory to optimize. But, it does not work as expected. Any thoughts on what the problem might be?</p>
  
  <pre><code>__global__ void subtractMean (const float *__restrict__ img, float *lineImg, int height, int width) {
  
    extern __shared__ float perRow[];
  
    int idx = threadIdx.x;    // set idx along x
    int stride = width/2; 
  
    while(idx &lt; width) { 
        perRow[idx] = 0; 
        idx += stride; 
    }
  
    __syncthreads(); 
  
    int tidx = threadIdx.x;   // set idx along x
    int tidy = blockIdx.x;    // set idx along y
  
    if(tidy &lt; height) { 
        while(tidx &lt; width) { 
            perRow[tidx] = img[tidy*width + tidx];
            tidx += stride; 
        }
    }
  
    __syncthreads(); 
  
    tidx = threadIdx.x;   // reset idx along x
    tidy = blockIdx.x;    // reset idx along y
  
    if(tidy &lt; height) { 
  
        float sumAllPixelsInRow = 0.0f; 
        float sumDiv = 0.0f; 
  
        while(tidx &lt; width) { 
            sumAllPixelsInRow += perRow[tidx];
            tidx += stride;
        }
        sumDiv = (sumAllPixelsInRow/width)/2;
  
        tidx = threadIdx.x;   // reset idx along x
  
        while(tidx &lt; width) { 
  
            lineImg[tidy*width + tidx] = img[tidy*width + tidx] - sumDiv; 
            tidx += stride;
        }
    }
  
    __syncthreads();  
  }
  </code></pre>
  
  <p>The shared memory function was called using: </p>
  
  <pre><code>subtractMean &lt;&lt;&lt; 1024, 256, sizeof(float)*512 &gt;&gt;&gt; (originalImage, rowMajorImage, actualImHeight, actualImWidth);
  </code></pre>