向量化优化入门
该文章主要为观看【张先轶】SIMD向量化编程入门(百校万人计划·高性能计算专才培养公开课系列)以及 AVX指令集_shang_ch的博客-CSDN博客_avx指令集的学习笔记
文章目录
1、什么是向量化优化?
向量化优化借助的是 CPU 的 SIMD 指令,即通过单条指令控制多组数据的运算。它被称为 CPU 指令级别的并行。
区别:标量指令与向量指令
●标量指令,每次处理一 个数据
●又称传统指令和通用指令
●向量指令,每次处理一批数据
●也叫SIMD指令(Single Instruction Multiple Data),单指令,多数据
●很多应用需要这种细粒度的并行性(比如,典型的多媒体应用)
2、为什么要有向量指令?
CPU依靠指令来计算和控制系统,指令集是指CPU能执行的所有指令的集合,每一类CPU都有其支持的指令集。比如说目前intel和AMD的绝大部分处理器都使用X86指令集,因为它们都源自于X86架构。
但无论CPU有多快,X86指令也只能一次处理一个数据,但这种单指令单数据(SISD)的指令集效率并不高,因此,为了提高CPU的工作效率,需要增加一些特殊的指令满足时代进步的需求,这些新增的指令就构成了扩展指令集。
现状,由于功耗限制,单核主频提高陷入瓶颈。需要进一步提升性能。
而提升性能更多依赖体系结构的改进,有以下几种提升性能的方法:
- 提高流水线效率
- 提高Cache效率
- 向量指令,提升向量长度
- 引入多核等等
此文章仅仅针对简明向量指令的入门概念
以Intel为例:
- MMX(64-bit)(意思是一条指令可以处理64比特的数据)
- SSE、SSE2(128-bit)
- AVX、AVX2(256-bit)
- AVX -512(512-bit)
可以发现硬件能力不断翻倍,性能也不断翻倍
当然除了Intel公司有向量指令外,ARM公司也有自己的向量指令集
3、向量化简单示意
eg.数组求和
float vsum (float *v, int n)
{
float s = 0.0; int i;
for (i=0; i<n; i++)
s+=v[i];
return s;
}
假设4向量长度的SIMD指令,其运算历程如下
①
②
③
④
⑤
其机理如下:(Not really SIMD;assumes 4-lanes SIMD)
float vsum (f1oat *v, int n)
{
float vs[4] = {0.0, 0.0, 0.0, 0.0};
float s = 0.0;
int i;
for (i=0; i<n-4; i+= 4) {
vs[0] += v[i ];
vs[1] += v[i+1];
vs[2] += v[i+2];
vs[3] += v[i+3]
//以上四行就是vs[0:3]+=v[i:i+3]
}
s = vs[0] + vs[1] + vs[2] + vs[3];
//★Handle leftover 若碰到n不是4的倍数时
for ( ; i<n; 1++){
s += v[i];
}
return s;
}
4、向量化编译方式
- 自动向量化
利用编译器分析串行程序中控制流和数据流的特征,识别程序中可以向量执行的部分,将标量语句自动转换为相应的SIMD向量语句。 - 手工向量化
通过内嵌汇编码或编译器提供的内函数(Intrinsic函数)来添加SIMD指令
虽然自动向量化方法减轻了程序员的负担,但是所生成
代码的性能仍然有待提高.
5、对数组求和进行向量化编程
采用ARMv8 Neon Intrinsic
#include <arm_ neon.h>
float vsum(float *V,int n){
float32x4_t vs,vv;
float s=0.0;
int i;
vs=vdupq_n_ f32(0.0);
for(i=0; i<n-4; i+=4){
vv = vld1q_f32(&v[i]);
vs += vv;
}
s=vs[0]+vs[1]+vs[2]+vs[3] ;
for(; i<n; i++){
s+=v[i] ;
}
return s;
}
细讲上述代码:
- include头文件arm_neon.h
- 基础数据结构——向量寄存器128-bit(如下图)
float32->float
float64->double
- 向量vs累加中间结果,需要初始vs为0.0,采用vdupq_n_ f32
- 向量vv读取数据,采用vld1q_f32
- 向量加法: ①Intrinsic函数,采用vaddq_f32 ②编译器重载运算符
- 获取向量vs中每个向量的值vs[0]……
- 标量累计s
- 处理剩余尾部
6、Intel的一类SIMD指令集——AVX指令集
AVX指令集是Sandy Bridge和Larrabee架构下的新指令集。AVX是在之前的128位扩展到和256位的单指令多数据流。而Sandy Bridge的单指令多数据流演算单元扩展到256位的同时数据传输也获得了提升,所以从理论上看CPU内核浮点运算性能提升到了2倍。
Intel AVX指令集,在单指令多数据流计算性能增强的同时也沿用了的MMX/SSE指令集。不过和MMX/SSE的不同点在于增强的AVX指令,从指令的格式上就发生了很大的变化。x86(IA-32/Intel 64)架构的基础上增加了prefix(Prefix),所以实现了新的命令,也使更加复杂的指令得以实现,从而提升了x86 CPU的性能。
AVX(Advanced Vector Extensions,高级矢量扩展)指令集借鉴了一些AMD SSE5的设计思路,进行扩展和加强,形成一套新一代的完整SIMD指令集规范。
6.1、查询cpu支持的指令集
①Linux
gcc -march=native -Q --help=target|grep march
或
cat /proc/cpuinfo
②windows
6.2、AVX编程+举例
AVX指令集_shang_ch的博客-CSDN博客_avx指令集已描述较为详细,自行查阅
这里用矩阵相乘的例子进行代码实现:
//#include <bits/types/struct_timeval.h>
#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;
}
运行结果:
虽然但是。。。在测试过程中发现运行时间很不稳定,向量化计算甚至有时比普通计算来得慢。初步推测是数据访问模式不规则(数据在内存中不连续储存,例如上述的__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]);
)或 矩阵维度过小(无法充分利用SIMD指令集的向量化寄存器)造成的。