目录
利用Jacobi迭代对二维泊松方程进行求解。采用OpenACC并行编程的方式加速迭代过程。
主要思想：用二阶中心差分来近似二阶偏导数，将差分近似代入泊松方程，得到五点差分离散格式，则将泊松方程的求解转换为稀疏线性方程组的求解。然后就可以利用经典的Jacobi算法来求解此方程组。从任意一个初始近似解开始出发，进行迭代计算，直至最终的近似解满足误差要求。
1. 搭建OpenACC运行环境
1.1 安装CUDA11.1
1.2 配置CUDA环境
配置.bashrc文件，添加环境变量，在.bashrc文件的最后添加如下路径。
配置路径如下：（注意，确保cuda的安装地址正确）
export PATH=/usr/local/cuda/bin:$PATH
export LD_LIBRARY_PATH=/usr/local/cuda/lib64:$LD_LIBRARY_PATH1.3 安装OpenACC 21.3
此处需要注意最好是安装的OpenACC应支持前一个版本的CUDA即，Bundled with the newest plus two previous CUDA versions
安装命令如下
$ wget https://developer.download.nvidia.com/hpc-sdk/21.3/nvhpc-21-3_21.3_amd64.deb \
https://developer.download.nvidia.com/hpc-sdk/21.3/nvhpc-2021_21.3_amd64.deb \
https://developer.download.nvidia.com/hpc-sdk/21.3/nvhpc-21-3-cuda-multi_21.3_amd64.deb
$ apt-get install ./nvhpc-21-3_21.3_amd64.deb ./nvhpc-2021_21.3_amd64.deb ./nvhpc-21-3-cuda-multi_21.3_amd64.deb1.4 配置openACC环境
同样是在.bashrc文件中，在文件最后添加路径，配置如下（注意，此处的OpenACC安装路径是自己实际的安装路径）：
export PATH=/opt/nvidia/hpc_sdk/Linux_x86_64/2021/compilers/bin:$PATH
export PGI_ACC_TIME=1
export PGI_ACC_NOTIFY=1执行pgcc --version，出现openacc版本内容，则说明安装成功，环境配置成功。2. 程序执行
2.1 GPU上并行，用kernels构件加速二重循环
// Jacobi_GPU_noerror.c
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#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))
typedef struct timeval timestruct;
#define Mx 8191 //定义网格尺寸Mx和Ny，此处为了数组的内存对齐，所以设为8192-1
#define Ny 1023 //1024-1
float uval(float x, float y){return (x*x+y*y);}
int main()
{
float Widthy = 2.0, Heightx = 1.0; //Widthy为Ω=(0,W)x(0,H)中的W和H
//添加restrict关键字，使得u0和u1指向的地址空间不交叉
float *restrict u0, *restrict u1; //定义两个一维数组u0和u1，来保存二维格点上的旧值和新值
float hx, hy, hx2, hy2, fij, c1, c2;
float *tprt;
int ix, jy, maxIter, iter;
int Nyp1 = Ny+1;
timestruct t1, t2;
long long telapsed;
u0 = (float*)malloc(sizeof(float)*(Mx+1)*(Ny+1));
u1 = (float*)malloc(sizeof(float)*(Mx+1)*(Ny+1));
maxIter = 100; //指定最大迭代次数
hx = Heightx/Mx;
hy = Widthy/Ny;
//初始化u0/u1的左右边界，这些边界再计算中保持不变
for(ix = 0; ix <= Mx; ix++)
{
u0[ix*Nyp1+0] = u1[ix*Nyp1+0] = uval(ix*hx, 0.0f);
u0[ix*Nyp1+Ny] = u1[ix*Nyp1+Ny] = uval(ix*hx, Ny*Ny);
}
// 初始化u0/u1的上下边界，固定值保持不变
for(jy = 0; jy <= Ny; jy++)
{
u0[jy]=u1[jy]=uval(0.0f, jy*hy);
u0[Mx*Nyp1+jy] = u1[Mx*Nyp1+jy] = uval(Mx*hx, jy*hy);
}
//为数组u0内部元素赋初值，这里指定为0
for(ix = 1; ix < Mx; ix ++)
for(jy = 1; jy < Ny; jy++)
{
u0[ix*Nyp1+jy]=0.0f;
}
fij = -4.0f;
c1 = hx*hx * hy*hy;
c2 = 1.0f/(2.0*(hx*hx+hy*hy));
hx2 = hx*hx;
hy2 = hy*hy;
gettime(&t1);
//开始迭代，每一轮迭代中用u0更新u1同时计算两者间的最大误差值uerr，然后将u1中的新值赋给u0，准备进行下一步迭代
for(iter = 1; iter <= maxIter; iter++)
{
#pragma acc kernels copy(u0[0:(Mx+1)*(Ny+1)],u1[0:(Mx+1)*(Ny+1)]) //应为u1,u0是动态数组，所以使用copy子语。
//copy的功能是数据传递，进入计算构件时，将数组u0,u1从主机内存复制到设备内存上，为并行计算做准备，计算完成后在将数组复制回主机内存。
{
#pragma acc loop independent
for(ix = 1; ix < Mx; ix++)
for(jy = 1; jy < Ny; jy++)
{
u1[ix*Nyp1+jy]=(c1*fij+hy2*(u0[(ix-1)*Nyp1+jy] + u0[(ix+1)*Nyp1+jy]) \
+ hx2*(u0[ix*Nyp1+jy-1] + u0[ix*Nyp1+jy+1]))*c2;
}
}
tprt = u0; u0 = u1; u1 = tprt;
}
gettime(&t2);
telapsed = usec(t1, t2);
printf("ElapsedTime = %13lld microseconds\n", telapsed); //打印循环主体消耗的时间
free(u0);
free(u1);
return 0;
}
执行代码：
pgcc -acc -Minfo Jacobi_GPU_noerror.c -o Jacobi_GPU_noerror.exe
./Jacobi_GPU_noerror.exe执行结果如下图
给出的结果表明编译器识别到了OpenACC导语100次，第62行和63行的for循环被成功并行化。
ElapsedTime = 6863439 microseconds //循环的墙上时间，从开始到自然结束的自然流逝的时间
Accelerator Kernel Timing data
/home/course_user1/openAcc/Jacobi_GPU_noerror.c
main NVIDIA devicenum=0
time(us): 4,127,249 //运行所耗总时间
62: compute region reached 100 times //进入计算区域100次
63: kernel launched 100 times //启动内核100次
grid: [65393] block: [128] //线程组织形式，内核使用的线程网格(grid)包含65393个线程块(block)，每个线程块包含128个线程。
elapsed time(us): total=893,758 max=9,052 min=8,832 avg=8,937 //内核墙上时间，内核运行总时间为893,758微秒，最大9052微秒，最小8832微秒，平均8937微秒。
62: data region reached 200 times //进入数据区域200次
62: data copyin transfers: 400 //向设备复制数据(H2D)次数200次
device time(us): total=2,125,445 max=5,348 min=5,293 avg=5,313 //详细时间，总的数据复制到设备的时间为2125445微秒，最大5348微秒，最小5293微秒，平均5313微秒。
67: data copyout transfers: 600 //向主机复制数据(D2H)次数600次
device time(us): total=2,001,804 max=5,012 min=5 avg=3,336 //详细时间，总的数据复制到主机的时间为2001804微秒，最大5012微秒，最小5微秒，平均3336微秒。
第62行是主体迭代的入口，内含100次迭代，而每次迭代都要将u0和u1分别复制到设备内存，各100次，共200次复制，所以进入数据区域200次。但是在运行时反馈信息中向设备复制数据（H2D）的次数为400次，是因为数组的长度较长，于是将每个数组分两次传输。进入数据区域的次数并不代表H2D或者D2H的复制数据的次数，因为运行时数组长度的变化会导致复制的次数不同。第67行是主体迭代的出口，将数据从设备内存复制回主机，复制次数600次，这就体现了运行时数据的长度的变化问题。虽然数据长度有变化，每次消耗的时间不同，但是总的复制次数乘以平均时间还是与总的复制数据的时间吻合的。
2.2 Data构件优化Jacobi迭代数据传输
//Jacobi_GPU_noerror_data.c
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#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))
typedef struct timeval timestruct;
#define Mx 8191 //定义网格尺寸Mx和Ny，此处为了数组的内存对齐，所以设为8192-1
#define Ny 1023 //1024-1
float uval(float x, float y){return (x*x+y*y);}
int main()
{
float Widthy = 2.0, Heightx = 1.0; //Widthy为Ω=(0,W)x(0,H)中的W和H
//添加restrict关键字，使得u0和u1指向的地址空间不交叉
float *restrict u0, *restrict u1; //定义两个一维数组u0和u1，来保存二维格点上的旧值和新值
float hx, hy, hx2, hy2, fij, c1, c2;
float *tprt;
int ix, jy, maxIter, iter;
int Nyp1 = Ny+1;
timestruct t1, t2;
long long telapsed;
u0 = (float*)malloc(sizeof(float)*(Mx+1)*(Ny+1));
u1 = (float*)malloc(sizeof(float)*(Mx+1)*(Ny+1));
maxIter = 100; //指定最大迭代次数
hx = Heightx/Mx;
hy = Widthy/Ny;
//初始化u0/u1的左右边界，这些边界再计算中保持不变
for(ix = 0; ix <= Mx; ix++)
{
u0[ix*Nyp1+0] = u1[ix*Nyp1+0] = uval(ix*hx, 0.0f);
u0[ix*Nyp1+Ny] = u1[ix*Nyp1+Ny] = uval(ix*hx, Ny*Ny);
}
// 初始化u0/u1的上下边界，固定值保持不变
for(jy = 0; jy <= Ny; jy++)
{
u0[jy]=u1[jy]=uval(0.0f, jy*hy);
u0[Mx*Nyp1+jy] = u1[Mx*Nyp1+jy] = uval(Mx*hx, jy*hy);
}
//为数组u0内部元素赋初值，这里指定为0
for(ix = 1; ix < Mx; ix ++)
for(jy = 1; jy < Ny; jy++)
{
u0[ix*Nyp1+jy]=0.0f;
}
fij = -4.0f;
c1 = hx*hx * hy*hy;
c2 = 1.0f/(2.0*(hx*hx+hy*hy));
hx2 = hx*hx;
hy2 = hy*hy;
gettime(&t1);
#pragma acc data copy(u0[0:(Mx+1)*(Ny+1)])\
copyin(u1[0:(Mx+1)*(Ny+1)]) //利用data构件创建一个大的数据区域，把计算构件上的数据传输操作剥离出来，统一管理。
//copy的功能是数据传递，进入计算构件时，将数组u0,u1从主机内存复制到设备内存上，为并行计算做准备，计算完成后在将数组复制回主机内存。
//copyin 当进入kernels构件区域时，如果设备上已经为u1变量开辟了内存，则不需要将数据传输过去,没有开辟内存，则开辟内存并将u1传输至设备内存;
//在离开kernels构件区域时，不进行任何复制操作，并释放设备内存
{
//开始迭代，每一轮迭代中用u0更新u1同时计算两者间的最大误差值uerr，然后将u1中的新值赋给u0，准备进行下一步迭代
for(iter = 1; iter <= maxIter; iter++)
{
#pragma acc kernels present(u0[0:(Mx+1)*(Ny+1)]) \
present(u1[0:(Mx+1)*(Ny+1)])
//present子语，表明在进入构件区域之前，变量u0和u1已在设备上开辟内存。
{
#pragma acc loop independent
for(ix = 1; ix < Mx; ix++)
for(jy = 1; jy < Ny; jy++)
{
u1[ix*Nyp1+jy]=(c1*fij+hy2*(u0[(ix-1)*Nyp1+jy] + u0[(ix+1)*Nyp1+jy]) \
+ hx2*(u0[ix*Nyp1+jy-1] + u0[ix*Nyp1+jy+1]))*c2;
}
}
tprt = u0; u0 = u1; u1 = tprt;
}
}
gettime(&t2);
telapsed = usec(t1, t2);
printf("ElapsedTime = %13lld microseconds\n", telapsed); //打印循环主体消耗的时间
free(u0);
free(u1);
return 0;
}
执行代码：
pgcc -acc -Minfo Jacobi_GPU_noerror_data.c -o Jacobi_GPU_noerror_data.exe
./Jacobi_GPU_noerror_data.c 执行结果：
如上图可以看到，使用了data构件，62行的数据区域入口处，数据复制次数只在内核启动前进入数据区域2次，将u0和u1复制到设备，copyin复制过程即H2D过程传输4次，理由同上，都是u0和u1数组较大，每个分2次传输。此时数据传输的时间大大减少。如下表所示。
表4. GPU运行阶段data构件优化耗时（不复制u1回主机内存）
| 程序运行阶段 | 耗时（微秒） |
| H2D | 21304 |
| D2H | 10011 |
| 内核计算（墙上时间），不计算误差 | 884030 |
69行，则是内核迭代100次时，进入data构件的数据区域200次，分别计算u0和u1，故需200次。此时数据已经在设备内存上了。
3. 总结
3.1 CPU并行也可以尝试，来比较，不过是属于openMP的内容，本文就没有展开讨论了。
3.2 基于指令的移植：
使用OpenACC或者OpenHMPP提供的指令集添加到你想使用GPU计算的源代码中的某个位置，让编译器来编译出GPU上执行的代码，这种方式就是基于指令的移植方式。
例如，在可以GPU并行的代码处添加OpenACC导语 “ #pragma acc kernels ” ，这个导语的意思就是利用kernels构件来加速执行。
3.3 网格化（gridification）
如上文中，编译器分析了#pragma acc kernels下面的那个二重循环，就会根据循环的次数来分配线程数量，这个过程就叫网格化。
OpenACC与CUDA对应的名词称谓，gang是指block，vector是指block内的线程数目。
参考文献
- Jacobi迭代法_阿桑的专栏-CSDN博客_jacobi迭代法
- 2.2.1 二维泊松方程 - 51CTO.COM
- OpenACC例子 - 青竹居士 - 博客园 (cnblogs.com)
- 《OpenACC并行编程实战》
- 基于指令的移植方式的几个重要概念的理解（OpenHMPP, OpenACC）-转载 - 青竹居士 - 博客园 (cnblogs.com)
扫一扫
457
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
