贾 迅 邬贵明 谢向辉 吴 东

(数学工程与先进计算国家重点实验室 江苏无锡 214125)

矩阵乘运算核心广泛应用于数值计算、数字信号分析、图像处理、人工智能等领域，其计算形式为

[i,j]+=A[i,k]×B[k,j],

其中，A∈M×K,B∈K×N,C∈M×N.循环k的计算存在写后读数据相关，无法并行；循环i和循环j的计算不存在真相关，可以并行.另外，计算过程中的存储访问均比较规整.由于计算复杂度为O(n3)，双精度浮点矩阵乘往往是大规模数值计算应用的性能瓶颈[1]，其性能直接反映了系统的计算能力.以性能测试程序Linpack[2]为例，程序实现了基于双精度浮点矩阵乘的线性方程求解，并以单位时间内系统执行浮点操作的次数(FLOPS)来衡量系统的计算能力.因此，为双精度浮点矩阵乘运算提供高性能且高效的计算平台对大规模数值计算应用的性能发挥至关重要.

目前，矩阵乘的计算平台主要包括CPU，GPGPU，FPGA.其中，CPU和GPGPU平台上矩阵乘的实现和应用较为成熟.处理器厂商会提供针对各自产品架构进行性能优化的计算函数接口库，如Intel MKL[3]，IBM ESSL[4]，NVIDIA cuBLAS[5]等.由于处理器架构面向通用计算，CPU和GPGPU上矩阵乘的计算效率无法达到最优，如Intel Knights Corner众核处理器矩阵乘的计算效率仅为89%[6]，NVIDIA P100的计算效率约为93%[7].矩阵乘效率对系统计算效率的影响几乎是线性的，矩阵乘效率的损失必然导致系统计算效率的降低[8].

FPGA计算平台具有灵活的可编程性、细粒度的并行能力和丰富的逻辑资源.同时，矩阵乘运算数据并行性好和存储访问规整的算法特点适合并行计算结构的设计[9]，因此基于FPGA平台实现针对矩阵乘的并行计算结构成为学术研究的热点，并且已经取得一些研究成果[10].以浮点乘加运算部件作为计算单元(processing elements， PE)的核心，多个计算单元组成广播[11]或阵列结构[12-14]进行加速计算是典型的实现方式.其中文献[12]提出的线性阵列计算结构可以处理任意规模的矩阵乘运算，且数据传输和存储的效率较高.由于矩阵乘的性能受限于FPGA的可用的逻辑资源和时钟频率，即使以资源利用和时钟频率为优化目标的结构设计，其计算性能也只能达到203.1 GFLOPS[15](252个PE，时钟频率为403 MHz).显然，目前基于FPGA的矩阵乘无法满足大规模数值计算应用对计算能力的需求.

近年来，面向特定应用的硬件定制加速成为高性能、高效计算结构设计的一个重要趋势[16-17].Google公开了TPU(tensor processing unit)芯片的设计[18]，其核心是整数矩阵乘脉动阵列计算结构.NVIDIA最新发布的V100 GPGPU[19]芯片中引入了称为“Tensor Core”的16 b，32 b混合精度矩阵乘阵列计算结构，单芯片集成了672个TensorCore，可提供120 TFLOPS的计算能力.由于TPU和V100芯片中矩阵乘定制加速结构的设计面向深度学习类应用，结构支持的计算精度较低[20-21]，因而无法直接应用于对计算精度要求较高的数值计算领域.

本文以线性阵列计算结构为基础，对其进行结构优化，并采用硬件定制的方法设计实现了支持双精度浮点运算的矩阵乘协处理器.同时，本文建立了性能模型并深入分析了各结构设计参数对协处理器实际计算性能和效率的影响.此外，本文还验证了矩阵乘协处理器的功能正确性并评估了其硬件实现的开销.本文探索了硬件定制结构设计在双精度浮点矩阵乘加速计算中的应用，研究成果对提升现有计算系统的性能和效率有一定的借鉴意义.

