高性能计算 第一次作业

张钧宇
18340211
计算机科学与技术
zhangjunyu@nelson-cheung.cn

运行方法

进入到code/build文件夹下,输入如下命令。

code/build会生成9个可执行文件,名字形如test_x.exex从1到9,分别对应下面每一节的测例,如下所示。

test_x.exe的测试代码是test_x.cpp

对于MKL库,进入文件夹code/mkl,然后输入如下命令。

即可运行使用MKL库的例子。

硬件数据

运行环境:VMware 16.x虚拟机

操作系统:Ubuntu 18.04

内存:8G

CPU:6个逻辑CPU

通用矩阵乘法

对于一个规模为M×K​的矩阵A​和一个规模为K×N的矩阵B,通用矩阵乘法(GEMM)定义为

C=AB

其中,AB的乘积C的每一项为

Ci,j=k=1K(Ai,kBk,j)

首先,为了方便地从代码中表示矩阵和矩阵之间的运算,定义一个矩阵类Matrix

矩阵类的定义中包含了矩阵的大部分操作,如下所示。

通用矩阵乘法的实现如下。

第3-6行,判断是否满足矩阵乘法的条件,如果不满足,则返回空矩阵。

第12-22行,按照通用矩阵乘法公式实现。

实现了矩阵乘法后,可以做一个简单的测试检验实现的正确性。首先,生成两个随机矩阵aba的大小是5×4b的大小是4×3

计算通用矩阵乘法,得到乘积c

运行结果如下。

1-1

可以看到计算的结果是正确的,完整代码保存在code/test/test_1.cpp下。

通用矩阵乘法优化

显然,优化的前提是保证正确性。而下面的优化方法是经过了正确性验证的,但限于篇幅,不再给出验证方法和结果。

计算时间的方法

优化的结果可以通过优化前后的矩阵乘法的运行时间长短来衡量。运行时间包括了两个主要方面,一个是算法的计算时间,另外一个是分配和释放存储中间计算结果的内存所消耗的时间。

由于C函数库的时间计算只能精确到秒级别,并不能满足本次实验的需要。所以,需要借助于Linux的函数库中的gettimeofday函数来得到微秒级别的时间,方法如下。

第6行,利用函数gettimeofday得到微秒级别的时间,获取的时间保存在struct timeval类型的变量tv中。

第7行,struct timeval包含两个域,tv_sectv_usectv_sec存储了当前时间的秒级别的部分,tv_usec存储了当前时间的微秒部分。因此,当前时间的计算公式为

time=tv_sec+tv_usec×106

计算出当前时间后,只需要在计算矩阵乘法前后分别调用Utils::current_time_us,得到的结果之差即为计算时间。

基于算法分析的方法

基于算法分析的方法只实现了strassen算法。

假设矩阵AB都是N×N的方阵,且N2n次幂。即存在整数n,使得

N=2n

对于矩阵乘法

C=AB

首先将A,B,C分成相等的4部分。

A=[A1,1A1,2A2,1A2,2], B=[B1,1B1,2B2,1B2,2], C=[C1,1C1,2C2,1C2,2]

利用分块矩阵乘法可以得到

C1,1=A1,1B1,1+A1,2B2,1C1,2=A1,1B1,2+A1,2B2,2C2,1=A2,1B1,1+A2,2B2,1C2,2=A2,1B1,2+A2,2B2,2

引入中间变量

M1=(A1,1+A2,2)(B1,1+B2,2)M2=(A2,1+A2,2)B1,1M3=A1,1(B1,2B2,2)M4=A2,2(B2,1B1,1)M5=(A1,1+A1,2)B2,2M6=(A2,1A1,1)(B1,1+B1,2)M7=(A1,2A2,2)(B2,1+B2,2)

可以将上述乘法化简为

C1,1=M1+M4M5+M7C1,2=M3+M5C2,1=M2+M4C2,2=M1M2+M3+M6

代码实现如下。

第3-8行,检查是否满足矩阵乘法和方阵的条件,如果不满足,则返回空矩阵。

第10-26行,将矩阵补齐至2的整数幂的形式,多余部分补零。

第18-37行,对矩阵ab分块。

第39-51行,实现strassen算法的计算公式。

