首页 > 编程语言 >2024.4.7 向量化编程AVX/NEON

2024.4.7 向量化编程AVX/NEON

时间:2024-04-07 16:37:51浏览次数:20  
标签:float32x4 matrix int NEON cols f32 2024.4 AVX struct

基本介绍
X86: Intel x86是英特尔公司于1978年推出的16位微处理器;而x86泛指一系列基于Intel 8086且向后兼容的中央处理器指令集架构

Intel ICC和开源的GCC编译器支持SSE/AVX指令的C语言接口(intrinsic,内置函数),在intrinsic.h头文件中(头文件可能有所不同)

  • 函数命名:
    第一部分:mm/mm256。mm表示SSE指令集,操作长度为64位或128位。mm256表示使用AVX指令集、操作位位256位。
    第二部分:操作函数名称——如 add、load、mul....
    第三部分:操作对象及数据类型— ps表示操作向量数据为单精度、pd表示操作的向量数据为双精度等。

  • 函数命名举例(AVX2):
    _mm256_add_ps:使用AVX 256位寄存器,进行加法操作,操作的向量数据位单精度。
    _mm256_mul_pd:使用AVX 256位寄存器,进行乘法操作,操作的向量数据类型为双精度。

gcc -march=native -Q --help=target|grep march
或
cat /proc/cpuinfo

AVX向量化编程
简单向量加:

#include<bits/stdc++.h>
#include<immintrin.h>
#include<omp.h>
using namespace std;

typedef chrono::high_resolution_clock Clock;

void out(int D[])
{
  for(int i=0;i<8;i++)
  {
    cout<<D[i]<<" ";
  }
  cout<<endl;
}
​
int main()​
{
  __m256i x; // 声明m256变量
  __m256i y;
  __m256i z;
​
  int A[8]={1,2,3,4,5,6,7,8};
  int B[8]={1,2,3,4,5,6,7,8};
  int ans[8];
​
  x = _mm256_loadu_si256((__m256i*)&A[0]); // 载入
  y = _mm256_loadu_si256((__m256i*)&B[0]);
  z = _mm256_add_epi64(x,y);  // 并行加法实现
​
_mm256_storeu_si256((__m256i*)&ans[0],z); // 回存
​
  out(A);
  out(B);
  out(ans);
​
  return 0;
}

以矩阵乘为例介绍AVX编程

#include <bits/stdc++.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <time.h>
#include <immintrin.h>
struct matrix {
	double* data;
	int rows;
	int cols;
};
struct matrix* matrix_create(int rows, int cols)
{
	struct matrix* m = (struct matrix*)malloc(sizeof(struct matrix));
	m->rows = rows;
	m->cols = cols;
	m->data = (double*)malloc(sizeof(double) * rows * cols);
	return m;
}
void matrix_destroy(struct matrix* m)
{
	free(m->data);
	free(m);
}
void matrix_random_init(struct matrix* m)
{
	for (int i = 0; i < m->rows * m->cols; i++)
	{
		m->data[i]=(double)rand();
		//m->data[i] = (double)rand() / RAND_MAX;
	}
}
void random_init()
{
	srand(time(NULL));
	//time(NULL)用来获取当前的系统时间,本质上得到的是一个大整数
	//srand()是一个设置随机数种子的函数,一般在调用rand()之前使用 
}
int check(struct matrix* a, struct matrix* b)
{
	if (a->rows != b->rows || a->cols != b->cols)
	{
		printf("A错误");
		return 0;
	}
	for (int i = 0; i < a->rows * a->cols; i++) 
	{
		if (a->data[i] != b->data[i]) 
		{
			printf("B错误:%d\na[i]=%d,b[i]=%d\n",i,a->data[i] , b->data[i]);
			return 0;
		}
	}
	return 1;
}
struct matrix* matrix_multiple(struct matrix* a, struct matrix* b)
{
	struct matrix* c = matrix_create(a->rows, b->cols);
	struct timeval start, end;
	gettimeofday(&start, NULL);
	for (int i = 0; i < a->rows; i++)
	{
		for (int j = 0; j < b->cols; j++)
		{
			double sum = 0;
			for (int k = 0; k < a->cols; k++) 
			{
				sum += a->data[i * a->cols + k] * b->data[k * b->cols + j];
			}
			c->data[i * c->cols + j] = sum;
		}
	}
	gettimeofday(&end, NULL);
	//gettimeofday是一个系统调用,用于获取当前时间戳。在C语言中,可以使用gettimeofday(&end, NULL)来获取当前时间戳并存储在变量end中。这个函数返回的时间戳是一个Timeval结构体,其中tv_sec表示时间戳的秒数部分,tv_usec表示时间戳的微秒部分。通过将时间戳转换为Timeval结构体,可以方便地比较两个时间戳之间的差异。
	//需要注意的是,gettimeofday()函数返回的时间戳是相对于系统启动时间的,因此可能不是真正的时间戳
	//总之,使用gettimeofday()函数可以方便地获取当前时间戳,但需要注意时间戳的精度和准确性
	printf("before optimize: %ld\n", (end.tv_sec - start.tv_sec) * 1000000 +end.tv_usec - start.tv_usec);
	return c;
}
struct matrix* matrix_multiple_optimize(struct matrix* a, struct matrix* b) 
{
    struct matrix* c = matrix_create(a->rows, b->cols);
    struct timeval start, end;
    
