title: LU Decomposition in CUDA
body: <p>I have been trying to implement LU Decomposition. This code works for matrices upto sizes <code>100x100</code> but fails for larger matrices. I have no clue as to what is happening. I changed the max finding function, but had similar results. Please help.</p>
  
  <pre><code> __device__ float max_value = 0;  
   __device__ int max_index;
   __device__ int lock = 0;
   __device__ int blockLock = 0;
  
  __global__ void maxIndex(float *LU, int n, int col)
   {
  max_value = 0;
  max_index = INT_MAX;
  
  __shared__ float inter[64];
  __shared__ int indexBlock;
  __shared__ int threadLast;
  
  indexBlock = INT_MAX;
  threadLast = INT_MAX;
  
  int index = n*(blockDim.x*blockIdx.x + threadIdx.x + col) + col;
  float myValue = fabsf(LU[index]);
  int noElem = n - col - blockIdx.x*blockDim.x &gt; blockDim.x ? blockDim.x : n - col - blockIdx.x*blockDim.x;
  
  if(index &lt; n*n)
      inter[threadIdx.x] = fabsf(LU[index]);
  else
      inter[threadIdx.x] = FLT_MIN;
  
  __syncthreads();
  
  int h = ceil(log2((float)noElem));
  
  for(int d=0; d&lt;h; d++)
  {
      if(threadIdx.x &lt; ceil(noElem/exp2((float)(d+1))))
      {
              int parent = threadIdx.x;
                  int left = 2*parent;
                  int right = 2*parent + 1;
  
                  if(right &lt; ceil(noElem/exp2((float)d)))
                      inter[parent] = (inter[left] &gt;= inter[right]) ? inter[left] : inter[right];
              else
                      inter[parent] = inter[left];
      }
      __syncthreads();
  }
  
  __syncthreads();
  
  int loop = 1;
  if(myValue == inter[0])
  {
      while(loop)
      {
          if( 0 == atomicCAS( &amp;blockLock, 0, 1 ) )
          {
              if(threadIdx.x &lt; threadLast)
              {
                  threadLast = threadIdx.x;
                  indexBlock = blockDim.x*blockIdx.x + threadIdx.x + col;
              }
              __threadfence_block();
  
              atomicExch( &amp;blockLock, 0);
              loop = 0;
          }
       }
  }   
  
  __syncthreads();
  
  if( threadIdx.x == 0 &amp;&amp; index &lt; n*n)
  {
      while(0 != atomicCAS(&amp;lock, 0, 1));
  
      if((max_value &lt; inter[0]) || (max_value == inter[0] &amp;&amp; indexBlock &lt; max_index))
      {
          max_value = inter[0];
          max_index = indexBlock;
      }
  
      __threadfence();
  
      atomicExch(&amp;lock, 0);   
  }
  }
  
  __global__ void swap(float *LU, int n, int row)
  {
  if(max_index == row)
      return;
  
  float temp;
  int index =  blockDim.x*blockIdx.x + threadIdx.x;
  
  if(index &lt; n)
  {
      temp = LU[index + n*row];
      LU[index + n*row] = LU[index + n*max_index];
      LU[index + n*max_index] = temp;
  }
  
    }
  
   __global__ void elimination(float *LU, int n, int row)
   {
  float factor;
  
  int indexX = row + 1 + blockDim.x*blockIdx.x + threadIdx.x;
  int indexY = row + 1 + blockIdx.y;
  
  
  if((indexX &lt; n) &amp;&amp; (indexY &lt; n))
  {
      factor = LU[n*indexY + row] / LU[n*row + row];
  
      LU[n*indexY + indexX] -= LU[n*row + indexX]*factor; 
  }
  
  __syncthreads();
  
  if(blockIdx.x == 0 &amp;&amp; threadIdx.x == 0)
  {
      LU[n*indexY + row] = factor;
      __threadfence();
  }
  
  }
  
  int luDecomposeP(float *LU, int n)
  {
  int i, noOfThreadsPerBlock = 64, noOfBlocks, sharedSize, pivotValue;
  float *dLU;
  
  cudaMalloc((void **)&amp;dLU, n*n*sizeof(float));
  
  cudaMemcpy(dLU, LU, n*n*sizeof(float), cudaMemcpyHostToDevice);
  
  dim3 gridDim(1,1,1);
  dim3 blockDim(noOfThreadsPerBlock,1,1);
  
  for(i=0; i&lt;n-1; i++)
  {
      noOfBlocks = ceil((float)(n-i) / (float)noOfThreadsPerBlock);
      sharedSize = ((noOfThreadsPerBlock &lt; (n-i) ? noOfThreadsPerBlock : n-i) + 3) * sizeof(float);
      maxIndex &lt;&lt;&lt; noOfBlocks, noOfThreadsPerBlock, sharedSize &gt;&gt;&gt; (dLU, n, i);
      //maxIndex &lt;&lt;&lt; 1, noOfThreadsPerBlock &gt;&gt;&gt; (dLU, n, i);
  
      cudaMemcpyFromSymbol(&amp;pivotValue, "max_value", sizeof(pivotValue), 0, cudaMemcpyDeviceToHost);      
      if(pivotValue &lt;= 1e-20F)
          return -1;
  
      noOfBlocks = ceil((float)n / (float)noOfThreadsPerBlock);
      swap &lt;&lt;&lt; noOfBlocks, noOfThreadsPerBlock &gt;&gt;&gt; ( dLU, n, i);
  
      gridDim.x = ceil((float)(n-i) / (float)noOfThreadsPerBlock);
      gridDim.y = n-i-1;
      elimination &lt;&lt;&lt; gridDim, blockDim &gt;&gt;&gt; ( dLU, n, i);
  
  }
  
  cudaMemcpy(LU, dLU, n*n*sizeof(float), cudaMemcpyDeviceToHost);
  
  cudaFree(dLU);
  
  return 0;
   }
  </code></pre>