通过 cuda 实现快速并行矩阵求逆算法但有些东西不起作用

Question

目前，我正在实施

用于 cuda 的下三角矩阵。我现在只使用一个线程在矩阵的左下角进行矩阵乘法。但是，我在代码中使用的临时矩阵似乎存在一些问题。我怀疑第24行存在临时矩阵无法存储矩阵乘法后的值的问题。这里是result。有谁知道如何解决这个问题？

void print_mat(int size_small,int size_large, float*m){
    for (int row = 0; row < size_small; row++){
        for(int col = 0; col<size_small;col++){
            printf("%.2f ",m[ col + row*size_large]);

            if (col % size_small == size_small-1)
                printf("\n");

        }
    }
}


__device__ void SMM(float *A1inv, float *A2, float *A3inv, float*temp, float * O,int small_size, int large_size){

    for (int i = 0; i  < small_size; i++){
        for (int j = 0; j  < small_size; j++){
            temp[small_size*i + j] = 0;
            for (int k = 0; k < small_size; k++){
                temp[small_size*i + j] += A3inv[i*large_size+k]*A2[k*small_size+j];
            }
            temp[small_size*i+j] = -temp[small_size*i+j];
        }
    }

    for (int i = 0; i  < small_size; i++){
        for (int j = 0; j  < small_size; j++){
            O[large_size*i + j] = 0;
            for (int k = 0; k < small_size; k++){
                O[large_size*i + j] += temp[i*small_size+k]*A1inv[k*large_size+j];
            }
        }
    }

}


__global__ void mat_inv(float*A,float*O,float* temp,int n){

    int ix = blockDim.x * blockIdx.x + threadIdx.x; // global thread ix

    if(ix < n* n){

        // Step 1: copy the element on diagonal
        for (int i = 0; i<n; i++)
            O[ix*(n+1)] = 1/A[ix*(n+1)];
        __syncthreads();
        // Generate n-series
        int k = n;

        int n_series = 0;

        while (k>0){

            k = k /2;

            n_series++;

        }


        for(int i = 0;i<(n_series-1);i++){
            int m_size = 1<<(i);//size_subm[i];
            int x_start = 0;
            int y_start = m_size;
            int jump = 2*m_size;

            // for each matrix
            int j = ix;

            float *A21 = &A[(x_start + j*jump) +  (y_start + j*jump)*n ];
            float *O21 = &O[(x_start + j*jump) +  (y_start + j*jump)*n];
            float *O11 = O21 - m_size*n;
            float *O22 = O21 + m_size;
            //      if(i ==1)
            //              print_mat_d(n,n,ix,A);

            SMM(O11,A21,O22,temp,O21,m_size,n);

            __syncthreads();


        }


    }
}






int main(){

    int n = 8;

    // Use only 1 block

    int num_blocks = 1;

    // Create a matrix

    float* hA = (float*) malloc(n*n*sizeof(float));
    float* hO = (float*) malloc(n*n*sizeof(float));
    float* dA;
    float* dO;
    float* temp;

    cudaMalloc(&dA,n*n*sizeof(float));

    cudaMalloc(&dO,n*n*sizeof(float));

    cudaMalloc(&temp,n*n*sizeof(float));




    for(int i = 0; i<n*n;i++){

        if(i%n>i/n)
            hA[i] = 0;
        else
            hA[i] = 1.1;

    }
    print_mat(n,n,hA);

    cudaMemcpy(dA,hA,n*n*sizeof(float),cudaMemcpyHostToDevice);

    int num_thread_per_block = n;

    mat_inv<<<num_blocks,num_thread_per_block>>>(dA,dO,temp,n);

    cudaMemcpy(hO,dO,n*n*sizeof(float),cudaMemcpyDeviceToHost);
    printf("\n");
    print_mat(n,n,hO);



    return 0;

}

Answer 1

一些初步说明：

如果您对“快速、并行”矩阵求逆感兴趣，您可能应该只使用 library。
我将假设矩阵的边维度是 2 的幂。
我只是想尝试解决您的代码中的一些缺陷，而不是尝试创建“最佳”的东西。我认为与其完善它，不如使用库。
我并不总是 assume 矩阵求逆是必要的。

你勾勒出的一般方法对我来说似乎对矩阵边维 2 的幂来说是合理的：

用倒数逐个元素替换主对角线。
沿着“连续的次对角线方向”工作，使用增加大小的矩阵乘法运算来填充下三角结果。

一般来说，我不清楚您是否真正了解 CUDA 开发的线程并行性质。我在这里看到了一个例子：

    for (int i = 0; i<n; i++)
        O[ix*(n+1)] = 1/A[ix*(n+1)];