1 线性阵列计算结构及大规模矩阵乘算法

线性阵列计算结构[12]如图1所示，多个计算单元线性互连，每个计算单元包含局部存储器c以存储矩阵C的分块数据，寄存器a和b分别存储矩阵A和矩阵B的数据元素，一个双精度浮点乘加运算部件和一个计数器CNT.

Fig. 1 Linear array computing architecture图1 线性阵列计算结构

数据b在PE之间进行传输，其有效时驱动PE执行一次乘加计算：c[CNT]=c[CNT]+a×b.

算法1给出了基于上述线性阵列计算结构的大规模分块矩阵乘算法[12].参数Sp，St分别是矩阵C行数M和列数N的分块粒度，同时也表示线性阵列计算结构中PE的数目和各PE集成的局部存储器c的深度.算法中索引为p和t的最内层循环对应线性阵列计算结构执行的浮点矩乘迭代计算，即Sp个PE实现矩阵A中规模为Sp×1的矩阵子块与矩阵B中规模为1×St的矩阵子块乘计算，并将计算结果与矩阵C中规模为Sp×St的子块相加.阵列计算结构在单次迭代计算中共执行2SpSt次双精度浮点运算.

Fig. 2 Data block mapping for the linear array图2 矩阵分块数据到阵列计算结构的映射

算法1. 基于线性阵列结构的大规模矩阵乘算法.

输入：双精度浮点矩阵A、矩阵B、矩阵C；

输出：双精度浮点矩阵C.

① ForTp=0 toM-1,Sp

② ForTt=0 toN-1,St

③ LoadC[Tp:Tp+Sp-1,Tt:Tt+St-1]；

④ Fork=0 toK-1

⑤ LoadA[Tp:Tp+Sp-1,k]；

⑥ LoadB[k,Tt:Tt+St-1]；

⑦ Forp=TptoTp+Sp-1

⑧ Fort=TttoTt+St-1

⑨C[p,t]+=A[p,k]×B[k,t]；

⑩ End For

算法1中索引为Tp，Tt，k的3层循环在线性阵列计算结构进行矩阵乘迭代计算的同时，处理矩阵分块数据的加载与写回，其功能在硬件上由数据传输控制部件DTC(data transfer controller)实现.各矩阵分块数据到线性阵列计算结构的映射如图2所示：

线性阵列结构经过K轮迭代计算才能得到矩阵C中一个分块数据的计算结果.因此，在K轮迭代计算开始前，数据传输控制部件需要先将矩阵C中规模为Sp×St的分块数据C[Tp:Tp+Sp-1,Tt:Tt+St-1]分别加载到Sp个计算单元的深度为St的局部存储器c中，并在K轮迭代计算完成后将对应的分块计算结果写回主存.另外，每次迭代计算中各PE需要访问矩阵A中相同的数据元素.因此，每次迭代计算开始前，数据传输控制部件还需要将矩阵A中规模为Sp×1的分块数据A[Tp:Tp+Sp-1,k]加载到Sp个计算单元的寄存器a中.

矩阵A和矩阵C的分块数据全部加载完成后，数据传输控制部件再加载矩阵B的数据元素以驱动阵列计算结构进行迭代计算.矩阵乘运算执行完成共需K×(MSp)×(NSt)次迭代计算.对计算过程中数据传输控制部件传输矩阵A，B，C分块数据的数据量、时间间隔以及次数的量化分析如表1所示：

Table 1 Analysis on Data Transfer for Matrix Multiplication表1 针对矩阵乘运算过程中数据传输的分析

文献[12]基于Xilinx公司Virtex-5系列的FPGA实现了上述线性阵列计算结构.PE个数Sp设置为96、局部存储器深度St设置为1 024的情况下，FPGA的最大运行频率为158.08 MHz.输入矩阵的规模M，N，K均为8 192时，阵列计算结构完成双精度浮点矩阵乘运算所需的时间为40.48 s，即结构的计算性能为27.16 GFLOPS，效率为89.49%.

