CUDA高斯-约旦矩阵求逆GPU比CPU慢10倍，求优化建议

我在CPU上实现了稳定的简易高斯-约旦矩阵求逆算法，移植到GPU后功能正常，但速度仅为CPU的1/10。测试数据如下：

matrix size = 88; CPU time = 6531 mcs; GPU time = 52806 mcs;

CPU实现代码 #

void GaussJordanInverse(float* original, float* temp, float* singular, float* bigmatrix, int size) {
    //create temp
    //create singular matrix
    for (int i = 0; i < size; i++)
    {
        for (int j = 0; j < size; j++) {
            temp[i * size + j] = original[i * size + j];
            singular[i * size + j] = 0; 
        }
        singular[i * size + i] = 1;
    }
    //create big matrix
    for (int i = 0; i < size; i++) {
        for (int j = 0; j < size; j++)
        {
            bigmatrix[i * (size * 2) + j] = temp[i * size + j];
            bigmatrix[i * (size * 2) + j + size] = singular[i * size + j];
        }
    }
    //direct
    for (int k = 0; k < size; k++)
    {
        for (int i = 0; i < 2 * size; i++) bigmatrix[k * (size * 2) + i] = bigmatrix[k * (size * 2) + i] / temp[k * size + k];
        for (int i = k + 1; i < size; i++) {
            float K = bigmatrix[i * (size * 2) + k] / bigmatrix[k * (size * 2) + k];
            for (int j = 0; j < 2 * size; j++)bigmatrix[i * (size * 2) + j] = bigmatrix[i * (size * 2) + j] - bigmatrix[k * (size * 2) + j] * K;
        }
        for (int i = 0; i < size; i++) {
            for (int j = 0; j < size; j++)
            {
                temp[i * size + j] = bigmatrix[i * (size * 2) + j];
            }
        }
    }
    //indirect
    for (int k = size - 1; k > -1; k--)
    {
        for (int i = 2 * size - 1; i > -1; i--) bigmatrix[k * (size * 2) + i] = bigmatrix[k * (size * 2) + i] / temp[k * size + k];
        for (int i = k - 1; i > -1; i--)
        {
            float K = bigmatrix[i * (size * 2) + k] / bigmatrix[k * (size * 2) + k];
            for (int j = 2 * size - 1; j > -1; j--)
                bigmatrix[i * (size * 2) + j] = bigmatrix[i * (size * 2) + j] - bigmatrix[k * (size * 2) + j] * K;
        }
    }
    //cut
    for (int i = 0; i < size; i++) {
        for (int j = 0; j < size; j++)
        {
            singular[i * size + j] = bigmatrix[i * (size * 2) + j + size];
            original[i * size + j] = singular[i * size + j];
        }
    }
}

GPU实现代码 #

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <cuda.h>
#include <iostream>
#ifndef __CUDACC__  
#define __CUDACC__
#endif
#include <device_functions.h>
#include <stdio.h>
#include <chrono>
#include "matrixconversion.h"
using namespace std;
using namespace std::chrono;

__global__ void create(float* original, float* temp, float* singular, float* bigmatrix, int size) {
    int x = blockDim.x * blockIdx.x + threadIdx.x;
    int y = blockDim.y * blockIdx.y + threadIdx.y;
    //create singular matrix
    if (x < size && y < size)
    {
        temp[x * size + y] = original[x * size + y];
        if (x == y) singular[x * size + y] = 1;
        //singular[x * size + y] = 0;
        //create big matrix
        bigmatrix[x * (2 * size) + y] = temp[x * size + y];
        bigmatrix[x * (2 * size) + y + size] = singular[x * size + y];
    }
}
__global__ void direct_kernel(float* temp, float* bigmatrix, int size, int k) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    if (x < 2 * size) {
        bigmatrix[k * (size * 2) + x] /= temp[k * size + k];
    }
    __syncthreads();
    if (x >= k + 1 && x < size * 2) {
        float K = bigmatrix[x * (size * 2) + k] / bigmatrix[k * (size * 2) + k];
        for (int j = 0; j < 2 * size; j++) {
            bigmatrix[x * (size * 2) + j] -= bigmatrix[k * (size * 2) + j] * K;
        }
    }
    __syncthreads();
    if (x < size) {
        for (int j = 0; j < size; j++) {
            temp[x * size + j] = bigmatrix[x * (size * 2) + j];
        }
    }
}
__global__ void indirect_kernel(float* temp, float* bigmatrix, int size, int k) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;

    if (x < 2 * size) {
        bigmatrix[k * (size * 2) + x] /= temp[k * size + k];
    }
    __syncthreads();
    if (x < k&& x >= 0) {
        float K = bigmatrix[x * (size * 2) + k] / bigmatrix[k * (size * 2) + k];
        for (int j = 2 * size - 1; j > -1; j--) {
            bigmatrix[x * (size * 2) + j] -= bigmatrix[k * (size * 2) + j] * K;
        }
    }
    __syncthreads();
}

