并行计算的 C++ 代码优化示例

Question

我正在尝试了解优化例程。我专注于我的代码中最关键的部分（代码有一些长度为“nc”的循环和一个长度为“np”的循环，其中数字“np”比“nc”大得多）。我在这里展示了部分代码。其余代码在计算时间的百分比上不是很重要，因此我更喜欢在算法的其余部分中进行代码纯化。但是，具有“np”长度的关键循环是一段非常简单的代码，并且可以并行化。因此，如果我将这部分重写为一些更有效且不太清晰的版本（可能是 SSE 指令），这不会有什么坏处。我正在使用 gcc 编译器、c++ 代码和 OpenMP 并行化。

这段代码是众所周知的particle-in-cell 算法的一部分（这也是基本的）。我正在尝试在这个版本上学习代码优化（所以我的目标不仅仅是拥有有效的 PIC 算法，因为它已经编写了数千种变体，但我也想为代码优化带来一些示范性示例）。我正在尝试做一些工作，但我不太确定我是否正确解决了所有优化属性。

const int NT = ...; // number of threads (in two versions: about 6 or about 30)
const int np = 10000000; // np is about 1000-10000 times larger than nc commonly
const int nc = 10000;
const int step = 1000;

float u[np], x[np];
float a[nc], a_lin[nc], rho_full[NT][nc], rho_diff[NT][nc] , weight[nc];

int p,num;

for ( i = 0 ; i<step ;  i++) {

// *** 
// *** some not very time consuming code for calculation 
// *** a, a_lin from values of rho_full and rho_diff

#pragma omp for private(p,num)
for ( k = np ; --k ; ) {

    num = omp_get_thread_num();

    p = (int) x[k];
    u[k] += a[p] + a_lin[p] * (x[k] - p);
    x[k] += u[k];

    if (x[k]<0 ) {x[k]+=nc;} else
    if (x[k]>nc) {x[k]-=nc;};

    p = (int) x[k];
    rho_full[num][p] += weight[k];
    rho_diff[num][p] += weight[k] * (x[k] - p);

    }
};

我意识到这有问题：

1)（主要问题）我使用一组数组rho_full[num][p]，其中num 是每个线程的索引。计算后我只是总结了这个数组（rho_full[0][p] + rho_full[1][p] + rho_full[2][p] ...）。原因是避免使用两个不同的线程写入数组的同一部分。我不太确定这种方式是否有效（注意数字“nc”相对较小，因此“np”的操作数量可能仍然是最重要的）

2)（也是一个重要问题）我需要多次阅读x[k]，它也被多次更改。也许最好将此值读入某个寄存器，然后忘记整个 x 数组或在此处修复一些指针。在所有计算之后，我可以再次调用x[k] 数组并存储获得的值。我相信编译器会为我完成这项工作，但我不太确定，因为我在算法中心使用了 x[k] 的修改。所以编译器可能会自己做一些有效的工作，但也许在这个版本中，它调用的次数比我调用和存储这个值的次数要多。

3)（可能不相关）代码使用整数部分和小数点部分以下的余数。它需要这两个值。我将整数部分识别为p = (int) x，将余数识别为x - p。我在循环内部开始和结束时计算这个例程。可以看到这种拆分可以存储在某处并在下一步使用（我的意思是i 索引处的步骤）。你觉得以下版本更好吗？我将整数和余数部分存储在 x 的数组中，而不是整个值 x。

int x_int[np];
float x_rem[np];
//...

for ( k = np ; --k ; ) {

    num = omp_get_thread_num();

    u[k] += a[x_int[k]] + a_lin[x_int[k]] * x_rem[k];
    x_rem[k] += u[k];

    p = (int) x_rem[k];   // *** This part is added into code for simplify the rest.
    x_int[k] += p;        // *** And maybe there is a better way how to realize
    x_rem[k] -= p;        // *** this "pushing correction".

    if (x_int[k]<0 ) {x_int[k]+=nc;} else
    if (x_int[k]>nc) {x_int[k]-=nc;};

    rho_full[num][x_int[k]] += weight[k];
    rho_diff[num][x_int[k]] += weight[k] * x_rem[k];

    }
};

Answer 1

您可以在 for 循环中使用 OMP reduction：

int result = 0;

#pragma omp for nowait reduction(+:result) 
for ( k = np ; --k ; ) {

    num = omp_get_thread_num();

    p = (int) x[k];
    u[k] += a[p] + a_lin[p] * (x[k] - p);
    x[k] += u[k];

    if (x[k]<0 ) {x[k]+=nc;} else
    if (x[k]>nc) {x[k]-=nc;};

    p = (int) x[k];
    result += weight[k] + weight[k] * (x[k] - p);

}