线性阵列计算结构虽然存储效率较高，但每次迭代计算中各PE需要等待计算所需的矩阵分块数据全部加载完成后才能开始计算，计算与存储访问在时间上是串行的，计算资源的闲置必然导致计算效率的损失，高效计算结构的设计必须对此进行改进.由于结构设计面向FPGA计算平台，文献[12]仅以FPGA可用的逻辑资源作为设计约束，并未考虑局部存储器的深度、输入矩阵的规模等参数对线性阵列计算结构实际计算性能和效率的影响，而这些参数在定制协处理器的设计中尤为重要.另外，本文协处理器的设计除了考虑核心的计算结构外，还需要考虑其与主机端的接口，以实现控制信息和计算数据的传输，从而构成完整的加速计算系统.

2 协处理器的设计

2.1 线性阵列计算结构的优化

根据算法1给出的大规模分块矩阵乘运算算法，线性阵列计算结构在进行矩阵乘迭代计算前，需要先从主存读入矩阵A和C的分块数据；K轮迭代计算完成后，需要先将矩阵C的分块计算结果写回主存，才能开始下一轮的迭代计算.矩阵分块数据在主存和线性阵列计算结构之间进行传输的过程中，结构中各计算单元无法进行计算，计算资源被闲置，从而导致线性阵列计算结构整体的效率较低.针对这个问题，本文采用双缓冲技术对阵列计算结构中各计算单元的设计进行了优化.优化后，计算单元的结构如图3所示：

Fig. 3 Structure of PE optimized with double buffering图3 采用双缓冲优化技术的PE结构

图3所示的PE结构中除了已有的寄存器a，b和局部存储器c之外，还包括双缓冲优化引入的存储矩阵A，B数据元素的寄存器a_buf，b_buf以及存储矩阵C分块数据的局部存储器c_buf.

应用双缓冲优化技术后，阵列计算结构可以在使用当前缓冲中的数据进行迭代计算的同时，提前将下一轮迭代计算所需的矩阵分块数据装入到另一缓冲，实现计算与访存的重叠，隐藏存储访问的开销.此时，高效的存储访问调度成为应用双缓冲优化提升计算结构资源利用率的关键.

2.2 存储访问调度

由表1列出的矩阵乘运算过程中矩阵分块数据的传输可知，线性阵列计算结构在进行矩阵乘运算时，每个时钟周期仅需访问矩阵B的一个数据元素，其对存储带宽的需求较低；而矩阵A，C的访问时间间隔分别为St和K×St个时钟周期，因此可以利用读取矩阵B数据元素剩余的存储访问带宽提前加载矩阵A，C的分块数据至PE的缓冲中，或将矩阵C的分块计算结果从缓冲写回主存.

算法2给出了基于双缓冲的存储访问调度算法.数据传输控制部件在硬件上实现了该算法，并支持双缓冲读写访问的控制.

算法2. 基于双缓冲的存储访问调度算法.

输入：双精度浮点矩阵A、矩阵B、矩阵C；

输出：双精度浮点矩阵C.

① LoadC[0:Sp-1,0:St-1]；

② LoadA[0:Sp-1,0]；

③ ForTp=0 toM-1,Sp

④ ForTt=0 toN-1,St

⑤ Fork=0 toK-1

⑥ LoadB[k,Tt:Tt+St-1]；

LoadA[Tp:Tp+Sp-1,k+1] or

StoreC[Tp:Tp+Sp-1,Tt-St:Tt-1] or

LoadC[Tp:Tp+Sp-1,Tt+St-1:Tt+2St-1]；

⑦ Forp=TptoTp+Sp-1

⑧ Fort=TttoTt+St-1

⑨C[p,t]+=A[p,k]×B[k,t]；

⑩ End For

