高性能计算-粒子状态模拟计算-性能优化(19)

1. 源码为对粒子移动状态模拟的项目。要求使用多种优化方法,对比串行优化、多线程优化、全部优化下的加速比。

2. 代码

项目代码地址:https://github.com/libo-0379/StellarSim_Optimize

以下为核心优化代码及分析

#include <stdlib.h>
#include <iostream>
#include <math.h>
#include "../headers/global.h"
#include "../headers/sphkernel.h"
using namespace std;
//计算距离
void getPairwiseSeparations(double** &pos)
{
// 1. 提出循环不变量
// 2. todo simd
// 3. 对 j 进行循环分块 影响dx dy dz的命中,不适用
// 4. omp
#if defined(OPT_BASE) && (defined(OPT_SIMD)||defined(OPT_OMP))
#ifdef OPT_OMP
#pragma omp parallel for schedule(guided) proc_bind(close)
#endif
for (int i = 0; i < N; i++)
{
#ifndef OPT_SIMD
double temp1 = pos[0][i];
double temp2 = pos[1][i];
double temp3 = pos[2][i];
for (int j = 0; j < N; j++)
{
// dx[i][j] = -dx[j][i] 粒子彼此计算相对距离
dx[i][j] = pos[0][j] - temp1;
dy[i][j] = pos[1][j] - temp2;
dz[i][j] = pos[2][j] - temp3;
}
#else
float64x2_t v0 = vld1q_dup_f64(&pos[0][i]);
float64x2_t v1 = vld1q_dup_f64(&pos[1][i]);
float64x2_t v2 = vld1q_dup_f64(&pos[2][i]);
for(int j=0;j<N/2*2;j+=2)
{
float64x2_t v0_0 = vld1q_f64(&pos[0][j]);
float64x2_t v1_0 = vld1q_f64(&pos[1][j]);
float64x2_t v2_0 = vld1q_f64(&pos[2][j]);
vst1q_f64(&dx[i][j],vsubq_f64(v0_0,v0));
vst1q_f64(&dy[i][j],vsubq_f64(v1_0,v1));
vst1q_f64(&dz[i][j],vsubq_f64(v2_0,v2));
}
for (int j = N/2*2; j < N; j++)
{
dx[i][j] = pos[0][j] - pos[0][i];
dy[i][j] = pos[1][j] - pos[1][i];
dz[i][j] = pos[2][j] - pos[2][i];
}
#endif
}
#else
#ifdef OPT_OMP
#pragma omp parallel for schedule(guided) proc_bind(close)
#endif
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
// dx[i][j] = -dx[j][i] 粒子彼此计算相对距离
dx[i][j] = pos[0][j] - pos[0][i];
dy[i][j] = pos[1][j] - pos[1][i];
dz[i][j] = pos[2][j] - pos[2][i];
//fprintf(stdout, "%12.6f", dz[i][j]);
//fflush(stdout);
}
//fprintf(stdout,"\n");
}
#endif
}
void getW(double** &dx, double** &dy, double** &dz, const double h)
{
// 1. 循环不变量提出
// 2. omp
#if defined(OPT_OMP) || defined(OPT_BASE)
double value1 = pow((1.0 / (h*sqrt(pi))), 3.0);
double value2 = pow(h,2);
#ifdef OPT_OMP
#pragma omp parallel for schedule(guided) proc_bind(close)
#endif
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
r[i][j] = sqrt(pow(dx[i][j],2.0) + pow(dy[i][j],2.0) + pow(dz[i][j],2.0));
W[i][j] = value1 * exp((-pow(r[i][j],2) / value2));
}
}
#else
#ifdef OPT_OMP
#pragma omp parallel for schedule(guided) proc_bind(close)
#endif
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
r[i][j] = sqrt(pow(dx[i][j],2.0) + pow(dy[i][j],2.0) + pow(dz[i][j],2.0));
W[i][j] = pow((1.0 / (h*sqrt(pi))), 3.0) * exp((-pow(r[i][j],2) / pow(h,2)));
//fprintf(stdout, "%12.6f", r[i][j]);
//fprintf(stdout, "%12.6f", W[i][j]);
//fflush(stdout);
}
}
#endif
//fprintf(stdout,"\n");
}
void getGradW(double** &dx, double** &dy, double** &dz, const double h)
{
// 1. 循环不变量提出
// 2. omp
#if defined(OPT_OMP) || defined(OPT_BASE)
double value1 = pow(h,2);
double value2 = -2/pow(h,5)/pow(pi,(3/2));
#ifdef OPT_OMP
#pragma omp parallel for schedule(guided) proc_bind(close)
#endif
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
r[i][j] = sqrt(pow(dx[i][j],2.0) + pow(dy[i][j],2.0) + pow(dz[i][j],2.0));
gradPara[i][j] = exp(-pow(r[i][j],2) / value1) * value2;
wx[i][j] = gradPara[i][j]*dx[i][j];
wy[i][j] = gradPara[i][j]*dy[i][j];
wz[i][j] = gradPara[i][j]*dz[i][j];
}
}
#else
#ifdef OPT_OMP
#pragma omp parallel for schedule(guided) proc_bind(close)
#endif
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
// r[i][j] = r[j][i]
r[i][j] = sqrt(pow(dx[i][j],2.0) + pow(dy[i][j],2.0) + pow(dz[i][j],2.0));
// gradPara[i][j] = gradPara[j][i]
gradPara[i][j] = -2 * exp(-pow(r[i][j],2) / pow(h,2)) / pow(h,5) / pow(pi,(3/2));
// wx[i][j] = -wx[j][i]
wx[i][j] = gradPara[i][j]*dx[i][j];
wy[i][j] = gradPara[i][j]*dy[i][j];
wz[i][j] = gradPara[i][j]*dz[i][j];
//fprintf(stdout, "%12.6f", wy[i][j]);
//fflush(stdout);
}
//fprintf(stdout,"\n");
}
#endif
}
void getDensity(double** &pos, double &m, const double h)
{
getPairwiseSeparations(pos);
getW(dx, dy, dz, h);
// 1. 改变访问顺序
// 2. 内层可 simd
// 3. 外层omp有数据竞争,内层omp伪共享(或者使用cacheline大小为步进), omp 不适用
// 4. 简化计算公式 rho[j] += W[i][j] rho[i] *= m;W 每一列计算出一个 rho[j] 此处不适用
#ifdef OPT_BASE
for(int i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
rho[j] += m * W[i][j];
}
#else
for (int j = 0; j < N; j++)
{
for (int i = 0; i < N; i++)
rho[j] += m * W[i][j];
//fprintf(stdout, "%12.6f", rho[j]);
//fflush(stdout);
//fprintf(stdout,"\n");
}
#endif
}
void getPressure(double* &rho, const double k, double &n)
{
// 1. 提出循环不变量,循环展开
// 2. omp 线程调度开销大,不合适
// 3. simd pow较复杂
#ifdef OPT_BASE
double value = 1+1/n;
for (int j = 0; j < N; j++)
P[j] = k * pow(rho[j], value);
#else
for (int j = 0; j < N; j++)
{
P[j] = k * pow(rho[j], (1+1/n));
//fprintf(stdout, "%12.6f\n", P[j]);
}
#endif
}
void getAcc(double** &pos, double** &vel, double &m, const double h, const double k, double &n, double lmbda, const double nu)
{
getDensity(pos, m, h);
getPressure(rho, k, n);
getPairwiseSeparations(pos);
getGradW(dx, dy, dz, h);
#if defined(OPT_BASE)
#ifdef OPT_OMP
#pragma omp parallel for schedule(guided) proc_bind(close)
#endif
for (int j = 0; j < N; j++)
{
// 1. wx[i][j] = -wx[j][i] 访问 w[j][i] 可以增加缓存命中,并且可以向量化
// 如果for循环交换 j i访问顺序,先访问i 后 j,内层循环无法做向量化优化(内循环每次计算不同的目标元素)
// 并且,如果for循环先访问j,计算 acc[0][j], P[j]/pow(rho[j],2)是常量可以提出最后计算
// 2. 简化计算 m* 放在求和之后
// 3. 循环不变量提出
double temp1 = P[j]/pow(rho[j],2);
double temp3 =0.0,temp4=0.0,temp5=0.0;
for (int i = 0; i < N; i++)
{
double temp2 = pow(P[i]/rho[i],2);
temp3 += (temp1 + temp2) * wx[j][i];
temp4 += (temp1 + temp2) * wy[j][i];
temp5 += (temp1 + temp2) * wz[j][i];
}
acc[0][j] += (temp3 *=m);
acc[1][j] += (temp4 *=m);
acc[2][j] += (temp5 *=m);
}
#else
#ifdef OPT_OMP
#pragma omp parallel for schedule(guided) proc_bind(close)
#endif
for (int j = 0; j < N; j++)
{
for (int i = 0; i < N; i++)
{
acc[0][j] -= m * ( P[j]/pow(rho[j],2) + pow(P[i]/rho[i],2) ) * wx[i][j];
acc[1][j] -= m * ( P[j]/pow(rho[j],2) + pow(P[i]/rho[i],2) ) * wy[i][j];
acc[2][j] -= m * ( P[j]/pow(rho[j],2) + pow(P[i]/rho[i],2) ) * wz[i][j];
}
}
#endif
// 1. simd
// 2. 循环合并
#ifdef OPT_BASE
for (int j = 0; j < N; j++)
{
acc[0][j] -= (lmbda * pos[0][j] + nu * vel[0][j]);
acc[1][j] -= (lmbda * pos[1][j] + nu * vel[1][j]);
acc[2][j] -= (lmbda * pos[2][j] + nu * vel[2][j]);
}
#else
for (int j = 0; j < N; j++)
{
acc[0][j] -= lmbda * pos[0][j];
acc[1][j] -= lmbda * pos[1][j];
acc[2][j] -= lmbda * pos[2][j];
}
for (int j = 0; j < N; j++)
{
acc[0][j] -= nu * vel[0][j];
acc[1][j] -= nu * vel[1][j];
acc[2][j] -= nu * vel[2][j];
}
#endif
}
#ifdef OPT_SIMD
// 需要用到 泰勒展开
float64_t exp_(float64_t x)
{
//初始化第一个值
int n = 0;
double prior = 1.0;
double sum = prior; //求和保存结果
while(1)
{
double cur = prior * x /++n;
sum += cur;
prior = cur;
if(cur<=EPSILON)
break;
}
return sum;
}
// a^b = e^(b*ln(a)); neon 未提供ln,设想采用 cmath ln函数,向量化对每个元素的计算用 omp task
// float64_t pow_(float64_t a,float64_t b)
// {
// logf()
// }
#endif

