在 CUDA 中应用 Gauss-Jordan 反演

问题描述

我正在尝试将矩阵求逆应用于给定矩阵，但内核仅适用于最大为 5x5 的矩阵。

如果我使用任何维度更大的矩阵，结果是不正确的。

mod1 = SourceModule("""
__global__ void invert(float* A,float* I,int n) {
    int tx = blockIdx.x * blockDim.x + threadIdx.x;
    int ty = blockIdx.y * blockDim.y + threadIdx.y;
    for(int i = 0; i < n; i++) {
        int col = i * n + ty;
        int row = tx * n + i; 
        if (tx == i) {
            I[tx * n + ty] /= A[i * n + i];
            A[tx * n + ty] /= A[i * n + i];
        }
        if (tx != i) {
            I[tx * n + ty] -= I[col] * A[row];
            A[tx * n + ty] -= A[col] * A[row];
        }
    }
}"""
)

解决方法

求逆矩阵的代码不正确，不能并行计算I的所有元素。我什至惊讶于它适用于高达 5x5 的矩阵。正确的解决方法是将矩阵的行从头到尾串联考虑，将两个矩阵中的行除以A中的对角元素，然后乘以后面的所有行减去它对角线元素下的每一行的元素，你可以并行执行这一步。最后，在对所有行完成此操作后，从最后一行到第一行向后执行相同的操作，这段代码可以向您说明这一点：

void inverse(float* A,float* I,int n)
{
    int i,j;
    size_t size;
    float v;
    float * d_A,* d_I,* d_v;
    size = (unsigned __int64)n * n * sizeof(float);
    cudaMalloc(&d_A,size);
    cudaMemcpy(d_A,A,size,cudaMemcpyHostToDevice);
    cudaMalloc(&d_I,size);
    cudaMalloc(&d_v,sizeof(float));
    for (i = 0; i < n; i++)
        I[i * n + i] = 1;
    for (i = 0; i < n; i++)
    {
        GetVal<<<1,1>>>(d_A,i * (n + 1),d_v); \\ Get value of diagonal element
        cudaMemcpy(&v,d_v,sizeof(float),cudaMemcpyDeviceToHost);
        if (i != n - 1)  \\ Divide row in A matrix starting from element after diagonal
            DivideRow<<<1,n - i - 1>>>(d_A,i * (n + 1) + 1,n - i - 1,v);
        DivideRow<<<1,n>>>(d_I,i * n,n,v);  \\ Divide row in I matrix
        cudaDeviceSynchronize();
        if (i != n - 1)  \\ Subtracting rows
        {
            dim3 GridA(1,1);
            dim3 BlockA(n - i - 1,n - i - 1);
            dim3 GridI(1,1);
            dim3 BlockI(n - i - 1,n);
            ModifyRow<<<GridA,BlockA>>>(d_A,i,i + 1,n - i - 1);
            ModifyRow<<<GridI,BlockI>>>(d_A,d_I,n);
            cudaDeviceSynchronize();
        }
    }
    cudaFree(d_v);
    for (i = n - 1; i > 0; i--)  \\ Backward subtraction
    {
        dim3 GridI(1,1);
        dim3 BlockI(i,n);
        ModifyRow<<<GridI,n);
        cudaDeviceSynchronize();
    }
    cudaMemcpy(I,cudaMemcpyDeviceToHost);
    cudaFree(d_A);
    cudaFree(d_I);
}

__global__ void GetVal(float* A,int p,float* v)
{
    v[0] = A[p];
}

__global__ void DivideRow(float* A,int s,int n,floatd)
{
    int c;
    c = blockIdx.x * blockDim.x + threadIdx.x;
    if (c < n)
        A[s + c] /= d;
}

__global__ void ModifyRow(float* MM,int fr,int fc,float* A,int sr,int sc,int nr,int nc)
{
    int r,c,nA;
    r = blockIdx.x * blockDim.x + threadIdx.x;
    if (r >= nr)
        return;
    c = blockIdx.y * blockDim.y + threadIdx.y;
    if (c >= nc)
        return;
    nA = sc + nc;
    A[(sr + r) * nA + sc + c] -= MM[(sr + r) * n + fc] * A[fr * nA + sc + c];
}

请注意最大块大小为 1024，因此如果您的矩阵大于 32x32，则必须修改网格和块大小。