算法2中，阵列计算结构进行第一次迭代计算所需的矩阵A，C的分块数据先被读入.读取矩阵B的数据元素进行迭代计算时，优先利用剩余的存储访问带宽提前读取下一次迭代计算所需的矩阵A的分块数据；待矩阵A的分块数据读入完成，再利用剩余带宽写回前K轮迭代计算得到的矩阵C的分块计算结果；待计算结果写回完成，最后利用剩余带宽装入后K轮迭代计算所需的矩阵C的分块数据.所有迭代计算完成后，矩阵C中最后一个分块计算结果被写回.

理想情况下，线性阵列计算结构在开始新一次的迭代计算时，各PE所需的矩阵数据均已在缓冲中，即存储访问的延迟可以有效地被计算隐藏，结构中的计算资源可以得到充分利用.而实际的存储访问调度受限于存储带宽、计算单元数、局部存储器的深度以及矩阵规模等参数，其对矩阵乘协处理器实际性能的影响将在2.3节进行分析.

2.3 系统结构

基于矩阵乘协处理器的加速计算系统的整体结构如图4所示.系统由主机端处理器芯片和双精度浮点矩阵乘协处理器芯片2部分组成.其中，矩阵乘协处理器通过IO总线与主机端处理器相连.系统采用PCIe-3.0×16作为芯片间的互连接口，双向通信带宽最高为31.5 GBps.

Fig. 4 System architecture based on matrix multiplication coprocessor图4 基于矩阵乘协处理器计算系统的架构

主机端和协处理器端存储空间相互独立，待计算的矩阵数据存储在协处理器端的内存.本文协处理器的设计面向大规模矩阵乘，支持的矩阵规模至少为16 384×16 384.因此，主机端与协处理器端的内存容量配置为8 GB，并采用DDR3-2133[22]规格的存储芯片，接口的数据宽度为128 b，带宽为34.13 GBps.

系统对主机端处理器的类型没有限制，只要支持PCIe传输协议即可，其指令集系统架构可以是x86，Power，ARM，RISC-V[23]等.矩阵乘协处理器由应用双缓冲优化的线性阵列计算结构、存储访问控制器、协处理器与主机端接口(coprocessor-host interface， CHI)3部分构成.乘加运算部件是计算单元的核心，本文采用全流水且兼容IEEE-754浮点数标准[24]的乘加运算部件设计，以保证协处理器的计算性能.CHI接口实现主机端与协处理器之间控制和数据的传输.主机端向矩阵乘协处理器发送的主要命令如下：

1)Initialization对协处理器的控制与状态寄存器进行初始化；

2)DataIn将待计算矩阵数据从主机端内存加载到协处理器端内存；

3)MatrixMultiply启动协处理器上的双精度浮点矩阵乘运算；

4)DataOut将矩阵乘运算的结果从协处理器端内存读出到主机端内存.

接收到来自主机端发送的命令，CHI接口对命令进行分析并生成线性阵列计算结构相关的控制信号.协处理器在进行矩阵乘运算的过程中，主机端与协处理器之间不需要进行数据传输，从而简化了协处理器的设计与调试.

应用程序在计算系统上加速矩阵乘运算的过程如图5所示.当矩阵乘协处理器连接至主机端时，主机端发送Initialization命令对协处理器进行初始化，如设置浮点计算的舍入模式、溢出控制等.应用程序运行在主机端，当其调用双精度浮点矩阵乘运算核心，如基础线性代数子程序库(basic linear algebra subroutines， BLAS)中的函数blas_dgemm(M,N,K,A,B,C,…)时，主机端通过CHI接口向协处理器发送DataIn命令和函数调用传入的矩阵规模、矩阵数据在内存中的起始地址等参数.协处理器的数据传输控制部件将待计算的矩阵数据从主机端内存以直接内存访问(direct memory access， DMA)方式传输至协处理器端内存.数据传输完成后，主机端发送MatrixMultiply命令以启动协处理器上的矩阵乘运算.运算过程中，数据传输控制部件按算法2设计的存储访问调度实现矩阵分块数据在协处理器端内存和线性阵列计算结构之间的传输.

Fig. 5 Computation based on matrix multiplication coprocessor图5 基于矩阵乘协处理器的计算过程