    gettimeofday(&start, NULL);
    for (int i = 0; i < a->rows; i++)
	{
        for (int j = 0; j < b->cols; j++)
		{
            __m256d sum = _mm256_setzero_pd();
            int k=0;
            for (; k < a->cols-a->cols%4; k += 4) {
                __m256d a_vec = _mm256_loadu_pd(&a->data[i * a->cols + k]);
                __m256d b_vec = _mm256_set_pd(b->data[(k+3)*b->cols+j],b->data[(k+2)*b->cols+j],b->data[(k+1)*b->cols+j],b->data[k*b->cols+j]);//注意这里的数组元素地址不连续,且顺序是相反的
                __m256d prod = _mm256_mul_pd(a_vec, b_vec);
                sum = _mm256_add_pd(sum, prod);
            }
            double result[4]={0.0};
            _mm256_storeu_pd(result, sum);
            for(;k<a->cols;k++)
            {
            	result[0]+=a->data[i * a->cols + k] * b->data[k * b->cols + j];
			}
			c->data[i * c->cols + j]=result[0] + result[1] + result[2] + result[3];
        }
    }
    gettimeofday(&end, NULL);
    printf("after optimize: %ld\n", (end.tv_sec - start.tv_sec) * 1000000 + end.tv_usec - start.tv_usec);
    return c;
}
int main() {
	random_init();
	struct matrix* m = matrix_create(100, 2000);
	matrix_random_init(m);
	struct matrix* n = matrix_create(2000, 100);
	matrix_random_init(n);
	struct matrix* r = matrix_multiple(m, n);
	struct matrix* r_optimize = matrix_multiple_optimize(m, n);
	if (check(r, r_optimize)) 
	{
		printf("check ok\n");
	} else 
	{
		printf("check failed\n");
	}
	matrix_destroy(m);
	matrix_destroy(n);
	
	matrix_destroy(r);
	matrix_destroy(r_optimize);
	return 0;
}

编译指令:g++ avx-add.cpp -o avx-add -mavx -mavx2 && ./avx-add
__mm256_loadu_si256:加载数据到寄存器上,起始坐标为数组的开头,即A[0];
__mm256_add_epi64:相加操作,即寄存器对应位置数据进行相加。
__mm256_storeu_si256:寄存器数据写入到结果数组Z内。
ARM Neon技术处理SIMD编程: https://blog.csdn.net/Selenitic_G/article/details/106565566

  • 在移动平台上进行一些复杂算法的开发,一般需要用到指令集来进行加速。NEON 技术是 ARM Cortex™-A 系列处理器的 128 位 SIMD(单指令,多数据)架构扩展,专门针对大规模并行运算设计的,是一个实现ARMv7-A或ARMv7-R架构的系列处理器。ARM NEON技术建立在SIMD的概念上,支持128位向量操作,也称为单指令多数据向量模式。 其本质上使用的是128位NEON SIMD寄存器,这意味着如果操作32位浮点数,可同时操作4个(变量可定义:float32x4_t);如果操作 16 位整数(short),可同时操作 8 个(变量可定义:int16x8_t);而如果操作 8 位整数,则可同时操作 16 个(变量可定义:int8x16_t)。
    ARMv7 NEON 指令集架构具有 16 个 128 位的向量寄存器,命名为 q0~q15。这 16 个寄存器又可以拆分成 32 个 64 位寄存器,命名为 d0~d31。其中qn和d2n,d2n+1是一样的,故使用汇编编写代码时要注意避免产生寄存器覆盖。

