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

Posted 2023-03-23

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

【英文标题】：Implementing a fast parallel Matrix Inversion algorithm by cuda but something is not working

【发布时间】：2022-01-19 02:52:21

【问题描述】：

目前，我正在实施

用于 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;

【问题讨论】：

我看不出您编写的代码在数学上与您引用的图像相同。我在该代码的任何地方都看不到任何子矩阵反转。你认为O[ix*(n+1)] = 1/A[ix*(n+1)]; 在做什么？

这是反转 1x1 矩阵以开始

但这在数学上是不可能的。如果这些块矩阵中的任何一个是 1 x 1，那么根据定义，其他块矩阵对于分解为 1x1 块的 2x2 矩阵以外的任何东西都不是正方形且不可逆转。

【参考方案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 线程执行矩阵-矩阵乘法的情况，因为我们通过“subdiagonal”（我们需要的最大数字是n/2，在第一个“subdiagonal”处）。因此，您必须进行过多的矩阵乘法运算，因为当您通过次对角线进行工作时，您在任何地方都没有使代码执行越来越少的矩阵乘法运算。因此，您必须进行非法索引，如果您使用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 的新手...感谢您的评论！

你的评论和代码对我学习cuda肯定很有帮助。真的让我大开眼界。再次感谢！