向量化编程
基本介绍
X86: Intel x86是英特尔公司于1978年推出的16位微处理器；而x86泛指一系列基于Intel 8086且向后兼容的中央处理器指令集架构
Intel ICC和开源的GCC编译器支持SSE/AVX指令的C语言接口（intrinsic，内置函数），在intrin.h头文件中（头文件可能有所不同）
immintrin.h头文件介绍：https://blog.csdn.net/fuxiaoxiaoyue/article/details/83153667
-
函数命名：
第一部分：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编程初识->https://blog.csdn.net/nbu_dahe/article/details/122157205
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 +bc
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
Pytorch初上手：
参考资料：
pytorch上手简要：https://www.elecfans.com/d/883656.html
pytorch中文官方文档：https://pytorch.apachecn.org/2.0/tutorials/beginner/basics/quickstart_tutorial/