协处理器完成矩阵乘运算后，CHI接口向主机端发出中断信号.主机端收到中断后，向CHI接口发送DataOut命令，协处理器的数据传输控制部件以DMA的方式将矩阵乘运算的结果从协处理器端内存传输回主机端内存.运算结果传输完成后，矩阵乘运算调用结束，应用程序继续执行.

3 性能模型

为了量化分析双精度浮点矩阵乘协处理器的各结构设计参数对其实际计算性能和效率的影响，本文基于算法2设计的存储访问调度算法构建了协处理器的性能模型.模型中的基本参数包括线性阵列计算结构中PE的个数Sp、PE集成的局部存储器的深度St以及输入矩阵的规模M，N，K.此外，本文将模型中协处理器的时钟频率记为Freq(单位为GHz)，协处理器端内存的访问带宽记为Band(单位为GBps)，矩阵乘运算核心在协处理器上执行完成所需的时间记为T(单位为s)，协处理器的实际计算性能和效率记为Perf(单位为GFLOPS)和Eff.

当计算过程中的存储访问延迟完全被计算隐藏，协处理器的计算性能可以充分发挥，此时的计算效率最高.就本文提出的存储访问调度算法，每次迭代计算提前读入矩阵A中规模为Sp的分块数据的延迟与读入矩阵B中规模为St的分块数据的延迟之和应当小于等于对应的计算延迟，即St个时钟周期.同理，K轮迭代计算的总存储访问延迟应当小于等于对应的总计算延迟，即K×St个时钟周期.存储访问的数据粒度为8 B的双精度浮点数，上述2个约束条件可表示为

(1)

协处理器端DDR3-2133存储芯片的数据传输率为2.133 GTs，单次传输的数据宽度为128 b.若协处理器的工作时钟频率Freq=1.5 GHz，则其每个时钟周期可访问2个双精度浮点数.此时，式(1)的约束条件可进一步改写为

(2)

若参数Sp，St，K可以满足式(2)中的约束条件，则计算过程中的存储访问延迟可以完全被计算隐藏.此时，矩阵乘运算核心在协处理器上的执行完成所需的时间T最小.其由3部分时间构成，包括读入矩阵A和C第1个子块数据的时间、迭代计算的时间以及将矩阵C中最后1个子块的计算结果写回内存的时间.计算时间T可表示为

(3)

若参数Sp，St，K不满足式(2)的约束条件，存在St1-SpSt这2种可能.下文将分别分析这2种情况下矩阵乘运算核心在协处理器上的执行时间.

对于第1种情况，提前读入矩阵A子块数据的存储访问延迟无法完全被计算时间隐藏，每次迭代计算的时间增加(Sp-St)2个时钟周期.根据算法2，当前迭代计算执行过程中，写回前K轮迭代计算中矩阵C子块计算结果和提前读入后K轮迭代计算中矩阵C子块计算数据的存储访问延迟均无法被计算隐藏，每K轮迭代计算的时间增加Sp×St个时钟周期.此时，矩阵乘运算核心在协处理器上执行完成所需的时间T可表示为

(4)

对于第2种情况，当前迭代计算的执行时间可以完全隐藏提前读入矩阵A子块数据的存储访问延迟，但只能部分隐藏写回前K轮迭代计算中矩阵C子块计算结果和提前读入后K轮迭代计算中矩阵C子块计算数据的存储访问延迟，其导致每K轮迭代计算的时间增加Sp×St-K×(St-Sp)2个时钟周期.此时，矩阵乘运算核心在协处理器上执行完成所需的时间T可表示为

(5)

协处理器完成矩阵乘核心的加速计算共需执行2MNK次双精度浮点操作.基于对3种情况下矩阵乘运算核心在协处理器上执行时间T的分析结果，可以构建式(6)所示的协处理器的性能模型.此时，协处理器的计算效率可表示为实际计算性能与峰值计算性能的比值，即Eff=Perf(2×Sp×Freq).