3. 测试数据

3.1 所有编译优化选项为 O2,不开启向量自动化优化。代码内函数计时

项目 耗时 s 相比源码加速比
original(源码) 141
BASE(循环优化) 58.7 2.4
BASE+SIMD 51 2.8
BASE+OMP 10 14.1
OMP 23.1 6.1
BASE+OMP+SIMD 6.47 21.8

3.2 gprof 耗时分析

original

image

BASE

image

BASE_SIMD

image

BASE_OMP

image

OMP

image

BASE_SIMD_OMP

image

4. 结果分析

4.1 程序内计时

(1) 单核串行耗时从 141s 优化为 51s,加速比为 2.8;

(2) 多核32线程耗时 23.1s(实测16线程耗时一致),加速比为 6.1;

(3) 综合优化后耗时 6.47,加速比为 21.8。

4.2 其他方面

(1) 多线程优化方面 schedule(guided) 策略拥有最高的效率,比 dynamic 略优。

(2) 基础优化和多线程有较大的加速比提升。

(3) 16线程与32线程的效率一致。

(4) 注意多线程数据竞争问题,可能由于此问题造成多线程效率下降。

(5) 注意cacheline 竞争问题,避免多线程效率下降和计算错误。

(6) 注意在不同测试项下确认计算结果的正确性,结算结果应保持一致。

4.3 gprof分析

(1) 使用openmp 会明显增加线程框架开销,每增加一个线程的并行使用会增加相应的线程开销;

(2) getAcc 函数从源码的时间从 36.9s 最终将为几乎为 0;

image

(3) 函数计算的耗时比例从 100% 降为 15.9%,其余 83.78% 为多线程和框架管理带来的开销。

posted @   安洛8  阅读(16)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· DeepSeek R1 简明指南:架构、训练、本地部署及硬件要求
· NetPad:一个.NET开源、跨平台的C#编辑器
· PowerShell开发游戏 · 打蜜蜂
点击右上角即可分享
微信分享提示