AVX 指令集并行技术优化积分计算圆周率 π

通过 AVX 指令集并行技术优化积分计算圆周率 π

完整代码和解释如下

  

   

    

     

    
    

     

      // AVX_PI.cpp : 定义控制台应用程序的入口点。
     
    
   

    

     

    
    

     

      //
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      #include "stdafx.h"
     
    
   

    

     

    
    

     

      #include <iostream>
     
    
   

    

     

    
    

     

      #include <immintrin.h>
     
    
   

    

     

    
    

     

      #include <time.h>
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      #include "timer.h"
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      double compute_pi_naive(size_t dt)
     
    
   

    

     

    
    

     

      {
     
    
   

    

     

    
    

     
	double pi = 0.0;
     
    
   

    

     

    
    

     
	double delta = 1.0 / dt;
     
    
   

    

     

    
    

     
	for (size_t i = 0; i < dt; i++)
     
    
   

    

     

    
    

     

      	{
     
    
   

    

     

    
    

     
		double x = (double)i / dt;
     
    
   

    

     

    
    

     

      		pi += delta / (1.0 + x*x);
     
    
   

    

     

    
    

     

      	}
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
	return pi*4.0;
     
    
   

    

     

    
    

     

      }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      double compute_pi_avx(size_t dt)
     
    
   

    

     

    
    

     

      {
     
    
   

    

     

    
    

     
	double pi = 0.0;
     
    
   

    

     

    
    

     
	double delta = 1.0 / dt;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      	__m256d ymm0, ymm1, ymm2, ymm3, ymm4;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      	ymm0 = _mm256_set1_pd(1.0);	// 注意并不是对应一条指令，而是多条指令混合而成
     
    
   

    

     

    
    

     

      	ymm1 = _mm256_set1_pd(delta);
     
    
   

    

     

    
    

     

      	ymm2 = _mm256_set_pd(delta * 3, delta * 2, delta, 0.0);	// 0.0， delta, 2*delta, 3*delta
     
    
   

    

     

    
    

     

      	ymm4 = _mm256_setzero_pd();
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
	for (int i = 0; i <= dt - 4; i += 4)
     
    
   

    

     

    
    

     

      	{
     
    
   

    

     

    
    

     

      		ymm3 = _mm256_set1_pd(i * delta);	// 归一化构造 x（积分区间为 0-1）
     
    
   

    

     

    
    

     

      		ymm3 = _mm256_add_pd(ymm3, ymm2);	// 分别累加四项（分别相差一个delta）
     
    
   

    

     

    
    

     

      		ymm3 = _mm256_mul_pd(ymm3, ymm3);	// 平方
     
    
   

    

     

    
    

     

      		ymm3 = _mm256_add_pd(ymm0, ymm3);	// + 1
     
    
   

    

     

    
    

     

      		ymm3 = _mm256_div_pd(ymm1, ymm3);   // delta 向量除
     
    
   

    

     

    
    

     

      		ymm4 = _mm256_add_pd(ymm4, ymm3);	// 累加
     
    
   

    

     

    
    

     

      	}
     
    
   

    

     

    
    

     

      #ifndef _WIN32
     
    
   

    

     

    
    

     
	double tmp[4] __attribute__((aligned(32)));
     
    
   

    

     

    
    

     

      #else
     
    
   

    

     

    
    

     

      //#pragma pack(8)
     
    
   

    

     

    
    

     

      	__declspec(align(32))	// arm 下的内存边界对齐使用 __align
     
    
   

    

     

    
    

     
	double tmp[4];
     
    
   

    

     

    
    

     

      //#pragma pack()
     
    
   

    

     

    
    

     

      #endif
     
    
   

    

     

    
    

     

      	_mm256_store_pd(tmp, ymm4);
     
    
   

    

     

    
    

     

      	pi = tmp[0] + tmp[1] + tmp[2] + tmp[3];
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
	return pi * 4.0;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      int _tmain(int argc, _TCHAR* argv[])
     
    
   

    

     

    
    

     

      {
     
    
   

    

     

    
    

     
	size_t dt = 1e10;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      	Timer timer;
     
    
   

    

     

    
    

     
	double pi_naive = compute_pi_naive(dt);
     
    
   

    

     

    
    

     
	printf("Pi = %lf, %lfms\n", pi_naive, timer.elapsed());
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      	timer.restart();
     
    
   

    

     

    
    

     
	double pi_avx = compute_pi_avx(dt);
     
    
   

    

     

    
    

     
	printf("Pi = %lf, %lfms\n", pi_avx, timer.elapsed());
     
    
   

    

     

    
    

     
	
     
    
   

    

     

    
    

     
	return 0;
     
    
   

    

     

    
    

     

      }
     
    
   

    

     

    
    

     
 
     
    
  
 

最后输出的结果，获得了 3 倍多的加速（耗时不到原来串行版本的 1/4 ）

文章来源: panda1234lee.blog.csdn.net，作者：panda1234lee，版权归原作者所有，如需转载，请联系作者。

原文链接：panda1234lee.blog.csdn.net/article/details/85016597