title: Matrix Multiplication giving wrong output
body: <p>What I am attempting to do is Multiply Matrix A &amp; Matrix B and then from the product matrix I get the index of the maximum value per column. But unfortunately, only the first 128*128 values of the matrix multiplication are correct while others are just garbage. I do not quite understand how this works. I request you to kindly guide me with this ..</p>
  
  <pre><code>#include&lt;stdio.h&gt;
  #include "cuda.h"
  #include&lt;stdlib.h&gt;
  
  #define blockD 32
  const int wA = 128;
  const int hA = 4096;    
  const int wB = 4096;
  const int hB = wA;
  
  main(void){
  
      void MatrixMultiplication(float *, float *, float *, float *);
  
      int size_A = wA * hA * sizeof(float);
      int size_B = wB * hB * sizeof(float);
      int size_C = wB * hA * sizeof(float);
      int size_max = 2 * wB * sizeof(float);
      float *M, *N, *P, *C;   
  
      // allocate memory on the CPU
      M = (float*)malloc(size_A);
      N = (float*)malloc(size_B);
      P = (float*)malloc(size_max);
      C = (float*)malloc(size_C);
  
      // initialize the matrices
      for (int y=0; y &lt; hA; y++) {
          for (int x=0; x &lt; wA; x++){
              M[y*wA + x] = 32; //x + y*wA; 
         }
      }
  
      for (int y=0; y&lt;hB; y++) {
          for (int x=0; x&lt;wB; x++){
              N[y*wB + x] = 21; //x + y*wB; 
         }
      }
  
  
      MatrixMultiplication(M, N, P, C);
  
      //Write
      FILE *f1;
      int i,j;
      f1 = fopen("C.txt","w");
      for(i = hA - 2 ; i &lt; hA; i ++){
      for(j = 0; j &lt; wB; j++){
          fprintf(f1,"%d\t",int(C[i*wB + j]));
      }
      fprintf(f1,"\n");
      }
      fclose(f1);
  
      // free the memory allocated on the CPU
      free( M );
      free( N );
      free( P ); 
      free( C );
      cudaDeviceReset();
      return 0;
  }
  
  
  __device__ void MaxFunction(float* Pd, float* max)
  {
   int x = (threadIdx.x + blockIdx.x * blockDim.x);  
   int y = (threadIdx.y + blockIdx.y * blockDim.y); 
  
   int k = 0;
  
   int temp = 0; int temp_idx = 0;
   for (k = 0; k &lt; wB; ++k) {
              if(Pd[x*wB + k] &gt; temp){
                  temp = Pd[x*wB + k];
                  temp_idx = x*wB + k;
              }
    }
    max[y*2 + 0] = temp;
    max[y*2 + 1] = temp_idx;
  }
  
  
  __global__ void MatrixMulKernel(float* Md, float* Nd, float* Pd, float* max)
  {
    // declare cache in the shared memory
    __shared__ float Mds[blockD][blockD];
    __shared__ float Nds[blockD][blockD];
  
    float Pvalue = 0;
    // Loop over the Md and Nd block dimension required to compute the Pd element
    for (int m = (wA * blockD * blockIdx.y), n = (blockD * blockIdx.x); 
                              m &lt; ((wA * blockD * blockIdx.y)+wA-1); 
                                          m += blockD, n += (blockD*hB)){
  
      // collaboratively loading of Md and Nd blocks into shared memory    
      Mds[threadIdx.y][threadIdx.x] = Md[m + wA * threadIdx.y + threadIdx.x];
      Nds[threadIdx.y][threadIdx.x] = Nd[n + wA * threadIdx.y + threadIdx.x];
      __syncthreads();
  
      // keep track of the running sum    
      for (int k = 0; k &lt; blockD; k++)
        Pvalue += Mds[threadIdx.y][k] * Nds[k][threadIdx.x];
      __syncthreads();
    }
  
    // write back to the global memory
    int p = hB * blockD * blockIdx.y + blockD * blockIdx.x;
    Pd[p + hB * threadIdx.y + threadIdx.x] = Pvalue;
    __syncthreads();
  
    MaxFunction(Pd, max);
  
  }
  
  void MatrixMultiplication(float *M, float *N, float *P, float *C) {
  
      int size_A = wA * hA * sizeof(float);
      int size_B = wB * hB * sizeof(float);
      int size_C = wB * hA * sizeof(float);
      int size_max = 2 * wB * sizeof(float);
      float *Md, *Nd, *Pd, *max; 
  
      // allocate memory on the GPU
      cudaMalloc((void**)&amp;Md, size_A);
      cudaMalloc((void**)&amp;Nd, size_B);
      cudaMalloc((void**)&amp;Pd, size_C);
      cudaMalloc((void**)&amp;max, size_max);
  
      // transfer M and N to device memory
      cudaMemcpy(Md, M, size_A, cudaMemcpyHostToDevice);
      cudaMemcpy(Nd, N, size_B, cudaMemcpyHostToDevice);
  
      // kernel invocation code
      dim3 dimBlock(blockD, blockD);
      dim3 dimGrid(wA/blockD, hB/blockD);
  
      //Execute Kernel
      MatrixMulKernel&lt;&lt;&lt;dimGrid, dimBlock&gt;&gt;&gt;( Md, Nd, Pd, max);
  
      // transfer P from device    
      cudaMemcpy(P, max, size_max, cudaMemcpyDeviceToHost);
      cudaMemcpy(C, Pd, size_C, cudaMemcpyDeviceToHost);
  
      // free the memory allocated on the GPU
      cudaFree(Md);
      cudaFree(Nd);
      cudaFree(Pd);
      cudaFree(max);
  }
  </code></pre>