▶ 使用Jacobi 迭代求泊松方程的数值解
● 原始串行版本
#include <stdio.h>
#include <stdlib.h>
#include <math.h> #if defined(_WIN32) || defined(_WIN64) // 统一计时器
#include <C:\Program Files\PGI\win64\19.4\include\wrap\sys\timeb.h>
#define gettime(a) _ftime(a)
#define usec(t1,t2) ((((t2).time - (t1).time) * 1000 + (t2).millitm - (t1).millitm)) // 单位 ms
typedef struct _timeb timestruct;
#else
#include <sys/time.h>
#define gettime(a) gettimeofday(a, NULL)
#define usec(t1,t2) (((t2).tv_sec - (t1).tv_sec) * 1000000 + (t2).tv_usec - (t1).tv_usec) // 单位 us
typedef struct timeval timestruct;
#endif #define IMPROV // 是否额外使用 “每次计算的修正量” 作为退出循环的条件 inline float uval(float x, float y) // 求该点到原点距离的平方
{
return x * x + y * y;
} int main()
{
const int row = , col = ; // 网格行数和列数，
const float height = 1.0, width = 2.0; // 实际高度和宽度，与网格行列数不成比例说明是矩形网格
const float hx = height / row, wy = width / col; // 每个网格的高度和宽度
const float fij = -4.0f; // 函数 f(x,y) = -4，此时方程的解为 z = x^2 + y^2
const float hx2 = hx * hx, wy2 = wy * wy, c1 = hx2 * wy2, c2 = 1.0f / (2.0 * (hx2 + wy2));// 其他用到的参数
const int maxIter = ; // 最大迭代次数
const int colPlus = col + ; // 实际列数
#ifdef IMPROV
const float errControl = 0.0f; // 修正量控制，取 0 表示无用
float err = 0.0f; // 修正量
#endif float *u0 = (float *)malloc(sizeof(float)*(row + )*colPlus); // 用来存放网格数据的两张表，行列数等于 row 和 col 各自加 1，
float *u1 = (float *)malloc(sizeof(float)*(row + )*colPlus);
float *utemp = NULL; // 用于交换 u1 和 u0 的临时指针 // 初始化边界为 g(x,y) = x^2+y^2
for (int ix = ; ix <= row; ix++) // 左右边界
{
u0[ix*colPlus + ] = u1[ix*colPlus + ] = uval(ix * hx, 0.0f);
u0[ix*colPlus + col] = u1[ix*colPlus + col] = uval(ix*hx, col * wy);
}
for (int jy = ; jy <= col; jy++) // 上下边界
{
u0[jy] = u1[jy] = uval(0.0f, jy * wy);
u0[row*colPlus + jy] = u1[row*colPlus + jy] = uval(row*hx, jy * wy);
}
for (int ix = ; ix < row; ix++) // 内部格点初始化为 0.0f
{
for (int jy = ; jy < col; jy++)
u0[ix*colPlus + jy] = 0.0f;
} // 计算
timestruct t1, t2;
gettime(&t1);
for (int iter = ; iter < maxIter; iter++)
{
for (int ix = ; ix < row; ix++)
{
for (int jy = ; jy < col; jy++)
{
u1[ix*colPlus + jy] = (c1*fij + wy2 * (u0[(ix - )*colPlus + jy] + u0[(ix + )*colPlus + jy]) + \
hx2 * (u0[ix*colPlus + jy - ] + u0[ix*colPlus + jy + ])) * c2;
#ifdef IMPROV
err = max(fabs(u0[ix*colPlus + jy] - u1[ix*colPlus + jy]), err); // 记录整张表上的最大修正量
#endif
}
}
#ifdef IMPROV
//printf("\niter = %d, err = %e\n", iter, err); // 逐次输出
if (err < errControl) // 修正量小于指定量就可以退出
break;
#endif
utemp = u0, u0 = u1, u1 = utemp; // 交换指针
}
gettime(&t2); long long timeElapse = usec(t1, t2);
printf("\nElapsed time: %13ld ms.\n", timeElapse);
free(u0);
free(u1);
getchar();
return ;
}
● 输出结果（使用 IMPROV），可以看到很多 not fused，这都是可以改进的地方
D:\Code\OpenACC>pgcc main.c -Minfo -o main.exe // 普通编译
main:
, FMA (fused multiply-add) instruction(s) generated // 使用乘加指令
uval:
, FMA (fused multiply-add) instruction(s) generated D:\Code\OpenACC>pgcc main.c -Minfo -o main-fast.exe -fast // 添加 fast 选项
main:
, uval inlined, size= (inline) file main.c () // 4 个内联函数
, Loop not fused: different loop trip count // 担心 for 中存在数据依赖，拒绝并行
Generated vector and scalar versions of the loop; pointer conflict tests determine which is executed // 担心 u0 和 u1是否重叠，拒绝并行
Loop not vectorized: data dependency
Loop unrolled times // 循环展开
Generated prefetches in scalar loop
, uval inlined, size= (inline) file main.c ()
, uval inlined, size= (inline) file main.c ()
, Loop not vectorized: data dependency
Loop unrolled times
, uval inlined, size= (inline) file main.c ()
, Memory zero idiom, loop replaced by call to __c_mzero4 // 使用 memcpy 来赋零值
, Loop not vectorized/parallelized: potential early exits // 有额外脱离循环的条件，拒绝并行
, Loop not vectorized: data dependency
Loop unrolled times
FMA (fused multiply-add) instruction(s) generated D:\Code\OpenACC>main.exe Elapsed time: ms. D:\Code\OpenACC>main-fast.exe Elapsed time: ms. // 加了 fast 反而更慢
● 输出结果（不用 IMPROV），发现变快了，可见提前跳出循环的 if 语句对并行化有很大影响。在本例中我们让 errControl = 0，每次循环多一个判断（实际绝对不会跳出），就严重干扰了编译
D:\Code\OpenACC>pgcc main.c -Minfo -o main.exe
main:
, FMA (fused multiply-add) instruction(s) generated
uval:
, FMA (fused multiply-add) instruction(s) generated D:\Code\OpenACC>pgcc main.c -Minfo -o main-fast.exe -fast
main:
, uval inlined, size= (inline) file main.c ()
, Loop not fused: different loop trip count
Generated vector and scalar versions of the loop; pointer conflict tests determine which is executed
Loop not vectorized: data dependency
Loop unrolled times
Generated prefetches in scalar loop
, uval inlined, size= (inline) file main.c ()
, uval inlined, size= (inline) file main.c ()
, Loop not vectorized: data dependency
Loop unrolled times
, uval inlined, size= (inline) file main.c ()
, Memory zero idiom, loop replaced by call to __c_mzero4
, Loop not fused : function call before adjacent loop // ？
, Loop not vectorized : data dependency
Loop unrolled times // 展开次数由 2 变成 4
FMA (fused multiply-add) instruction(s) generated D:\Code\OpenACC>main.exe Elapsed time: ms. // 变快了 1 倍 D:\Code\OpenACC>main-fast.exe Elapsed time: ms. // 再变快 1 倍
● 使用 OpenMP 优化（就一句导语）
// #include <math.h> 下面
#include <omp.h> //for (int iter = 1; iter < maxIter; iter++){ 下面
#ifdef IMPROV
#pragma omp parallel for reduction(max:err) default(none) shared(u0, u1, c1, c2, hx2, wy2, colPlus) private(err)
#else
#pragma omp parallel for default(none) shared(u0, u1, c1, c2, hx2, wy2, colPlus)
#endif
● 输出结果
D:\Code\OpenACC>set OMP_NUM_THREADS= // 使用 4 个线程 D:\Code\OpenACC>pgcc main.c -Minfo -o main4I.exe -fast -mp // 用 IMPROV
main:
, uval inlined, size= (inline) file main.c ()
, Loop not fused: different loop trip count
Generated vector and scalar versions of the loop; pointer conflict tests determine which is executed
Loop not vectorized: data dependency
Loop unrolled times
Generated prefetches in scalar loop
, uval inlined, size= (inline) file main.c ()
, uval inlined, size= (inline) file main.c ()
, Loop not vectorized: data dependency
Loop unrolled times
, uval inlined, size= (inline) file main.c ()
, Memory zero idiom, loop replaced by call to __c_mzero4
, Loop not vectorized/parallelized: potential early exits
, Parallel region activated // OpenMP 并行区
Parallel loop activated with static block schedule
, Loop not vectorized: data dependency
Loop unrolled times
FMA (fused multiply-add) instruction(s) generated
, Begin critical section // 脱出循环的判断导致的串行区
End critical section
Barrier // 栅栏
Parallel region terminated D:\Code\OpenACC>main4I.exe Elapsed time: ms. // 还是快了 3.8 倍 D:\Code\OpenACC>pgcc main.c -Minfo -o main4.exe -fast -mp // 不用 IMPROV
main:
, uval inlined, size= (inline) file main.c ()
, Loop not fused: different loop trip count
Generated vector and scalar versions of the loop; pointer conflict tests determine which is executed
Loop not vectorized: data dependency
Loop unrolled times
Generated prefetches in scalar loop
, uval inlined, size= (inline) file main.c ()
, uval inlined, size= (inline) file main.c ()
, Loop not vectorized: data dependency
Loop unrolled times
, uval inlined, size= (inline) file main.c ()
, Memory zero idiom, loop replaced by call to __c_mzero4
, Loop not vectorized/parallelized: contains a parallel region // 有 OpenMP的并行区，拒绝并行
, Parallel region activated
Parallel loop activated with static block schedule
, Loop not vectorized: data dependency
Loop unrolled times
FMA (fused multiply-add) instruction(s) generated
, Barrier // 没有了串行区
Parallel region terminated D:\Code\OpenACC>main4.exe Elapsed time: ms. // 还能再快点，加速比 1.4 D:\Code\OpenACC>set OMP_NUM_THREADS= // 使用 8 线程 D:\Code\OpenACC>pgcc main.c -Minfo -o main8.exe -fast -mp ...// 跟 4 线程时一模一样 D:\Code\OpenACC>main8.exe Elapsed time: ms. // 不宰线性加速，加速比 1.5
▶ 在 Ubuntu 下跑的结果，加速前比 win10 慢很多，关闭 IMPROV 并开启 OpenMP 和 fast 选项后速度接近
mainI.exe us mainI-fast.exe us // 极速比 3.1 main.exe us // 加速比 2.1 main-fast.exe us // 加速比 6.4 cuan@CUAN:~$ pgcc mainI.c -Minfo -o main4I-fast.exe -fast -mp // 要求我将 row，col，fij 放入 OpenMP 的 shared 导语中，在 win10 下没有显式放入也行
PGC-S--row must appear in a proper data sharing clause (e.g., PRIVATE) (mainI.c: )
PGC-S--col must appear in a proper data sharing clause (e.g., PRIVATE) (mainI.c: )
PGC-S--fij must appear in a proper data sharing clause (e.g., PRIVATE) (mainI.c: )
PGC/x86- Linux 19.4-: compilation completed with severe errors main4I-fast.exe us main4-fast.exe us // 加速比 8.8 main8-fast.exe us // 不能继续线性加速