在代码中,对矩阵中的块进行处理的函数如下。

实现了strassen算法后,测试strassen算法的运算时间。测试的用例分别是128,256,512,1024,2048的方阵。

结果如下。

2-1

从结果中可以看到,5种不同大小的方阵的运行时间分别是2ms23ms187ms1220ms10448ms​。

完整代码在code/test/test_2.cpp中。

基于软件优化的方法

内存重排

对于通用矩阵乘法公式

Ci,j=k=1K(Ai,kBk,j)

在实际计算的过程当中,需要使用三重循环来分别枚举i,j,k,且i,j,k的枚举顺序并不会影响最终结果的计算。前面的通用矩阵乘法的实现方法是先枚举i,然后枚举j,最后枚举k。这样的实现方法是最为直观的,也是直接反映了矩阵乘法的计算公式。

但是,对于一个使用了行优先存储的CPU来说,之前的实现方法并不是效率最高的。因为对于i,j,k的枚举顺序,对每一个i,j来说,都需要枚举全部的k。此时,对于矩阵B​的访问是按照下面的序列的顺序访问的。

B1,j, B2,j, BK,j

显然,对矩阵B的访问是列优先的,在访问过程中会造成较多的cache失效的情况。为了提高内存访问的效率,对矩阵B的访问也需要遵循行优先的方式。注意到对矩阵A,B,C的访问方式分别是(i,j),(i,k),(k,j),为了体现行优先,只需要ij之前枚举,ik之前枚举,kj之前枚举即可,即遵循i,k,j的枚举方式。

实现代码如下。

Matrix::dot相比,这里只是简单地调整了第14行和第16行的顺序,其余均保持不变。

使用128,256,512,1024,2048的方阵进行测试,结果如下。

2-2

从结果中可以看到,5种不同大小的方阵的运行时间分别是1ms​​,3ms​​,27ms​​,194ms​​,1817ms​​​。

完整代码在code/test/test_3.cpp中。

AVX指令集

SIMD(Single Instruction Multiple Data,单指令多数据流),是一种实现空间上的并行性的技术。这种技术使用一个控制器控制多个处理单元,同时对一组数据中的每一个数据执行相同的操作。在 SIMD 指令执行期间,任意时刻都只有一个进程在运行,即 SIMD 没有并发性,仅仅只是同时进行计算。

在 Intel 的 x86 微架构处理器中,SIMD 指令集有 MMX、SSE、SSE2、SSE3、SSSE3、SSE4.1、SSE4.2、AVX、AVX2、AVX512。

AVX 是 SSE 架构的延伸,将 SSE 的 XMM 128bit 寄存器升级成了 YMM 256bit 寄存器,同时浮点运算命令扩展至 256 位,运算效率提升了一倍。另外,AVX 还添加了三操作数指令,以减少在编码时先复制再运算的动作。

AVX2 将大多数整数运算命令扩展至 256 位,同时支持 FMA(Fused Multiply-Accumulate,融合乘法累加)运算,可以在提高运算效率的同时减少运算时的精度损失。

AVX512 将 AVX 指令进一步扩展至 512 位。

本次实验使用的是AV2指令集。

由于AVX指令集进行的是向量运算,需要将原本的矩阵乘法向量化。为了简述矩阵乘法向量化的思想,这里首先使用4×4的方阵来解释。

对于两个4×4的方阵A,B的矩阵乘法

[a1,1a1,2a2,1a2,2][b1,1b1,2b2,1b2,2]=[c1,1c1,2c2,1c2,2]

对于矩阵C的第一行,计算过程为

c1,1=a1,1b1,1+a1,2b2,1c1,2=a1,1b1,2+a1,2b2,2

可向量化为

[c1,1c1,2]=a1,1[b1,1b1,2]+a1,2[b2,1c2,2]

同理

[c2,1c2,2]=a2,1[b1,1b1,2]+a2,2[b2,1b2,2]

上面的情况不难推广到一般情况。

假设一个AVX向量只能够存储T​个矩阵的元素,对于矩阵乘法

C=AB

矩阵C的第i行第j列开始的、按行遍历的T个元素可以按如下公式计算。

[ci,jci,j+1ci,j+T1]=k=1Kai,k[bk,jck,j+1ck,j+T1]

