Параллельная матричная инверсия в CUDA / GPU с использованием разделяемой памяти

Question

Я написал программу cuda для взятия инверсии заданной матрицы A (nxn) с помощью шагов:

1. увеличение с помощью матрицы I (nxn)
2. сокращение обеих матриц с использованиемМетод исключения Гаусса-Джордана. В следующем цикле я использовал два ядра для него

for (int i = 0; i<n; i++){ normalize << <numBlocks, threadsPerBlock >> >(d_A, dI, n, i); gaussjordon << <numBlocks, threadsPerBlock >> >(d_A, dI, n, i); }

, где i = индекс столбца

Подробности ядер:

1. normalize: переводит i-ю строку A (nxn) и I (nxn) i матрицы A и I в общую память и делит на Aii. (i-й элемент i-й строки)
2. guassjordan: перенести строки! = i в общую память и выполнить элементарную операцию над строкой, чтобы уменьшить количество элементов i-го столбца до 0.

Но мой кодработает только тогда, когда число потоков в блоке: (n, n) я хочу использовать код для больших размеров матрицы до (1000x1000).

Я думаю, что есть проблема с индексацией разделяемой памяти. (не уверен насчет этого): Код ядер следующий:

__global__ void normalize(cuDoubleComplex *A, cuDoubleComplex *I, int i, int n){

__shared__ cuDoubleComplex RA[512];
    __shared__ cuDoubleComplex RI[512];
      __shared__ cuDoubleComplex Aii;

           Aii  = A[i*n + i];

      int rowId = threadIdx.x;

    RA[rowId] = A[i*n + rowId];
    RI[rowId] = I[i*n + rowId];


  RA[rowId] = cuCdiv(RA[rowId],Aii);
  RI[rowId] = cuCdiv(RI[rowId],Aii);

         A[i*n + rowId] = RA[rowId];
         I[i*n + rowId] = RI[rowId];
}

и

__global__ void gaussjordon(cuDoubleComplex *A, cuDoubleComplex *I, int n, int i){  
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;

    int tx = threadIdx.x; int ty=threadIdx.y; 
  __shared__ cuDoubleComplex rowA[512];
  __shared__ cuDoubleComplex rowI[512];

  __shared__ cuDoubleComplex rowAj[512];
  __shared__ cuDoubleComplex rowIj[512];

          cuDoubleComplex Aii;

             rowA[tx] = A[i*n + col];
             rowI[tx] = I[i*n + col];

            rowAj[ty*n+tx] = A[row*n + col];
            rowIj[ty*n+tx] = I[row*n + col];
                      Aii  =  A[ty*n + i]; 
if(ty != i){
         if (rowA[tx].x != 0 || rowI[tx].x != 0){
            rowAj[ty*n+tx] = cuCsub(rowAj[ty*n+tx],(cuCmul(Aii,rowA[tx])));
            rowIj[ty*n+tx] = cuCsub(rowIj[ty*n+tx],(cuCmul(Aii,rowI[tx])));
         } 
             }
             A[row*n + col] = rowAj[ty*n+tx];
             I[row*n + col] = rowIj[ty*n+tx];      
        }