StellarSim项目性能优化

已于 2024-06-17 10:45:21 修改
阅读量244 收藏 5
点赞数 10
文章标签： c语言 算法
于 2024-03-12 20:48:50 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Begoniaish/article/details/136662265
版权
StellerSim 是一个基于离散单元数值计算的物理演化程序，该 算法最早从天体动力学星云演化计算中产生，被广泛用于游戏开 发引擎中的物理演化、碰撞事件、特效生成等等模块。

1.下载解压StellarSim压缩文件，进入安装目录，编辑Makefile，在LDFLAGS参数后添加-pg

CC = g++
CFLAGS = -Wall -Wextra -O1 -Iheaders 


SRC = $(wildcard src/*.cpp)
OBJ = $(SRC:.cpp=.o)
EXEC = stellarSim
LDFLAGS = -lm -pg

all: $(EXEC)

$(EXEC): $(OBJ)
	$(CC)  $^ -o $@  $(LDFLAGS)        
    $(CC) $(CFLAGS) $(OBJECTS) -o $@

%.o: %.cpp
	$(CC) $(CFLAGS) -c $< -o $@ 
    $(CC) $(CFLAGS) $(INC) -c $< -o $@

clean:
	rm -f $(OBJ) $(EXEC)

2.运行make命令，编译StellarSim程序

3.使用如下命令运行程序。

yhrun -p thcp1 ./stellarSim 1024

 

4.使用gprof工具分析StellarSim程序，找到时间开销最大的核心子函数。

可以看到耗时最长的是getAcc函数，时间45.88秒，占比43.88%，getDensity第二，占比17.47....

5.找到这些函数

void getPairwiseSeparations(double** &pos)
{
  for (int i = 0; i < N; i++) {
    for (int j = 0; 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];
      //fprintf(stdout, "%12.6f", dz[i][j]);
      //fflush(stdout);
    }
    //fprintf(stdout,"\n");
  }
}

void getW(double** &dx, double** &dy, double** &dz, const double h)
{
  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);
    }
    //fprintf(stdout,"\n");
  }
}

void getGradW(double** &dx, double** &dy, double** &dz, const double h)
{
  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] = -2 * exp(-pow(r[i][j],2) / pow(h,2)) / pow(h,5) / pow(pi,(3/2));

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

void getDensity(double** &pos, double &m, const double h)
{
  getPairwiseSeparations(pos);
  getW(dx, dy, dz, h);
  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[i]);
    //fflush(stdout);
    //fprintf(stdout,"\n");
  }
}


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

  for (int i = 0; i < N; i++) {
	for (int j = 0; j < N; j++) {
      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];
    }
  }
  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];
  }
}

6.要优化getAcc函数的性能，可以考虑以下几个方面：

  1. 减少重复计算：在循环中存在大量重复的计算，例如重复计算 P[j]/pow(rho[j],2) 和 P[i]/rho[i]。您可以在循环外部进行这些计算，并将结果存储起来，以避免重复计算。

  2. 使用并行化：考虑循环的并行化处理，特别是第一个和第二个循环。您可以使用并行计算库（如OpenMP）来加速计算过程。

  3. 数据结构优化：检查数据结构的使用方式，确保内存访问模式是高效的。例如，可以尝试重新组织数据结构，以便更好地利用缓存和减少内存访问次数。

  4. 函数调用开销：如果 getDensity、getPressure、getPairwiseSeparations、getGradW 等函数具有较大的开销，可以考虑将它们的计算结果缓存起来，避免重复调用。

  5. 参数传递优化：考虑对函数参数进行优化，避免不必要的参数传递，尤其是对于数组和指针类型的参数。

  6. 算法优化：根据具体的物理模型和计算需求，考虑是否存在更高效的数值计算算法。

  7. 编译器优化：确保使用适当的编译器选项进行优化，例如启用优化标志 -O3

7.接下来采用 OpenMP多线程并行优化

#include <omp.h>

void getPairwiseSeparations(double** &pos)
{
  #pragma omp parallel for collapse(2)
  for (int i = 0; i < N; i++) {
    for (int j = 0; 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];
      //fprintf(stdout, "%12.6f", dz[i][j]);
      //fflush(stdout);
    }
    //fprintf(stdout,"\n");
  }
}

void getW(double** &dx, double** &dy, double** &dz, const double h)
{
  #pragma omp parallel for collapse(2)
  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);
    }
    //fprintf(stdout,"\n");
  }
}

void getGradW(double** &dx, double** &dy, double** &dz, const double h)
{
  #pragma omp parallel for collapse(2)
  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] = -2 * exp(-pow(r[i][j],2) / pow(h,2)) / pow(h,5) / pow(pi,(3/2));

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

void getDensity(double** &pos, double &m, const double h)
{
  getPairwiseSeparations(pos);
  getW(dx, dy, dz, h);
  #pragma omp parallel
  for (int i = 0; i < N; i++) {
	#pragma omp  for //reduction(+:rho[j])
    for (int j = 0; j < N; j++) {
      rho[j] += m * W[i][j];
    }
    //fprintf(stdout, "%12.6f", rho[j]);
   // fflush(stdout);
    //fprintf(stdout,"\n");
  }
}


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);
  #pragma omp parallel 
  {
  for (int i = 0; i < N; i++) {
	#pragma omp for 
	for (int j = 0; j < N; j++) {
      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];
    }
  }
  #pragma omp barrier
  #pragma omp for
  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];

    acc[0][j] -= nu * vel[0][j];
    acc[1][j] -= nu * vel[1][j];
    acc[2][j] -= nu * vel[2][j];
  }
  }
}

 8.接下来采用 SIMD向量化并行计算优化

#include <arm_neon.h>

// 使用C/C++函数计算exp()函数
float64x2_t neon_exp(float64x2_t input)
{
    // 计算指数值
    float64x2_t result;
    result[0] = exp(input[0]);
    result[1] = exp(input[1]);

    return result;
}

void getPairwiseSeparations(double** &pos)
{
	float64x2_t dx_i_j,dy_i_j,dz_i_j;
	float64x2_t pos_0_j, pos_1_j, pos_2_j;
	float64x2_t pos_0_i, pos_1_i, pos_2_i;
	for (int i = 0; i < N ; i ++) {
		pos_0_i = vdupq_n_f64(pos[0][i]);
		pos_1_i = vdupq_n_f64(pos[1][i]);
		pos_2_i = vdupq_n_f64(pos[2][i]);
		for (int j = 0; j < N/ 2 * 2; j += 2) {
			pos_0_j = vld1q_f64(&pos[0][j]);
			pos_1_j = vld1q_f64(&pos[1][j]);
			pos_2_j = vld1q_f64(&pos[2][j]);
			
			dx_i_j = vld1q_f64(&dx[i][j]);
			dy_i_j = vld1q_f64(&dy[i][j]);
			dz_i_j = vld1q_f64(&dz[i][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];
			dx_i_j = vsubq_f64(pos_0_j,pos_0_i);
			dy_i_j = vsubq_f64(pos_1_j,pos_1_i);
			dz_i_j = vsubq_f64(pos_2_j,pos_2_i);			
			
			vst1q_f64(&dx[i][j], dx_i_j);
            vst1q_f64(&dy[i][j], dy_i_j);
            vst1q_f64(&dz[i][j], dz_i_j);
			//fprintf(stdout, "%12.6f", dz[i][j]);
			//fflush(stdout);
		}
    //fprintf(stdout,"\n");
	}
}

void getW(double** &dx, double** &dy, double** &dz, const double h)
{
	float tmp1 = pow((1.0 / (h*sqrt(pi))), 3.0);
	float tmp2 = -1 / pow(h,2);
	float64x2_t	tmp1Vec = vdupq_n_f64(tmp1);
	float64x2_t tmp2Vec = vdupq_n_f64(tmp2);
	float64x2_t r_i_j,W_i_j; 
	float64x2_t dx_i_j,dy_i_j,dz_i_j;
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N/ 2 * 2; j += 2)  {
			dx_i_j = vld1q_f64(&dx[i][j]);
			dy_i_j = vld1q_f64(&dy[i][j]);
			dz_i_j = vld1q_f64(&dz[i][j]);
			
			r_i_j  = vld1q_f64(&r[i][j]);
			W_i_j = vld1q_f64(&W[i][j]);
			//r[i][j] = sqrt(pow(dx[i][j],2.0) + pow(dy[i][j],2.0) + pow(dz[i][j],2.0));
			r_i_j = vsqrtq_f64(vaddq_f64(vmulq_f64(dx_i_j, dx_i_j),vaddq_f64(vmulq_f64(dy_i_j, dy_i_j),vmulq_f64(dz_i_j, dz_i_j))));
			//W[i][j] = pow((1.0 / (h*sqrt(pi))), 3.0) * exp((-pow(r[i][j],2) / pow(h,2))); 
			W_i_j = vmulq_f64(tmp1Vec,neon_exp(vmulq_f64(vmulq_f64(r_i_j,r_i_j),tmp2Vec)));
			
			vst1q_f64(&r[i][j], r_i_j);
			vst1q_f64(&W[i][j], W_i_j);
		//fprintf(stdout, "%12.6f", r[i][j]);
		//fprintf(stdout, "%12.6f", W[i][j]);
		//fflush(stdout);
		}
    //fprintf(stdout,"\n");
	}
}


void getGradW(double** &dx, double** &dy, double** &dz, const double h)
{
	
	float tmp1 = -1 / pow(h,2);
	float tmp2 = -2 / pow(h,5) / pow(pi,(3/2));
	float64x2_t	tmp1Vec = vdupq_n_f64(tmp1);
	float64x2_t tmp2Vec = vdupq_n_f64(tmp2);
	float64x2_t dx_i_j,dy_i_j,dz_i_j;
	float64x2_t r_i_j,gradPara_i_j;
	float64x2_t wx_i_j,wy_i_j,wz_i_j;
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N/ 2 * 2; j += 2) {
			dx_i_j = vld1q_f64(&dx[i][j]);
			dy_i_j = vld1q_f64(&dy[i][j]);
			dz_i_j = vld1q_f64(&dz[i][j]);
			
			r_i_j  = vld1q_f64(&r[i][j]);
			gradPara_i_j = vld1q_f64(&gradPara[i][j]);
			
			wx_i_j = vld1q_f64(&wx[i][j]);
            wy_i_j = vld1q_f64(&wy[i][j]);
            wz_i_j = vld1q_f64(&wz[i][j]);
			
			//r[i][j]  = sqrt(pow(dx[i][j],2.0) + pow(dy[i][j],2.0) + pow(dz[i][j],2.0));
			r_i_j = vsqrtq_f64(vaddq_f64(vmulq_f64(dx_i_j, dx_i_j),vaddq_f64(vmulq_f64(dy_i_j, dy_i_j),vmulq_f64(dz_i_j, dz_i_j))));
			//gradPara[i][j] = -2 * exp(-pow(r[i][j],2) / pow(h,2)) / pow(h,5) / pow(pi,(3/2));
			gradPara_i_j = vmulq_f64(neon_exp(vmulq_f64(vmulq_f64(r_i_j,r_i_j),tmp1Vec)),tmp2Vec);
			//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];		
			wx_i_j = vmulq_f64(gradPara_i_j,dx_i_j);
			wy_i_j = vmulq_f64(gradPara_i_j,dy_i_j);
			wz_i_j = vmulq_f64(gradPara_i_j,dz_i_j);
			
			vst1q_f64(&wx[i][j], wx_i_j);
            vst1q_f64(&wy[i][j], wy_i_j);
            vst1q_f64(&wz[i][j], wz_i_j);
			//fprintf(stdout, "%12.6f", wy[i][j]);
			//fflush(stdout);
		}
		//fprintf(stdout,"\n");
	}
}

void getDensity(double** &pos, double &m, const double h)
{
	float64x2_t	mVec = vdupq_n_f64(m);
	float64x2_t	W_i_j,rho_j;
	getPairwiseSeparations(pos);
	getW(dx, dy, dz, h);
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N/ 2 * 2; j += 2) {
			rho_j = vld1q_f64(&rho[j]);
			W_i_j = vld1q_f64(&W[i][j]);
			rho_j = vfmaq_f64(rho_j,mVec,W_i_j);
			vst1q_f64(&rho[j], rho_j);
		}
    //fprintf(stdout, "%12.6f", rho[i]);
    //fflush(stdout);
    //fprintf(stdout,"\n");
	}
}


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

	float64x2_t mVec = vdupq_n_f64(m);
    float64x2_t lmbdaVec = vdupq_n_f64(lmbda);
    float64x2_t nuVec = vdupq_n_f64(nu);
	float64x2_t acc_0_j, acc_1_j, acc_2_j;
	float64x2_t P_i,rho_i,P_i_div_rho_i_sq;
	float64x2_t P_j,rho_j,P_j_div_rho_j_sq;
	float64x2_t wx_i_j,wy_i_j,wz_i_j;
    for (int i = 0; i < N; i++) {
        P_i = vdupq_n_f64(P[i]);
        rho_i = vdupq_n_f64(rho[i]);
        P_i_div_rho_i_sq = vmulq_f64(vmulq_f64(P_i, P_i), vrecpeq_f64(vmulq_f64(rho_i, rho_i)));
        for (int j = 0; j < N/ 2 * 2; j += 2) {
			P_j = vld1q_f64(&P[j]);
            rho_j = vld1q_f64(&rho[j]);
            P_j_div_rho_j_sq = vmulq_f64(P_j, vrecpeq_f64(vmulq_f64(rho_j, rho_j)));
            acc_0_j = vld1q_f64(&acc[0][j]);
            acc_1_j = vld1q_f64(&acc[1][j]);
            acc_2_j = vld1q_f64(&acc[2][j]);

			
            wx_i_j = vld1q_f64(&wx[i][j]);
            wy_i_j = vld1q_f64(&wy[i][j]);
            wz_i_j = vld1q_f64(&wz[i][j]);
			//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];
            acc_0_j = vfmsq_f64(acc_0_j, vmulq_f64(mVec,wx_i_j), vaddq_f64(P_j_div_rho_j_sq, P_i_div_rho_i_sq));
            acc_1_j = vfmsq_f64(acc_1_j, vmulq_f64(mVec,wy_i_j), vaddq_f64(P_j_div_rho_j_sq, P_i_div_rho_i_sq));
            acc_2_j = vfmsq_f64(acc_2_j, vmulq_f64(mVec,wz_i_j), vaddq_f64(P_j_div_rho_j_sq, P_i_div_rho_i_sq));
			
			vst1q_f64(&acc[0][j], acc_0_j);
            vst1q_f64(&acc[1][j], acc_1_j);
            vst1q_f64(&acc[2][j], acc_2_j);
			

        }
    }
	for (int j = 0; j < N /2 * 2; j += 2) {
		acc_0_j = vld1q_f64(&acc[0][j]);
        acc_1_j = vld1q_f64(&acc[1][j]);
        acc_2_j = vld1q_f64(&acc[2][j]);
		//acc[0][j] -= lmbda * pos[0][j]; 
		//acc[1][j] -= lmbda * pos[1][j];
		//acc[2][j] -= lmbda * pos[2][j];
	    acc_0_j = vfmsq_f64(acc_0_j, lmbdaVec, vld1q_f64(&pos[0][j]));
        acc_1_j = vfmsq_f64(acc_1_j, lmbdaVec, vld1q_f64(&pos[1][j]));
        acc_2_j = vfmsq_f64(acc_2_j, lmbdaVec, vld1q_f64(&pos[2][j]));
		//acc[0][j] -= nu * vel[0][j];
		//acc[1][j] -= nu * vel[1][j];
		//acc[2][j] -= nu * vel[2][j];
        acc_0_j = vfmsq_f64(acc_0_j, nuVec, vld1q_f64(&vel[0][j]));
        acc_1_j = vfmsq_f64(acc_1_j, nuVec, vld1q_f64(&vel[1][j]));
        acc_2_j = vfmsq_f64(acc_2_j, nuVec, vld1q_f64(&vel[2][j]));

        vst1q_f64(&acc[0][j], acc_0_j);
        vst1q_f64(&acc[1][j], acc_1_j);
        vst1q_f64(&acc[2][j], acc_2_j);
	}
}

9.最后我们将两种优化方式合并

#include <arm_neon.h>
#include <omp.h>

// 使用C/C++函数计算exp()函数
float64x2_t neon_exp(float64x2_t input)
{
    // 计算指数值
    float64x2_t result;
    result[0] = exp(input[0]);
    result[1] = exp(input[1]);

    return result;
}

void getPairwiseSeparations(double** &pos)
{
	float64x2_t dx_i_j,dy_i_j,dz_i_j;
	float64x2_t pos_0_j, pos_1_j, pos_2_j;
	float64x2_t pos_0_i, pos_1_i, pos_2_i;
	#pragma omp parallel for
	for (int i = 0; i < N ; i ++) {
		pos_0_i = vdupq_n_f64(pos[0][i]);
		pos_1_i = vdupq_n_f64(pos[1][i]);
		pos_2_i = vdupq_n_f64(pos[2][i]);		
		for (int j = 0; j < N/ 2 * 2; j += 2) {
			pos_0_j = vld1q_f64(&pos[0][j]);
			pos_1_j = vld1q_f64(&pos[1][j]);
			pos_2_j = vld1q_f64(&pos[2][j]);
			
			dx_i_j = vld1q_f64(&dx[i][j]);
			dy_i_j = vld1q_f64(&dy[i][j]);
			dz_i_j = vld1q_f64(&dz[i][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];
			dx_i_j = vsubq_f64(pos_0_j,pos_0_i);
			dy_i_j = vsubq_f64(pos_1_j,pos_1_i);
			dz_i_j = vsubq_f64(pos_2_j,pos_2_i);			
			
			vst1q_f64(&dx[i][j], dx_i_j);
            vst1q_f64(&dy[i][j], dy_i_j);
            vst1q_f64(&dz[i][j], dz_i_j);
			//fprintf(stdout, "%12.6f", dz[i][j]);
			//fflush(stdout);
		}
    //fprintf(stdout,"\n");
	}
}

void getW(double** &dx, double** &dy, double** &dz, const double h)
{
	float tmp1 = pow((1.0 / (h*sqrt(pi))), 3.0);
	float tmp2 = -1 / pow(h,2);
	float64x2_t	tmp1Vec = vdupq_n_f64(tmp1);
	float64x2_t tmp2Vec = vdupq_n_f64(tmp2);
	float64x2_t r_i_j,W_i_j; 
	float64x2_t dx_i_j,dy_i_j,dz_i_j;
	#pragma omp parallel for 
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N/ 2 * 2; j += 2)  {
			dx_i_j = vld1q_f64(&dx[i][j]);
			dy_i_j = vld1q_f64(&dy[i][j]);
			dz_i_j = vld1q_f64(&dz[i][j]);
			
			r_i_j  = vld1q_f64(&r[i][j]);
			W_i_j = vld1q_f64(&W[i][j]);
			//r[i][j] = sqrt(pow(dx[i][j],2.0) + pow(dy[i][j],2.0) + pow(dz[i][j],2.0));
			r_i_j = vsqrtq_f64(vaddq_f64(vmulq_f64(dx_i_j, dx_i_j),vaddq_f64(vmulq_f64(dy_i_j, dy_i_j),vmulq_f64(dz_i_j, dz_i_j))));
			//W[i][j] = pow((1.0 / (h*sqrt(pi))), 3.0) * exp((-pow(r[i][j],2) / pow(h,2))); 
			W_i_j = vmulq_f64(tmp1Vec,neon_exp(vmulq_f64(vmulq_f64(r_i_j,r_i_j),tmp2Vec)));
			
			vst1q_f64(&r[i][j], r_i_j);
			vst1q_f64(&W[i][j], W_i_j);
		//fprintf(stdout, "%12.6f", r[i][j]);
		//fprintf(stdout, "%12.6f", W[i][j]);
		//fflush(stdout);
		}
    //fprintf(stdout,"\n");
	}
}


void getGradW(double** &dx, double** &dy, double** &dz, const double h)
{
	
	float tmp1 = -1 / pow(h,2);
	float tmp2 = -2 / pow(h,5) / pow(pi,(3/2));
	float64x2_t	tmp1Vec = vdupq_n_f64(tmp1);
	float64x2_t tmp2Vec = vdupq_n_f64(tmp2);
	float64x2_t dx_i_j,dy_i_j,dz_i_j;
	float64x2_t r_i_j,gradPara_i_j;
	float64x2_t wx_i_j,wy_i_j,wz_i_j;
	#pragma omp parallel for 
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N/ 2 * 2; j += 2) {
			dx_i_j = vld1q_f64(&dx[i][j]);
			dy_i_j = vld1q_f64(&dy[i][j]);
			dz_i_j = vld1q_f64(&dz[i][j]);
			
			r_i_j  = vld1q_f64(&r[i][j]);
			gradPara_i_j = vld1q_f64(&gradPara[i][j]);
			
			wx_i_j = vld1q_f64(&wx[i][j]);
            wy_i_j = vld1q_f64(&wy[i][j]);
            wz_i_j = vld1q_f64(&wz[i][j]);
			
			//r[i][j]  = sqrt(pow(dx[i][j],2.0) + pow(dy[i][j],2.0) + pow(dz[i][j],2.0));
			r_i_j = vsqrtq_f64(vaddq_f64(vmulq_f64(dx_i_j, dx_i_j),vaddq_f64(vmulq_f64(dy_i_j, dy_i_j),vmulq_f64(dz_i_j, dz_i_j))));
			//gradPara[i][j] = -2 * exp(-pow(r[i][j],2) / pow(h,2)) / pow(h,5) / pow(pi,(3/2));
			gradPara_i_j = vmulq_f64(neon_exp(vmulq_f64(vmulq_f64(r_i_j,r_i_j),tmp1Vec)),tmp2Vec);
			//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];		
			wx_i_j = vmulq_f64(gradPara_i_j,dx_i_j);
			wy_i_j = vmulq_f64(gradPara_i_j,dy_i_j);
			wz_i_j = vmulq_f64(gradPara_i_j,dz_i_j);
			
			vst1q_f64(&wx[i][j], wx_i_j);
            vst1q_f64(&wy[i][j], wy_i_j);
            vst1q_f64(&wz[i][j], wz_i_j);
			//fprintf(stdout, "%12.6f", wy[i][j]);
			//fflush(stdout);
		}
		//fprintf(stdout,"\n");
	}
}

void getDensity(double** &pos, double &m, const double h)
{
	float64x2_t	mVec = vdupq_n_f64(m);
	float64x2_t	W_i_j,rho_j;
	getPairwiseSeparations(pos);
	getW(dx, dy, dz, h);
	//#pragma omp parallel for 这里不能进行omp优化会影响结果
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N/ 2 * 2; j += 2) {
			rho_j = vld1q_f64(&rho[j]);
			W_i_j = vld1q_f64(&W[i][j]);
			rho_j = vfmaq_f64(rho_j,mVec,W_i_j);
			vst1q_f64(&rho[j], rho_j);
		}
    //fprintf(stdout, "%12.6f", rho[i]);
    //fflush(stdout);
    //fprintf(stdout,"\n");
	}
}


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

	float64x2_t mVec = vdupq_n_f64(m);
    float64x2_t lmbdaVec = vdupq_n_f64(lmbda);
    float64x2_t nuVec = vdupq_n_f64(nu);
	float64x2_t acc_0_j, acc_1_j, acc_2_j;
	float64x2_t P_i,rho_i,P_i_div_rho_i_sq;
	float64x2_t P_j,rho_j,P_j_div_rho_j_sq;
	float64x2_t wx_i_j,wy_i_j,wz_i_j;
	//#pragma omp parallel for 这里不能进行omp优化会影响结果
    for (int i = 0; i < N; i++) {
        P_i = vdupq_n_f64(P[i]);
        rho_i = vdupq_n_f64(rho[i]);
        P_i_div_rho_i_sq = vmulq_f64(vmulq_f64(P_i, P_i), vrecpeq_f64(vmulq_f64(rho_i, rho_i)));
		
        for (int j = 0; j < N/ 2 * 2; j += 2) {
			P_j = vld1q_f64(&P[j]);
            rho_j = vld1q_f64(&rho[j]);
            P_j_div_rho_j_sq = vmulq_f64(P_j, vrecpeq_f64(vmulq_f64(rho_j, rho_j)));
            acc_0_j = vld1q_f64(&acc[0][j]);
            acc_1_j = vld1q_f64(&acc[1][j]);
            acc_2_j = vld1q_f64(&acc[2][j]);

			
            wx_i_j = vld1q_f64(&wx[i][j]);
            wy_i_j = vld1q_f64(&wy[i][j]);
            wz_i_j = vld1q_f64(&wz[i][j]);
			//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];
            acc_0_j = vfmsq_f64(acc_0_j, vmulq_f64(mVec,wx_i_j), vaddq_f64(P_j_div_rho_j_sq, P_i_div_rho_i_sq));
            acc_1_j = vfmsq_f64(acc_1_j, vmulq_f64(mVec,wy_i_j), vaddq_f64(P_j_div_rho_j_sq, P_i_div_rho_i_sq));
            acc_2_j = vfmsq_f64(acc_2_j, vmulq_f64(mVec,wz_i_j), vaddq_f64(P_j_div_rho_j_sq, P_i_div_rho_i_sq));
			
			vst1q_f64(&acc[0][j], acc_0_j);
            vst1q_f64(&acc[1][j], acc_1_j);
            vst1q_f64(&acc[2][j], acc_2_j);
			

        }
    }
	for (int j = 0; j < N /2 * 2; j += 2) {
		acc_0_j = vld1q_f64(&acc[0][j]);
        acc_1_j = vld1q_f64(&acc[1][j]);
        acc_2_j = vld1q_f64(&acc[2][j]);
		//acc[0][j] -= lmbda * pos[0][j]; 
		//acc[1][j] -= lmbda * pos[1][j];
		//acc[2][j] -= lmbda * pos[2][j];
	    acc_0_j = vfmsq_f64(acc_0_j, lmbdaVec, vld1q_f64(&pos[0][j]));
        acc_1_j = vfmsq_f64(acc_1_j, lmbdaVec, vld1q_f64(&pos[1][j]));
        acc_2_j = vfmsq_f64(acc_2_j, lmbdaVec, vld1q_f64(&pos[2][j]));
		//acc[0][j] -= nu * vel[0][j];
		//acc[1][j] -= nu * vel[1][j];
		//acc[2][j] -= nu * vel[2][j];
        acc_0_j = vfmsq_f64(acc_0_j, nuVec, vld1q_f64(&vel[0][j]));
        acc_1_j = vfmsq_f64(acc_1_j, nuVec, vld1q_f64(&vel[1][j]));
        acc_2_j = vfmsq_f64(acc_2_j, nuVec, vld1q_f64(&vel[2][j]));

        vst1q_f64(&acc[0][j], acc_0_j);
        vst1q_f64(&acc[1][j], acc_1_j);
        vst1q_f64(&acc[2][j], acc_2_j);
	}
  
}

Begoniaish
关注
  • 10
    点赞
  • 踩
  • 5
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
HPC高性能计算可以提升算力的速度，猿代码为此行业不断培养人才
YDM6211的博客
10-09 369
2021年中国信息通信研究院发表的《中国算力发展指数白皮书》显示，通过国家投入产出表模型计算，2020年以计算机为代表的算力产业规模达2万亿元，直接带动经济总产出1.7万亿元，间接带动经济总产出6.3万亿元，即在算力中每投入1元，平均将带动3—4元经济产出。相关机构预测，到2035年5G、互联网、人工智能将为全球经济带来40万亿美元增长，算力正是其背后重要的资源依托和基础设施。
“先进计算”领域仍待开拓，猿代码做先进计算人才的引路人
YDM6211的博客
09-16 569
“先进计算”领域仍待开拓，猿代码做先进计算人才的引路人 “先进计算”是对传统计算机技术发展的一种总结，是面向未来的全新技术。从技术发展的角度来看，先进计算将具有时代特征的多样化的计算技术及方式进行融合，进一步完善传统计算机技术的可靠性、稳定性。从程序架构来看，先进计算拥有更快的速度，更强的计算能力，甚至在融合学习其他计算技术后，随着资源库的发展能够进行自适应性运算。
大环境不确定，高端IT人才独有的大趋势在哪？|猿代码科技
YDM6211的博客
09-06 269
为什么企业招聘程序员越来越看重学历？ 为什么人形机器人被认为是下一个新能源汽车级的风口？ 为什么先进计算即将在中国迎来大发展？
深入理解计算机系统 CSAPP 家庭作业11.10
apple_50569014的博客
07-25 247
【代码】深入理解计算机系统 CSAPP 家庭作业11.10。
打卡第22天------回溯算法
GodLoveMe的博客
07-24 742
开始学习了，希望我可以尽快成功上岸！
【算法/训练】：动态规划（线性DP）
island1314的博客
07-24 884
先创建一个 dp[ n + 2 ][ m + 2 ]，然后让dp[ 0 ][ 1 ] = 1 或者 dp[ 1 ][ 0 ] = 1。注意这样初始化的时候，x需要+1,y也需要+1.dp[i][j] 表示：从[0, 0]出发，到底[i, j]位置，一共有多少种方法。当走到马可以走的地方，dp[ i ][ j ] = 0;对输入的字符矩阵我们按照要求将其。故（i，j）位置的积分应该为dp。对于第一行和第一列我们应该特殊处理，利用。，然后我们就可以按照dp。
代码随想录算法训练营第二十二天|回溯算法part01
qq_43456971的博客
07-25 693
这一题与上题一样，除了终止条件中需要判断和是否为target。多了一个剪枝操作，就是当值已经大于target，直接返回。本题还有一个建立数字和字母之间映射的问题，使用string数组来进行映射。因此回溯算法的宽度是每个集合里面的个数，回溯算法的深度是集合的个数。参数是digits与index，index是用来记录遍历第几个集合。在单个集合1-9中，寻找所有个数为k且和为target的组合。在单个集合1-n中，寻找所有个数为k的组合。这一题是在多个集合中找所有可能的组合。和所有递归一样，都有三部曲。
数据结构与算法--【数组2】力扣练习 || 双指针 / 移除元素 / 数组排序
qq_44578108的博客
07-25 860
简单删除元素 / 指定某元素位置，其余元素位置相对不变双指针指向数组头。快指针扫描数组并判断，慢指针收集变化后的数组元素。不创建新数组，只覆盖原数组。对数组元素排序一个指针指开头，一个指针指结尾。判断两指针指向值大小，循环结束条件first <= last；for循环，while循环均可，本人更习惯while循环。使用malloc为新数组分配空间，注意，创建了新数组。
【C#学习笔记】数组
qq_39150957的博客
07-23 341
本篇文章由浅入深介绍了C#中的数组，并且有代码举例
河南萌新联赛2024第（二）场：南阳理工学院
a13997460860的博客
07-25 1287
这场比赛里面出到了很多的题，都是需要很多的时间去写，主要考察到了代码的实现能力，与codeforces里面的题有些差距。codeforces打多了还是有点不习惯的。题挺好的，尤其是H题感觉很考代码实现能力与STL的使用熟练度。多写写这样的题还是对我们很有好处的。只不过前三题都太水了，没有一点难度。
2024 暑假友谊赛 2
qq_42643660的博客
07-23 834
LCA，排序dp，四维dp，交互题构造，思维构造，思维
大模型算法备案流程最详细说明【流程+附件】
cybtec的博客
07-25 946
大模型备案,模型服务协议 1.大模型备案 2.大模型算法备案 3.大模型网信办备案 4.网信办大模型备案 5.AI大模型备案制度 6.大模型 备案 7.大模型备案要求 8.大模型备案管理制度 9.算法备案和大模型备案区别
蓝桥杯 2024 年第十五届省赛真题 —— 最大异或结点
执梗的博客
07-23 1638
蓝桥杯 2024 年第十五届省赛真题 —— 最大异或结点
C语言100道基础拔高题（1）
2303_78660417的博客
07-23 2842
最后按最小值，中间值，最大值的顺序打印即可。解题思路：首先输出由这几个数字所组成的所有三位数，接着再设置条件，使其输出的三位数不重复，下面我们来看下源代码。值得注意的是：所以题目的代码都是作者自行编写，如有更好的思路或者代码的优化，还请指出更正。解题思路：首先找到能整除这个正整数的最小正数，接着重置循环初值，让它继续查找最小正数，最后将得到的所有正数放到一个数组里，最后设置好格式之后循环打印出来。解题思路：首先确定好区间，接着判断是否是素数，素数就是一个大于1的数，除了1和它本身之外，不能被其他数整除。
2024秋招算法
weixin_43364551的博客
07-22 2100
根据代码随想录和力扣进行算法训练的笔记
算法 —— 暴力枚举
最新发布
fen_0108的博客
07-25 3384
第二看纵向宽发生改变，长重置为1，长为1，宽为2的时候，第一行个数为7个，总共有5 x 7 个，长为2，宽为2的时候 第一行个数为6个，共有5 x 6个……首先是横着的长方形，宽始终为1，长不断发生改变，可以看出长为2的时候，第一行个数为6个，总共有6 x 6个，长为3的时候，总共有6 x 5个……这是一道简单的模拟题，枚举出所有可能情况，不会超过规定时间的，以下附上k
代码随想录算法训练营Day 63| 图论 part03 | 417.太平洋大西洋水流问题、827.最大人工岛、127. 单词接龙
m0_51759998的博客
07-25 454
用两个visited来表示。
2024.7.24 作业
qq_64434282的博客
07-24 2047
【代码】2024.7.24 作业。
DP 整数拆分&不同的二叉搜索树 DAY21
小木公的博客
07-25 915
没完全消化，特别是递推公式…我觉得是,另外还需要取最大值如何理解？//为什么是n + 1?//初始化dp[2] = 1;i <= n;j++){这道题想到用动规的方法来解决，就不太好想，需要举例，画图，分析，才能找到递推的关系。
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值