4 实验与结果

本节以第3节提出的性能模型为基础，分析矩阵乘协处理器结构设计中双缓冲优化技术、局部存储器的深度和输入矩阵规模等参数对协处理器实际计算性能和效率的影响.实验假设线性阵列计算结构的PE数为1 024，协处理器的时钟频率为1.5 GHz，待计算的矩阵均为方阵，即M=N=K=n.此外，本节还验证了矩阵乘协处理器功能实现的正确性；在主流工艺下评估了其硬件实现的开销，并与同等工艺的GPGPU芯片进行了对比，以探索其硬件可实现性和在实际应用中的优势.

4.1 双缓冲优化技术对矩阵乘协处理器性能的影响

实验在不同的矩阵规模下，分别统计了应用双缓冲优化技术前后矩阵乘协处理器的性能，以分析应用双缓冲优化技术对协处理器性能的影响.实验中，各PE局部存储器的深度均设置为2 048，实验的结果如图6所示.

从图6可以看出，不同矩阵规模下，双缓冲优化技术的应用均可以提升矩阵乘协处理器的计算性能，协处理器的性能平均提升30.06%.矩阵规模为4 096时，性能提升最为明显，达到了45.45%.可见，本文设计的存储访问调度算法可以在计算过程中利用双缓冲优化有效隐藏存储访问的延迟.

同时，双缓冲优化对矩阵乘协处理器的性能提升受矩阵规模的影响.对此分析如下：当矩阵规模较小时，计算无法有效隐藏存储访问的延迟，应用双缓冲优化对性能的提升作用有限；随着矩阵规模逐渐增加，计算可以隐藏更多的访存延迟，此时应用双缓冲优化可以显著提升协处理器的计算性能；当矩阵规模较大时，计算可以很好隐藏存储访问的延迟，而实验中局部存储器的深度保持不变，双缓冲优化对性能的提升作用逐渐降低.

Fig. 6 Effect of double-buffering optimization technique on coprocessor performance图6 双缓冲优化技术对矩阵乘协处理器性能的提升

4.2 局部存储器深度对矩阵乘协处理器性能的影响

Fig. 7 Effect of local memory depth on coprocessor performance图7 局部存储器深度对矩阵乘协处理器性能的影响

局部存储器的深度直接影响协处理器硬件实现的开销，因而是矩阵乘协处理器结构设计中一个重要的参数.实验在128，1 024，8 192，16 384这4种矩阵规模下，分别将计算单元局部存储器的深度h从128增加至8 192，统计矩阵乘协处理器实际计算性能的变化情况，实验结果如图7所示:

从图7可以看出，4种矩阵规模下，矩阵乘协处理器的实际计算性能均随局部存储器深度的增加而提升，局部存储器深度增加大于2 048时，矩阵乘协处理器的性能基本保持不变.协处理器实际计算性能随局部存储器深度的变化可以根据式(2)中计算完全隐藏存储访问延迟的约束条件进行分析.当St

4.3 矩阵规模对矩阵乘协处理器计算效率的影响

尽管本文矩阵乘协处理器的设计面向数值计算领域的大规模应用，矩阵规模对协处理器实际计算性能和效率的影响仍然值得关注.实验在应用双缓冲优化、局部存储器深度设置为2 048时，统计了不同矩阵规模下，矩阵乘协处理器实际的计算性能和效率，以及存储访问延迟占矩阵乘运算核心总执行时间的比例，实验结果如表2所示:

Table 2 Effect of Matrix Size on Coprocessor Performance

从表2中的数据可以看出，随着矩阵规模的增加，矩阵乘协处理器的实际计算效率不断提高：矩阵规模为128时，协处理器的计算效率仅为6.25%；矩阵规模增加至4 096时，协处理器的计算效率达到96.97%；矩阵规模大于等于8 192时，协处理器的计算效率可以达到99%.

