高斯消除并行答案

【问题标题】：Gaussian Elimination Parallelism高斯消除并行
【发布时间】：2017-12-16 15:40:30
【问题描述】：

我已经成功地在 CUDA 中实现了一个单线程程序以进行高斯消除，并希望实现并行性。到目前为止，并行代码如下所示：

__global__ void ParallelGaussian(double* A)
{
    int index = threadIdx.x;
    int stride = blockDim.x;

    if (index < ROWS) //Skip additional threads
    {
        for (unsigned int r = index; r < ROWS; r += stride)
        {
            //Forward elimination to reduce to row echelon form
            for (unsigned int k = r + 1; k < ROWS; ++k)
            {
                double c = -A[(ROWS + 1) * k + r] / A[(ROWS + 1) * r + r];
                for (unsigned int j = r; j < ROWS + 1; ++j)
                {
                    if (r == j)
                        A[(ROWS + 1) * k + j] = 0.0;
                    else
                        A[(ROWS + 1) * k + j] += c * A[(ROWS + 1) * r + j];
                }
            }
        }
    }
}

正如我们所见，GPU 上的代码会将一维数组（矩阵）转换为下三角矩阵，然后在 CPU 上我将继续进行反向替换以获得最终结果。这种方法没有进行旋转，因为它不是完全需要的，但确实提高了算法的数值稳定性。

使用单线程和块启动内核可以将矩阵转换为行梯形：

ParallelGaussian << < 1, 1 >> >(dev_a);

但是，如果我想增加线程数，比如

ParallelGaussian << < 1, 32 >> >(dev_a);

它将无法生成下三角矩阵。现在将 __syncthreads() 调用添加到代码中以同步块中的线程并不能改善这种情况，我不知道为什么。

【问题讨论】：

标签： math parallel-processing cuda gpu

【解决方案1】：

考虑您的内部循环。每个线程都访问A，并且由于k 和j 从r 运行到矩阵的末尾，因此多个线程有可能修改相同的A[(ROWS + 1) * k + j] 值。

您还可能有一些线程访问A[(ROWS + 1) * r + j]，而其他线程正在更新该值。

一种可能的解决方案是让每个线程累积到单独的结果数组中，然后在最后将它们组合起来。这是内存密集型的。

另一种方法是对其进行重组，以便只有一个线程将写入特定值，并将这些值存储在一个新矩阵中（这样您就不会更改其他线程可能需要的任何值）。

【讨论】：