float32x4_t vdupq_n_f32 (float32_t value)
将value复制4分存到返回的寄存器中
float32x4_t vld1q_f32 (float32_t const * ptr)
从数组中依次Load4个元素存到寄存器中
float32x4x2_t vld2q_f32 (float32_t const * ptr)
从数组中依次Load8个元素存到两个寄存器中
相应的 有void vst1q_f32 (float32_t * ptr, float32x4_t val)
将寄存器中的值写入数组中
float32x4_t vaddq_f32 (float32x4_t a, float32x4_t b)
返回两个寄存器对应元素之和 r = a+b
float32_t vaddvq_f32 (float32x4_t a)
返回a向量寄存器中元素之和给一个新的标量
相应的 有float32x4_t vsubq_f32 (float32x4_t a, float32x4_t b)
返回两个寄存器对应元素之差 r = a-b
float32x4_t vmulq_f32 (float32x4_t a, float32x4_t b)
返回两个寄存器对应元素之积 r = ab
float32x4_t vmlaq_f32 (float32x4_t a, float32x4_t b, float32x4_t c)
返回乘加值r = a +b
c
float32x2_t vmaxq_f32 (float32x4_t a, float32x4_t b)
返回向量a、b对应通道的最大值给存入到一个新的向量中
float32_t vmaxvq_f32 (float32x4_t a)
求出向量a中的最大值返回给一个新的标量

  • 通过一个示例来解释如何利用NEON内置函数来加速实现统计一个数组内的元素之和。
#include <iostream>
#include <arm_neon.h> //需包含的头文件
using namespace std;
/*
本质上是对寄存器进行处理(128位NEON SIMD寄存器)
*/
float sum_array(float *arr, int len)
{
    if(NULL == arr || len < 1)
    {
        cout<<"input error\n";
        return 0;
    }
 
    int dim4 = len >> 2; // 数组长度除4整数
    int left4 = len & 3; // 数组长度除4余数
    float32x4_t sum_vec = vdupq_n_f32(0.0);//定义用于暂存累加结果的寄存器且初始化为0
    for (; dim4>0; dim4--, arr+=4) //每次同时访问4个数组元素
    {
        float32x4_t data_vec = vld1q_f32(arr); //依次取4个元素存入寄存器vec
        sum_vec = vaddq_f32(sum_vec, data_vec);//ri = ai + bi 计算两组寄存器对应元素之和并存放到相应结果
    }
    //将累加结果寄存器中的所有元素相加得到最终累加值
    float sum = vgetq_lane_f32(sum_vec, 0)+vgetq_lane_f32(sum_vec, 1)+vgetq_lane_f32(sum_vec, 2)+vgetq_lane_f32(sum_vec, 3);
    for (; left4>0; left4--, arr++)
        sum += (*arr) ;   //对于剩下的少于4的数字,依次计算累加即可
    return sum;
}

程序中函数解释:

float32x4_t vdupq_n_f32(float32_t val):将val复制四份放入返回的寄存器中。
float32x4_t vld1q_f32(float32_t const * ptr):从地址ptr依次向后加载四个元素放入返回的寄存器中。
float32x4_t vaddq_f32(float32x4_t a, float32x4_t b):返回a+b的值,向量运算,四个值同时相加。
float32_t vgetq_lane_f32(float32x4_t v, const int lane):返回v中某一个lane的值

向量点乘实现:

void neonax(float A[][1000],float *x,int len)
{
	double t2,t1;
	double y[len];
	float32x4_t vec1,vec2,sum_vec = vdupq_n_f32(0);
	for(int i=0;i<len;i++)
	{
		for(int j=0;j<len;j+=4)
		{
			vec1 = vld1q_f32(A[i]+j); // 寄存器vec1和vec2每次可以取4个float值
			vec2 = vld1q_f32(x+j);
			sum_vec = vmlaq_f32(sum_vec,vec1,vec2); // 使用vmlaq_f32函数进行向量乘
		}
		y[i]=vgetq_lane_f32(sum_vec,0)+vgetq_lane_f32(sum_vec,1)+vgetq_lane_f32(sum_vec,2)+vgetq_lane_f32(sum_vec,3);
		sum_vec=vdupq_n_f32(0);
	}
}

向量加实现:

void add_neon1(float *c, float *a, float *b,int count)
{
	int i;
	float32x4_t in1,in2,out;
	for(i=0;i<count;i += 4)
	{
		in1 = vld1q_f32(a);
		a += 4;
		in2 = vld1q_f32(b);
		b += 4;
		out = vaddq_f32(in1, in2);
		vst1q_f32(c,out);
		c +=4;
	}	
}

自动向量化:

  • 在GCC中启用自动向量化使用命令行选项:
    -ftree-vectorize
    -mfpu=neon
    -mcpu 来指定核心或架构。
    在优化级编译-O3意味着-ftree向量化。如果没有指定-mcpu选项,那么GCC将使用其内置的默认内核。

手动向量化:(主要为使用下图中的intrinsics内联函数方式):