实现代码如下。

第5-8行,判断是否满足矩阵乘法的条件,如果不满足,则返回空矩阵。

第16-18行,由于这里使用了256位的AVX向量,矩阵的每个元素是float格式,32位。每个AVX向量只能够存下8个矩阵的元素。对于矩阵乘法

C=AB

为了计算方便,需要将矩阵B的行对齐到8的倍数,多余部分补0即可。

第22-35行,AVX矩阵乘法的实现。这里使用了i,k,j​的遍历方式,首先需要将矩阵的元素加载到AVX向量,然后计算,最后写回到结果矩阵res即可。

代码使用了AVX的四个函数。

第37行,由于之前进行了内存对齐,在返回结果之前需要将对齐部分去除。

最后的编译参数需要加上-mfma-mavx

使用128,256,512,1024,2048​的方阵进行测试,结果如下。

2-3

从结果中可以看到,5种不同大小的方阵的运行时间分别是0.5ms​​​​,6.1ms​​​​,58.0ms​​​​,249.2ms​​​​,2511.9ms​​​​​。

完整代码在code/test/test_4.cpp中。

多线程

注意到矩阵乘法的计算过程中

Ci,j=k=1K(Ai,kBk,j)

只有对矩阵C的计算才涉及写操作,对矩阵AB的运算都是读操作,并且矩阵C的各个元素的计算是相互独立的。因此,矩阵C的计算不会产生race condition的现象,可以容易地进行并行化处理。

可以使用OpenMP来进行并行化处理,实现代码如下。

实现代码只是在原来的内存重排的优化方法的基础上加上了一条语句即可。

最后的编译参数需要加上-fopenmp

使用128,256,512,1024,2048​的方阵进行测试,结果如下。

2-4

从结果中可以看到,5种不同大小的方阵的运行时间分别是0.4ms​​​​,3.1ms​​​​,20.2ms​​​​,99.1ms​​​​,424.5ms​​​​​。

完整代码在code/test/test_5.cpp中。

编译器优化

g++编译选项中有4个编译优化选项。

分别使用上面4个优化选项测试在通用矩阵乘法上的性能

-O0的性能如下。

2-5-1

从结果中可以看到,5种不同大小的方阵的运行时间分别是10ms​,127ms​,1056ms​,77323ms​,131530ms​。

 

-O1的性能如下。

2-5-1

从结果中可以看到,5种不同大小的方阵的运行时间分别是8ms​,70ms​,592ms​,4211ms​,100477ms​。

 

-O2的性能如下。

2-5-1

从结果中可以看到,5种不同大小的方阵的运行时间分别是2ms28ms271ms2160ms100289ms

 

-O3的性能如下。

2-5-1

从结果中可以看到,5种不同大小的方阵的运行时间分别是3ms22ms314s2109ms100259ms

综合方法

可以将多线程、内存重排、编译器优化综合起来优化,即在strassen算法中的矩阵乘法的实现使用了多线程+内存重排的方式,最后的代码编译使用-O3优化。

实现代码如下。

使用128,256,512,1024,2048​的方阵进行测试,结果如下。

2-3

从结果中可以看到,5种不同大小的方阵的运行时间分别是1ms​​​​,25ms​​​​,46ms​​​​,100ms​​​​,422ms​​​​​。

完整代码在code/test/test_6.cpp中。

MKL

安装了Intel的oneAPI Base Toolkit后,使用oneapi-cli命令创建一个MKL矩阵乘法的例子,无需自己编写。

使用128,256,512,1024,2048​的方阵进行测试,结果如下。

2-3

从结果中可以看到,5种不同大小的方阵的运行时间分别是0.09ms​​​​​,0.07ms​​​​​,0.11ms​​​​​,0.10ms​​​​​,0.09ms​​​​​​。

MKL的代码保存在code/mkl下。

综合比较

综合比较结果如下。

Figure_2

可以看到MKL库的加速比远超于其他的所有优化方法,达到了105106的数量级,原因可能是MKL库可能利用了包括上面的所有优化方法,且利用了许多CPU特性和矩阵本身的特征。

去掉MKL项后的比较结果如下。

Figure_2