对协处理器计算效率与矩阵规模相关性的分析如下：当PE数Sp=1 024、局部存储器深度St=2 048时，式(3)表示的计算完全隐藏访存延迟的约束条件是否满足仅取决于计算中的矩阵规模K.随着矩阵规模的增大，计算可以更好地隐藏存储访问的延迟，协处理器的计算效率也进一步提升.表2中存储访问延迟占矩阵乘总执行时间的比例验证了上述分析，即较大的矩阵规模对应较低的存储访问延迟占比.

4.4 功能正确性验证与硬件实现开销的评估

本文基于Verilog硬件描述语言实现了双精度浮点矩阵乘协处理器.其中，线性阵列计算结构中各计算单元采用双缓冲优化技术，PE数设置为1 024、局部存储器的深度设置为2 048.

对于矩阵乘协处理器功能正确性的验证，本文以Linpack性能测试程序HPL[25]随机生成的双精度浮点矩阵数据作为测试数据集，通过调用OpenBLAS[26-27]线性代数计算函数库实现了双精度浮点矩阵乘运算核心.本文将OpenBLAS的计算结果与矩阵乘协处理器的仿真计算结果进行了比对，从而验证了协处理器功能实现的正确性.

对于矩阵乘协处理器硬件开销的评估，本文基于Synopsys Design Compiler设计工具在1.5 GHz、主流工艺下对协处理器的硬件设计进行了逻辑综合，综合过程中的参数配置如表3所示：

Table 3 Configurations for Synopsys Design Compiler表3 逻辑综合的参数设置

根据Design Compiler输出的综合结果，矩阵乘协处理器的硬件设计完全满足1.5 GHz频率下的时序要求.其功耗为38.57 W，总面积为116.26 mm2.其中，组合逻辑的面积为16.78 mm2，非组合逻辑的面积为99.48 mm2.

本文将矩阵乘协处理器与采用同等工艺的NVIDIA K40 GPGPU[19]芯片进行了对比，比较的内容包括时钟频率、峰值计算性能，进行矩阵乘运算的实际性能和效率以及硬件开销等，结果如表4所示.针对双精度浮点矩阵乘运算，本文协处理器的实际性能是K40 GPGPU芯片的1.95倍，而面积开销仅为后者的21.05%(不考虑GPGPU芯片中PCIe、图形显示和存储控制器等接口的硬件开销).

Table 4 Coprocessor versus K40 GPGPU on Performance

5 结束语

矩阵乘运算广泛应用于科学与工程领域.因计算复杂度较高，且在CPU，GPGPU，FPGA等现有计算平台上的性能和效率受限，双精度浮点矩阵乘运算往往成为大规模数值计算应用的性能瓶颈.近年来，面向应用的硬件定制成为高性能、高效计算结构设计的重要趋势.本文将这一结构设计方法应用到矩阵乘的加速计算中，设计了双精度浮点矩阵乘定制协处理器和基于协处理器的加速计算系统.

本文双精度浮点矩阵乘协处理器的设计以线性阵列计算结构为基础，应用了双缓冲技术并对存储访问调度进行了优化.同时，本文构建了针对协处理器的性能模型，并基于性能模型深入分析了协处理器的各结构设计参数对其计算性能和效率的影响.此外，本文还验证了矩阵乘协处理器功能实现的正确性，并在主流工艺下评估了其硬件开销.本文双精度浮点矩阵乘协处理器的实际计算性能可达3 TFLOPS，计算效率最高可达99%，有效实现了大规模数值计算应用的加速.与同等工艺的GPGPU计算平台相比，本文设计的矩阵乘协处理器在计算性能和硬件开销方面均具有明显优势.

本文探索了硬件定制的结构设计在高性能计算中的应用.后续研究将基于BLAS库和Linpack性能测试程序评估本文双精度浮点矩阵乘协处理器的应用对大规模计算系统性能和效率的提升效果；同时分析协处理器实际性能随阵列计算结构中计算单元数增加的可扩展性，从而进一步优化协处理器的结构设计.另外，协处理器的应用适应性也是未来研究工作的重点.