另附参考资料:Neon内置函数:https://blog.csdn.net/billbliss/article/details/78488013
Neon常用函数总结:https://blog.csdn.net/may0324/article/details/72847800
Neon并行化技术简介:https://blog.csdn.net/qq_39589936/article/details/109154895?spm=1001.2101.3001.6650.1&utm_medium=distribute.pc_relevant.none-task-blog-2~default~CTRLIST~Rate-1-109154895-blog-129000993.235^v43^pc_blog_bottom_relevance_base9&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2~default~CTRLIST~Rate-1-109154895-blog-129000993.235^v43^pc_blog_bottom_relevance_base9&utm_relevant_index=2
intel指令集doc官网:https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html
[x86]SIMD指令集发展历程表(MMX、SSE、AVX等):https://www.cnblogs.com/zyl910/archive/2012/02/26/x86_simd_table.html
AVX指令集基本介绍:https://blog.csdn.net/m0_55063425/article/details/128603137

标签:float32x4,matrix,int,NEON,cols,f32,2024.4,AVX,struct
From: https://www.cnblogs.com/jibinghu/p/18119328

相关文章

  • 2024.4.6 组合数学补题
    CF128CGameswithRectangle个人认为突破点是“严格包含”,一开始没注意严格不知道怎么处理。严格的话就是横竖分别在若干条边中,分别选出2k条边。横竖互不影响可以乘法原理,只考虑一个方向即可。#include<iostream>#include<cstdio>#include<algorithm>#definemaxn10000......
  • 2024.4.6练习笔记
    浙江理工大学2024年程序设计竞赛(同步赛)Fleetcode题目要求:求出一个序列中对于每个位置\(i\),以\(i\)为起点第一个\(\text{leetcode}\)子序列的终止位置。需要注意的是不要求子序列连续。不存在则答案为零。容易想到双指针,但是会TLE,考虑一些优化。扫描序列,字母是属于......
  • 2024.4.6 - 4.12
    SatJOI2023Final宣传2\(n\)个人,每个人有住所位置\(X_i\)与影响力\(E_i\),一个人\(i\)拿到书后会号召另一个人\(j\)买书仅当\(|X-i-X_j|\leqE_i-E_j\),你最少送多少个人书才能使得所有人都会有书(送的或者被号召买书)。\(n\leq5\times10^5\)。拆一下绝对值,得:\[......
  • 2024.4.5 RMQ补题
    P2048[NOI2010]超级钢琴前缀和处理连续子段的和弦,然后rmq求最大值运用堆存储最优答案每次取出堆头统计一次后,除掉统计点再分成两段加入即可,共k次#include<bits/stdc++.h>#definemaxn500010#defineINF0x3f3f3f3f#defineintlonglongusingnamespacestd;int......
  • 1.5 警惕和揭秘伪创新(1)《软件方法》2024.4更新
    DDD领域驱动设计批评文集做强化自测题获得“软件方法建模师”称号《软件方法》各章合集1.5警惕和揭秘伪创新初中数学里要学习全等三角形、相似三角形、SSS、SAS……,到了高中以后学了正弦定理、余弦定理等解三角形的知识……就不会再回去用初中的方法解题了。但是,不是所......
  • 2024.4 做题记录
    299.CF1534ELostArray难崩。题意转化为每次翻转\(m\)个\(01\)序列的元素,要把全\(0\)翻成全\(1\)。不想分讨。考虑直接最短路求最小步数,转移就枚举选多少个原本已经有的数。交互就记录方案就行了。300.P9537[YsOI2023]CF1764B很棒的题。考察终态,可以发现最后输......
  • 2024.4 第一周做题记录
    \(2024.4.2\)CF1336DYuiandMahjongSet题意:初始有一个值域在\([1,n]\)的多重整数集\(S\),且每个元素重复次数最多为\(n\),定义\(\operatorname{triple}(S)\)表示\(S\)中相同无序三元组数量,\(\operatorname{straight}(S)\)表示\(S\)中连续整数的无序三元组数量,告诉......
  • #19 2024.4.3
    694.pjudge21633【PER#2】2048695.loj3483「USACO2021.2Platinum」CountingGraphs696.loj2468「2018集训队互测Day2」神秘货币史。697.cf1935fAndrey'sTree反思。考虑一个\(mx\rightarrowmx+1\)的构造。那么它挺赢的。考虑一些cornercase,即\(u=mx......
  • 随堂练习2024.4.3
    建立规则,仪式,流程,模式:  定义代码编写和审查的标准,确保开发质量。  实施敏捷开发的仪式,如日常站会,迭代评审和回顾会议,以提高团队协作和项目透明度。 建立清晰的开发流程和里程碑,确保项目按时推进。给好行为正面的反馈:  对于按时完成任务且代码质量高的开发人员......
  • 2024.4.3每日一题
    mysql1.创建表在里面加备注createtablexxx(idintprimarykeycomment'编号',namevarchar(15)notnullcomment'姓名')2.date和timestamp的区别Date类型只包含日期部分,没有时间部分,一般格式为'YYYY-MM-DD'。Timestamp类型包含日期和时间部分,可以精确到毫秒......