从结果中可以看到,每一种加速方法均有明显的加速效果。

从平均加速效果来看,利用多线程加速的效果最好,内存重排的效果次之,AVX的效果稍逊,而利用strassen算法的效果稍差,O3优化的效果最差。可能的原因是AVX和strassen算法中都涉及到许多的内存复制操作,消耗了部分时间。内存重排并不会涉及到内存复制操作,因此加速效果比前者显著。多线程本身使用了内存重排又利用了CPU的6个核心,因此效果最好。代码默认使用O3优化,因此前面的优化效果是叠加了O3优化的效果的,相比之下O3优化的效果最差。

对于综合方法,strassen算法涉及的内存复制操作削减了综合方法内部多线程的加速效果,因此当矩阵规模较小时,加速效果接近于单纯的strassen算法。但当矩阵规模增大时,矩阵的计算成为性能的主要影响因素后,strassen算法的作用便可以发挥出来,使得综合方法的加速效果优于单纯的多线程加速。

大规模矩阵计算优化

在给出大规模矩阵算法与实现之前,需要作出如下假设。

  1. 矩阵很大,内存中不仅无法保存任何一个矩阵,还无法保存矩阵的任意一行或任意一列。

  2. 外存很大,外存可以保存任意多个矩阵。

  3. 矩阵的行和列数分别使用一个4字节的整数表示。

  4. 矩阵的每个元素使用一个4字节的单精度浮点数表示。

  5. 矩阵保存在二进制文件中,保存格式为先存储行数和列数,然后按行优先的方式存储矩阵的元素,即。

     

大规模通用矩阵乘法

基于假设1,假设内存中只可以存储下两个行向量和少量的变量。每个行向量可以保存ROW_MAX个矩阵的元素。基于假设5,采用行优先的方式读取矩阵可以加快矩阵的读取速度,并且可以一次读入ROW_MAX。此时不难想到,大规模通用矩阵乘法应该类似于AVX的向量化计算方法。在AVX的计算方法中,每次计算前需要将行向量从内存中载入AVX向量。类似地,只需要在每次计算前将行向量从外存中载入内存中,然后使用和AVX计算相同的计算过程即可实现大规模通用矩阵乘法。

实现代码如下。

假设矩阵乘法是

C=AB

矩阵A,B,C的保存路径分别是path_apath_bpath_c

第3-18行,打开文件进行读写。

第20-33行,读入矩阵A,B的行数和列数,判断是否满足矩阵乘法的条件。

第35-39行,计算矩阵C的行数和列数并写入结果文件。

第43-49行,初始化矩阵C​的结果文件为零,写入的数量和矩阵C​的大小相同,为后面矩阵计算时随机访问矩阵C的结果文件​保留空间。

第51-53行,向内存申请两个大小为ROW_MAX的向量。基于前面的假设,内存中只能存储下两个固定大小的行向量和若干变量。此后,不会向内存中申请任何内存。

第65-88行,执行向量化计算,计算公式如下。

[ci,jci,j+1ci,j+T1]=k=1Kai,k[bk,jck,j+1ck,j+T1]

第69行,读入矩阵A​的元素ai,k​。

第73-74行,计算矩阵B,C​的读取位置,读取行向量的函数如下。

第75行,由于矩阵B的列数并不一定是8的倍数,需要做边缘处理。

第77-78行,读入B,C的行向量。

第80-83行,执行向量化计算。

第85行,将计算结果写回矩阵C​的对应部分,写入行向量的函数如下。

为了测试上述大规模矩阵计算方法的正确性,需要将ROW_MAX设置为一个较小的值。

同时,矩阵乘法可以通过之前的一般方法计算。通过比较一般方法和大规模计算方法的结果即可验证结果的正确性。

完整代码保存在code/src/huge_matrix.cpp中。

分别使用三个测例进行测试,下面的矩阵A,B都是随机矩阵。

然后比较大规模矩阵乘法的结果和一般方法的结果,结果如下。

3-1

可以看到三个测例中,大规模矩阵乘法的结果和一般方法的结果差异为0,均通过正确性验证。由于算法并不依赖于具体的ROW_MAX,在实际过程中将ROW_MAX调大后,即可推广到一般情况。

测试代码如下。