__global__ void cut(float* original, float* singular, float* bigmatrix, int size) {
    int x = blockDim.x * blockIdx.x + threadIdx.x;
    int y = blockDim.y * blockIdx.y + threadIdx.y;

    if (x < size && y < size) {
        singular[x * size + y] = bigmatrix[x * (size * 2) + y + size];
        original[x * size + y] = singular[x * size + y];
    }
}

void GPUInverse(float* copy, float* original, int size, long int* time) {
    Copy1DFloat(copy, original, size);
    //set kernels
    int tr = 16;
    int bl = (int)ceil(size / tr);
    dim3 grid_size(16, 16);
    dim3 block_size(bl, bl);
    //create temp
    float* d_copy, * d_temp, * d_singular, * d_bigmatrix;
    //memory
    cudaMalloc((void**)&d_copy, size * size * sizeof(float));
    cudaMalloc((void**)&d_temp, size * size * sizeof(float));
    cudaMalloc((void**)&d_singular, size * size * sizeof(float));
    cudaMalloc((void**)&d_bigmatrix, size * (size * 2) * sizeof(float));
    cudaMemset(d_copy, 0, size * size * sizeof(float));
    cudaMemset(d_temp, 0, size * size * sizeof(float));
    cudaMemset(d_singular, 0, size * size * sizeof(float));
    cudaMemset(d_bigmatrix, 0, size * (2 * size) * sizeof(float));
    //to device
    cudaMemcpy(d_copy, copy, size * size * sizeof(float), cudaMemcpyHostToDevice);

    auto start = high_resolution_clock::now();
    //fill temp
    create <<<grid_size, block_size >>> (d_copy, d_temp, d_singular, d_bigmatrix, size);
    cudaDeviceSynchronize();
    //direct
    for (int k = 0; k < size; k++) {
        direct_kernel <<<(2 * size + 255) / 256, 256 >>> (d_temp, d_bigmatrix, size, k);
        cudaDeviceSynchronize();
    }
    //indirect
    for (int k = size - 1; k > -1; k--) {
        indirect_kernel <<<(2 * size + 255) / 256, 256 >>> (d_temp, d_bigmatrix, size, k);
        cudaDeviceSynchronize();
    }
    //cut
    cut <<<grid_size, block_size >>> (d_copy, d_singular, d_bigmatrix, size);
    cudaDeviceSynchronize();

    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(stop - start);
    *time = duration.count();
    //from device
    cudaMemcpy(copy, d_copy, size * size * sizeof(float), cudaMemcpyDeviceToHost);
    //clean

    cudaFree(d_copy);
    cudaFree(d_temp);
    cudaFree(d_singular);
    cudaFree(d_bigmatrix);
}

本人对C++和CUDA了解不深，希望得到优化该GPU实现的建议。

优化建议 #

1. 消除不必要的同步开销 #

• 目前每个核函数调用后都使用cudaDeviceSynchronize()，强制CPU等待GPU完成任务，完全抵消了GPU异步执行的优势。仅在需要读取GPU计算结果（如计时结束后、内存拷贝前）保留一次同步即可。
• 核函数内部的__syncthreads()并非必要：比如direct_kernel中三个阶段的线程工作范围互不重叠，不需要全局同步，移除这些指令可减少线程等待时间。

2. 重构核函数，充分利用GPU并行性 #

当前direct_kernel和indirect_kernel中，单个线程通过for循环处理一整行的所有元素，完全没利用GPU多线程并行能力。应改为每个线程负责一个元素的计算：

• 将主行归一化拆分为单独核函数，每个线程处理主行的一个元素；
• 消去其他行的操作，让每个线程负责某一行的某一个元素，避免单线程串行循环。

3. 优化线程映射与资源利用率 #

• create和cut核函数的grid_size设置为固定16x16，对于88x88的矩阵会产生大量空闲线程。应根据矩阵大小计算合适的网格尺寸：

  dim3 block_size(16, 16); // 常用256线程块（16*16）
  dim3 grid_size((size + block_size.x - 1)/block_size.x, (size + block_size.y - 1)/block_size.y);
  

• direct_kernel和indirect_kernel的网格大小计算应更准确，避免创建过多冗余线程块。

4. 减少冗余内存操作 #

• CPU实现中temp矩阵用于存储中间结果，但GPU中每次direct_kernel都复制整个bigmatrix左半部分到temp，完全不必要。可直接使用bigmatrix左半部分替代temp，节省内存带宽和拷贝时间。
• 移除不必要的cudaMemset：d_copy会被后续cudaMemcpy覆盖，无需提前清零；d_temp可在create核函数中直接赋值，无需提前清零。

5. 拆分多功能核函数 #

当前direct_kernel包含主行归一化、消去其他行、复制temp三个独立步骤，导致不同阶段大量线程空闲。将这三个步骤拆分为三个独立核函数，让每个核函数的线程专注于单一任务，提升线程利用率。

6. 优化全局内存访问模式 #

GPU全局内存需要合并访问才能达到最高带宽。当前bigmatrix采用行优先存储，处理行数据时是合并访问，可进一步确保线程访问的内存地址连续，例如让线程索引对应元素的全局线性索引而非行号。

内容的提问来源于stack exchange，提问作者Gagy Krayper