LU 使用 openmp 分解

Question

我有一个问题：LU 分解算法的并行版本与序列同时运行：

void lup_od_omp(double* a, int n){
int i,j,k;
for(k = 0; k < n - 1; ++k)
{
    #pragma omp parallel for shared(a,n,k) private(i,j)
    for(i = k + 1; i < n; i++)
    {
        a[i*n + k] /= a[k*n + k];
        for(j = k + 1; j < n; j++)
        {
            a[i*n + j] -= a[i*n + k]*a[k*n + j];
        }
    }
}}

也许我做错了什么？

Answer 1

由于您只在两个内核上工作，因此并行化实际上可能会妨碍矢量化器。SSE2 上的矢量化将为您提供每次操作 2 次双精度的数据带宽，在 AVX 上为 4 次。

双线程有很多同步开销，这可能会减慢您的速度，尤其是在您松散矢量化的情况下。此外，出于某种原因，您的#pragma omp不会启动任何线程，除非调用omp_set_num_threads以实际使其使用线程。

与矢量化有关的另一件事是，并非所有编译器都明白a[i*n + j]旨在处理二维数组，因此最好首先将其声明为这样。

以下是对您的代码的轻微优化，在我的Xeon上运行得相当不错：

void lup_od_omp(int n, double (*a)[n]){
    int i,k;
    for(k = 0; k < n - 1; ++k) {
        // for the vectoriser
        for(i = k + 1; i < n; i++) {
            a[i][k] /= a[k][k];
        }
        #pragma omp parallel for shared(a,n,k) private(i) schedule(static, 64)
        for(i = k + 1; i < n; i++) {
            int j;
            const double aik = a[i][k]; // some compilers will do this automatically
            for(j = k + 1; j < n; j++) {
                a[i][j] -= aik * a[k][j];
            }
        }
    }
}

3000x3000 icc -O2数组的运行时：

Your code sequential:  0:24.61 99%  CPU
Your code 8 threads :  0:05.21 753% CPU
My   code sequential:  0:18.53 99%  CPU
My   code 8 threads :  0:05.42 766% CPU

在另一台机器上，我在 AVX 上对其进行了测试（256 位矢量，每次操作 4 个双精度）：

My code on AVX sequential :  0:09.45 99%  CPU
My code on AVX 8 threads  :  0:03.92 766% CPU

如您所见，我对矢量器进行了一些改进，但对并行部分没有做太多工作。

Answer 2

代码的主要问题是你以不好的方式分解了工作负载。

对于单个 LU 分解，可以并行调用 n-1 次。每次，并行 for 都会做线程分叉和连接，这会带来很多开销。特别是当k很大时，内环（for(i){for(j){...}}）只包含很少的工作。并行它将非常低效。

您可以考虑使用适当的聚合方案来减少开销。有关详细信息，请参阅此幻灯片。

http://courses.engr.illinois.edu/cs554/notes/06_lu_8up.pdf

另一方面，您可以使用现有的性能库来获得 LU 分解的最大性能，例如英特尔 MKL

http://software.intel.com/en-us/node/468682