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

Question

我试着理解优化例程。我的重点是代码中最关键的部分(代码有一些长度为"nc“的循环和一个长度"np”的循环，其中数字"np“比”nc“大得多)。我在这里展示了部分代码。其余的代码在计算时间的百分比中并不是很重要，所以在算法的其余部分中，我更喜欢代码purify。然而，具有"np“长度的临界循环是一段相当简单的代码，可以并行化。所以，如果我把这个部分改写成一些更有效和不那么清晰的版本(可能是SSE指令)，这不会有什么影响。我使用gcc编译器、c++代码和OpenMP并行化。

这段代码是众所周知的细胞内粒子算法的一部分(这个算法也是基本的)。我正在尝试学习这个版本上的代码优化(所以我的目标不仅仅是拥有有效的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还原：

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);

}