让我们研究一下。首先，i 上的循环让我感到困惑，因为 i 没有在循环主体中使用（并且输出与输入不相交）。循环的目的是什么？主体功能是循环不变的。其次，让我们考虑一下您的索引。如果我们信任前面的 if 语句，您的线程的 ix 变量取值范围从 0 到 n*n-1。那是不对的。当然，您在这里只启动 n 个线程，所以它可以工作。但是循环是不必要的。

另请注意，__syncthreads() 在条件代码中是非法的，除非所有线程都可以访问它。因此，只要我们不启动超过 n 个线程（并且我们修改初始 if 语句），我们就可以了。

因为您选择了在每种情况下都由单个线程执行矩阵-矩阵乘法，所以永远不会出现我们需要所有 n 个线程执行矩阵-矩阵乘法的情况，因为我们通过“次对角线”前进（我们需要的最大数字是 n/2，在第一个“次对角线”）。因此，您一定进行了太多的矩阵乘法运算，因为随着您继续通过子对角线来完成工作，您在任何地方都无法使代码进行越来越少的矩阵乘法运算。所以你一定是在做非法索引，如果你运行你的代码带有 cuda-memcheck 或 compute-sanitizer，这是显而易见的，因为我总是建议你在使用 CUDA 代码时遇到问题。换句话说，需要在第二个 i 循环的每次迭代中完成的矩阵乘法不是 n。但这正是您的代码试图做的。

我注意到的另一件事是您没有正确处理 temp。所有线程都使用相同的 temp 数组，有效起始位置为 0。temp 每个线程的偏移量不正确，因此每个线程都使用一个单独的区域。

可能还有其他问题。由于上述问题中的至少 2 个，需要进行重大重写。因此我展示了我自己的代码，并没有试图完全匹配你的代码：

$ cat t12.cu
#include <cstdio>
// matrix-matrix multiply, with optional negation, allowing output to overwrite input
template <typename T>
__host__ __device__
void mm(T *A, T *B, T *D,  int stride, int n, bool negate){

  T *C = new T[n*n];  // avoids temp indexing madness
  for (int i = 0; i < n; i++)
    for (int j = 0; j < n; j++) {
      C[i*n+j] = 0;
      for (int k = 0; k < n; k++)  C[i*n+j] += A[i*stride+k]*B[k*stride+j];}
  for (int i = 0; i < n; i++)
    for (int j = 0; j < n; j++) D[i*stride+j] = (negate)?(-C[i*n+j]):(C[i*n+j]);
  delete[] C;
}

// assumes one block only, n == threads per block, n power of 2
// lower triangular matrix invert
template <typename T>
__global__
void ltmi(T *A, int n){
  int i = threadIdx.x;
    // "invert" main diagonal
  A[i*(n+1)] = 1/A[i*(n+1)];

// process "subdiagonals" in order
  int mask = 1;
  for (int s = 1; s < n; s *= 2){
    __syncthreads();
    if (!(i & mask)){ // select needed threads at each "subdiagonal"
      T *mA = A + i*(n+1); // pointer to thread's A sub matrix
      mm(mA+s*(n+1), mA+s*n, mA+s*n, n, s, true); // -A3i*A2
      mm(mA+s*n, mA, mA+s*n, n, s, false);       // A2*A1i
      mask = 2*mask+1;}
    }
}


typedef float mt;
const int n = 8;  // must be power-of-2, 1024 or less
int main(){
// setup matrix
  mt *A;
  cudaMallocManaged(&A, sizeof(*A)*n*n);
  memset(A, 0, sizeof(*A)*n*n);
  for (int i = 0; i < n; i++)
    for (int j = 0; j < i+1; j++)
      A[i*n+j] = 1.1f;

  ltmi<<<1,n>>>(A, n);
  cudaDeviceSynchronize();

  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++)
      printf("%.6f ",  A[i*n+j]);
    printf("\n");}
}
$ nvcc -o t12 t12.cu -lineinfo
$ compute-sanitizer ./t12
========= COMPUTE-SANITIZER
0.909091 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
-0.909091 0.909091 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 -0.909091 0.909091 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 -0.909091 0.909091 0.000000 0.000000 0.000000 0.000000
-0.000000 -0.000000 0.000000 -0.909091 0.909091 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 -0.909091 0.909091 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 -0.909091 0.909091 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.909091 0.909091
========= ERROR SUMMARY: 0 errors
$

我还没有彻底测试过这个。我真的不推荐这种方法。

通过 cuda 实现快速并行矩阵求逆算法但有些东西不起作用

Implementing a fast parallel Matrix Inversion algorithm by cuda but something is not working

matrix

cuda