CN113987414B - 基于ARMv8多核处理器的小型和不规则矩阵乘优化方法 - Google Patents

基于ARMv8多核处理器的小型和不规则矩阵乘优化方法 Download PDF

Info

Publication number
CN113987414B
CN113987414B CN202111296475.6A CN202111296475A CN113987414B CN 113987414 B CN113987414 B CN 113987414B CN 202111296475 A CN202111296475 A CN 202111296475A CN 113987414 B CN113987414 B CN 113987414B
Authority
CN
China
Prior art keywords
matrix
armv8
multiplication
core processor
mode
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202111296475.6A
Other languages
English (en)
Other versions
CN113987414A (zh
Inventor
董德尊
方建滨
杨维玲
苏醒
庞征斌
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202111296475.6A priority Critical patent/CN113987414B/zh
Publication of CN113987414A publication Critical patent/CN113987414A/zh
Application granted granted Critical
Publication of CN113987414B publication Critical patent/CN113987414B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F7/00—Methods or arrangements for processing data by operating upon the order or content of the data handled
    • G06F7/38—Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation
    • G06F7/48—Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation using non-contact-making devices, e.g. tube, solid state device; using unspecified devices
    • G06F7/57—Arithmetic logic units [ALU], i.e. arrangements or devices for performing two or more of the operations covered by groups G06F7/483 – G06F7/556 or for performing logical operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F9/00—Arrangements for program control, e.g. control units
    • G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
    • G06F9/30—Arrangements for executing machine instructions, e.g. instruction decode
    • G06F9/30098—Register arrangements
    • G06F9/30101—Special purpose registers
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D10/00—Energy efficient computing, e.g. low power processors, power management or thermal management

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Computing Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开了一种基于ARMv8多核处理器的小型和不规则矩阵乘优化方法，利用ARMv8多核处理器实现,其步骤包括：建立矩阵存储空间，用于存放矩阵A与矩阵B相乘得到的结果矩阵C；ARMv8多核处理器对矩阵B执行打包操作，并将打包操作和小型矩阵乘的计算操作同步进行。本方法对于不同模式的矩阵乘选择不同打包策略来节省打包开销，使用更加高效的边缘微内核来处理边缘案例，此外还采用更加合理的并行化方法来并行化矩阵乘，这些对ARMv8多核处理器中的小型和不规则矩阵乘的性能进行了大幅优化，这能促进ARMV8多核处理器上其他实际应用的发展。

Description

基于ARMv8多核处理器的小型和不规则矩阵乘优化方法
技术领域
本发明涉及高性能计算领域，尤其涉及一种基于ARMv8多核处理器的小型和不规则矩阵乘优化方法。
背景技术
通用矩阵乘法(GEMM)是从传统科学模拟到新兴的深度学习高性能计算(HPC)应用的基本构建块。如何对GEMM进行优化，是一个被大量研究的领域，但现有的线性代数库方法主要针对大型和规则形状的GEMM(即当矩阵的两个维度或多或少相同时)。
由于HPC工作负载的多样性和不断发展，GEMM内核的输入矩阵的大小和形状可能会因所使用的应用程序算法和输入数据而异。例如，计算流体动力学(CFD)中，如有限元方法和波动方程，通常采用在小矩阵上运行的GEMM实现，以在现代多核系统上实现良好的可扩展性能。例如，流行的分子动力学模拟器CP2K在大小为5×5和23×23的矩阵上执行GEMM。又例如，用于CFD的Nek5000高阶求解器的内核严重依赖于8×8矩阵的GEMMs计算。除了这些传统的HPC应用程序之外，实现深度学习和机器学习方法等HPC应用程序建立在小的GEMM内核上，其中一些数据分析算法还需要对不规则形状的矩阵进行操作，而其中不同的矩阵维度的大小有显著差异。例如,ResNet深度神经网络的卷积核使用的GEMM计算一维等于64而另一维大于3000的矩阵。
这些利用HPC求解的矩阵特性使得我们需要对GEMM计算方式进行优化。尽管像OpenBLAS和BLIS这样的传统线性代数库可以在大型和规则形状的GEMM上提供接近最佳的性能，但它们在小型GEMM上的性能通常很差。虽然现有方法在x86和GPU架构上提供了可喜的结果，但现有的解决方案不足以优化基于ARMv8的CPU架构上的小型和不规则形状GEMM。
发明内容
GEMM是矩阵的乘加运算操作，定义为C＝αA·B+βC，其中A和B是矩阵输入，α和β是标量输入，C是预先存在的矩阵被输出覆盖的矩阵。矩阵A表示为M×K的输入矩阵，即M行K列，矩阵B的大小为K×N,即即K行N列,C的大小是M×N,为输出矩阵。针对现有的BLAS(BasicLinear Algebra Subprograms)库难以适应小型和不规则GEMM，且在基于ARMv8的CPU架构上求解小型和不规则GEMM效率偏低的问题，本发明公开了一种基于ARMv8多核处理器的小型和不规则矩阵乘优化方法。
本发明公开了一种基于ARMv8多核处理器的小型矩阵乘优化方法，利用ARMv8多核处理器实现,其步骤包括：
S1，对于参与相乘运算的两个矩阵，乘号左侧的矩阵称为矩阵A，乘号右侧的矩阵称为矩阵B，判断矩阵A和矩阵B是否被转置，在相乘运算中，若矩阵没被转置称为N模式，矩阵被转置称为T模式，若矩阵A为N模式，矩阵B为T模式，则将该小型矩阵乘法称为NT模式，若矩阵A为N模式，矩阵B为N模式，则将该小型矩阵乘法称为NN模式；建立矩阵存储空间，用于存放矩阵A与矩阵B相乘得到的结果矩阵C；
S2,对于NT模式的小型矩阵乘，ARMv8多核处理器对矩阵B执行打包操作，并将打包操作和小型矩阵乘的计算操作同步进行。将ARMv8多核处理器中负责打包的微内核称为打包微内核，将ARMv8多核处理器中负责计算的微内核称为主微内核。为了将打包的访存开销隐藏在计算过程中，实现矩阵乘的微内核要在ARMv8处理器器中需要具有足够高的计算访存比CMR，其中CMR的计算公式和约束条件为：
Figure BDA0003336782700000021
其中mr和nr分别为ARMv8多核处理器上的向量寄存器所能存储的最大矩阵的行数和列数，同时也是主微内核在向量寄存器上所占用的存储空间矩阵的行数和列数。C0为计算访存比CMR的下限。ARMv8多核处理器提供了32个128bit的向量寄存器，该32个128bit的向量寄存器命名为寄存器V0、V1、V2、…、V31，因此mr和nr还需要满足的约束条件为：
Figure BDA0003336782700000031
其中j表示ARMv8多核处理器的一个向量寄存器可以加载数据元素的个数，％表示求模。以计算访存比的约束条件、mr和nr的约束条件为目标函数，使用拉格朗日乘子法，求解得到mr0和nr0，作为主微内核在向量寄存器上所占用的存储空间矩阵的行数和列数的取值。
在NT模式下实现小型矩阵乘时，用于对矩阵B进行打包的微内核，所占用向量寄存器的存储空间矩阵的维度为mp×np，其中mp＝mr0,np＝nr0/N0,经过N0次调用打包微内核，得到一个主微内核实现小型矩阵乘计算所需的数据。在对矩阵B进行打包的过程中，ARMv8多核处理器沿着矩阵A的行方向访问矩阵A，沿着矩阵B的列方向访问矩阵B，使用内积公式更新矩阵C。ARMv8多核处理器在沿上述方向对矩阵A和矩阵B分别进行访问的过程中，使用mp次加载指令将矩阵A的元素存储到从V0开始依次至V(mp-1)的mp个向量寄存器中，并使用np次加载指令将矩阵B的元素存储到从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中。
ARMv8多核处理器在微内核打包过程中执行mp×np次向量乘加融合指令FMA，以产生矩阵C的N0×mp×np个中间结果元素，这些元素存储在从V(mp+np)依次至V31的向量寄存器中，同时将从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中的前N0个元素存储到内存空间Bc，Bc为临时开辟的用于存放矩阵B的部分元素的连续内存空间，该N0个元素之间在内存空间Bc上的物理距离为nr0(N0×np)个元素，从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中的相同位置的元素存储到内存空间Bc的相邻位置上。所述的用于打包的微内核中实现小型矩阵乘计算使用内积公式，主微内核实现小型矩阵乘计算使用外积公式。
对于NN模式的小型矩阵乘，如果矩阵B的大小超过ARMv8多核处理器的一级缓存后，ARMv8多核处理器对矩阵B执行打包操作，将其打包到内存空间Bc中，并将打包操作和小型矩阵乘的计算操作同步进行，该打包操作和计算操作均通过迭代来实现，否则直接执行步骤S3。每次迭代中，将矩阵A的元素存储到从V0开始依次至V(mp-1)的mp个向量寄存器中，将矩阵B的元素存储到从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中。从V0开始依次至V(mp-1)的mp个向量寄存器，在每次迭代中，每个向量寄存器都只有一个元素用来参与小型矩阵乘的计算操作，从而需要N0次迭代，才能使得每个向量寄存器中的元素均参与到小型矩阵乘的计算操作中。而从V(mp)开始依次至V(mp-1+np)的np个向量寄存器的所有元素在一次迭代中就使用完毕，在下一轮的迭代中又需要重新从矩阵B中取数据。当所有迭代完成之后，矩阵C的1行至mr0行已经被更新，内存空间Bc已经被矩阵B的打包数据填充完整。
S3,ARMv8多核处理器使用内存空间Bc中的已经打包的数据来更新矩阵C的mr0行到M行的元素。步骤S2得到的内存空间Bc来提供计算矩阵C的mr0+1行至M行所需的数据，此步骤对于NN模式与NT模式，其实现流程是相同的。矩阵C的mr0+1行至M行的计算操作通过迭代来完成。ARMv8多核处理器从矩阵A中选择mr0行，对其每一行分别取N0个元素，存放到从V0开始依次至V(mr0-1)的mr0个向量寄存器中，该操作称为取数操作，该mr0行的地址是不连续的，相互之间的距离为N0×K个字节，在N0次迭代结束后再次执行该取数操作。在每次迭代中，内存空间Bc中的元素被取出并存储到从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中，将矩阵A中取出的每个元素作为标量，每次迭代中从内存空间Bc中取出nr0个元素作为向量，二者依次相乘，得到N0×mp×np个元素，存储在从V(mp+np)依次至V31的向量寄存器中。在矩阵C的mr0+1行至M行的计算过程中，有部分执行数据预取操作的指令被插入计算与访存指令的执行时间间隙中。数据预取操作是指在计算指令执行前，将数据从内存预先取出到缓存中。
S4,当M％mr0不为0或N％nr0不为0时，调用ARMv8多核处理器中处理边缘数据的微内核，简称边缘微内核。对于边缘微内核的实现，其计算指令和访存指令需要按照计算需求进行排布，以确保计算指令所需的数据在执行计算指令前，已存放至向量寄存器中。计算指令FMA所需的数据提前通过取数据指令加载到向量寄存器中。
S5,在沿着K维度进行的循环迭代结束后，将ARMv8多核处理器中存储矩阵C的相关向量寄存器中的数据存回到内存中的矩阵C。
本发明公开了一种基于ARMv8多核处理器的不规则矩阵乘优化方法，利用ARMv8多核处理器实现,其步骤包括：
S1，对于参与相乘运算的两个矩阵，乘号左侧的矩阵称为矩阵A，乘号右侧的矩阵称为矩阵B，判断矩阵A和矩阵B是否被转置，在相乘运算中，若矩阵没被转置称为N模式，矩阵被转置称为T模式，若矩阵A为N模式，矩阵B为T模式，则将该不规则矩阵乘法称为NT模式，若矩阵A为N模式，矩阵B为N模式，则将该不规则矩阵乘法称为NN模式；建立矩阵存储空间，用于存放矩阵A与矩阵B相乘得到的结果矩阵C；
S2,确定对不规则形状矩阵乘的行方向和列方向的数据进行并行化处理所需要的并行线程数，为了保证并行线程之间的工作平衡，对矩阵C的行方向和列方向的数据分别进行并行化。将矩阵C划分为二维的子块网格，每个并行线程更新一个子块网格。当使用T个并行线程对矩阵C的计算负载进行划分时，每个并行线程将执行
Figure BDA0003336782700000051
次计算矩阵乘法的操作。每个并行线程所需的内存访问次数为
Figure BDA0003336782700000052
其中Tm和Tn分别是分配给M维和N维度数据的并行线程数，Tm×Tn＝T。对于一个子块网格的CMR，其计算公式为：
Figure BDA0003336782700000053
为了最大化CMR，使用算术-几何平均不等式，得到该CMR的取值范围为：
Figure BDA0003336782700000054
考虑S2中执行微内核打包的开销，取Tn的上界值，即
Figure BDA0003336782700000055
以最大化CMR，
Figure BDA0003336782700000056
表示向上取整。
S3，在NN模式的不规则形状矩阵乘中，如果矩阵B的大小超过最后一级缓存的容量，将t×nr0个元素打包到内存空间Bc中，其中t为自然数。打包微内核中的mr和nr的取值与主内核中对应值的取值相同。小型和不规则形状的矩阵乘的t值分别设置为0和1。
S4,当M％mr0不为0或N％nr0不为0时，调用ARMv8多核处理器中处理边缘数据的微内核，简称边缘微内核。对于边缘微内核的实现，其计算指令和访存指令需要按照计算需求进行排布，以确保计算指令所需的数据在执行计算指令前，已存放至向量寄存器中。计算指令FMA所需的数据提前通过取数据指令加载到向量寄存器中。
S5,在沿着K维度进行的循环迭代结束后，将ARMv8多核处理器中存储矩阵C的相关向量寄存器中的数据存回到内存中的矩阵C。
本发明的有益效果为：
本方法对于不同模式的矩阵乘选择不同打包策略来节省打包开销，使用更加高效的边缘微内核来处理边缘案例，此外还采用更加合理的并行化方法来并行化矩阵乘，这些对ARMv8多核处理器中的小型和不规则矩阵乘的性能进行了大幅优化，这能促进ARMV8多核处理器上其他实际应用的发展。
附图说明
图1为NN模式下小型矩阵乘的内核设计；
图2为NT模式下矩阵乘的打包微内核设计；
图3为边缘微内核的设计；
图4为NT模式下不规则矩阵乘的微内核流程图；
图5为单线程小型矩阵乘的性能(热cache)；
图6为单线程小型矩阵乘的性能(冷cache)；
图7为多线程不规则矩阵乘在Phytium 2000+上的性能；
图8为多线程不规则矩阵乘在KP920和Thunder X2上的性能；
图9为LibShalom对于CP2K中使用矩阵的性能；
图10为LibShalom对于VGG中使用矩阵的性能。
具体实施方式
为了更好的了解本发明内容，这里给出一个实施例。
图1为NN模式下小型矩阵乘的内核设计；图2为NT模式下矩阵乘的打包微内核设计；图3为边缘微内核的设计；图4为NT模式下不规则矩阵乘的微内核流程图；图5为单线程小型矩阵乘的性能(热cache)；图6为单线程小型矩阵乘的性能(冷cache)；图7为多线程不规则矩阵乘在Phytium2000+上的性能；图8为多线程不规则矩阵乘在KP920和Thunder X2上的性能；图9为LibShalom对于CP2K中使用矩阵的性能；图10为LibShalom对于VGG中使用矩阵的性能。
针对现有的BLAS(Basic Linear Algebra Subprograms)库难以适应小型和不规则GEMM，且在基于ARMv8的CPU架构上求解小型和不规则GEMM效率偏低的问题，本发明公开了一种ARMv8多核处理器上小型和不规则形状矩阵乘的优化方法。
本发明公开了一种基于ARMv8多核处理器的小型矩阵乘优化方法，利用ARMv8多核处理器实现,其步骤包括：
S1，对于参与相乘运算的两个矩阵，乘号左侧的矩阵称为矩阵A，乘号右侧的矩阵称为矩阵B，判断矩阵A和矩阵B是否被转置，在相乘运算中，若矩阵没被转置称为N模式，矩阵被转置称为T模式，若矩阵A为N模式，矩阵B为T模式，则将该小型矩阵乘法称为NT模式，若矩阵A为N模式，矩阵B为N模式，则将该小型矩阵乘法称为NN模式；建立矩阵存储空间，用于存放矩阵A与矩阵B相乘得到的结果矩阵C；
S2,对于NT模式的小型矩阵乘，ARMv8多核处理器对矩阵B执行打包操作，并将打包操作和小型矩阵乘的计算操作同步进行。将ARMv8多核处理器中负责打包的微内核称为打包微内核，将ARMv8多核处理器中负责计算的微内核称为主微内核。为了将打包的访存开销隐藏在计算过程中，实现矩阵乘的微内核要在ARMv8处理器器中需要具有足够高的计算访存比CMR，其中CMR的计算公式和约束条件为：
Figure BDA0003336782700000071
其中mr和nr分别为ARMv8多核处理器上的向量寄存器所能存储的最大矩阵的行数和列数，同时也是主微内核在向量寄存器上所占用的存储空间矩阵的行数和列数。C0为计算访存比CMR的下限，其与ARMv8多核处理器的向量寄存器容量有关。ARMv8多核处理器提供了32个128bit的向量寄存器，该32个128bit的向量寄存器命名为寄存器V0、V1、V2、…、V31，因此mr和nr还需要满足的约束条件为：
Figure BDA0003336782700000081
其中j表示ARMv8多核处理器的一个向量寄存器可以加载数据元素的个数，％表示求模。对于单精度浮点数而言，j＝4。以计算访存比的约束条件、mr和nr的约束条件为目标函数，使用拉格朗日乘子法，求解得到mr0和nr0，作为主微内核在向量寄存器上所占用的存储空间矩阵的行数和列数的取值。
为了方便在后续计算中使用该主微内核，在NT模式下实现小型矩阵乘时，用于对矩阵B进行打包的微内核，所占用向量寄存器的存储空间矩阵的维度为mp×np，其中mp＝mr0,np＝nr0/N0,经过N0次调用打包微内核，得到一个主微内核实现小型矩阵乘计算所需的数据，其流程如图1所示。在对矩阵B进行打包的过程中，ARMv8多核处理器沿着矩阵A的行方向访问矩阵A，沿着矩阵B的列方向访问矩阵B，使用内积公式更新矩阵C。ARMv8多核处理器在沿上述方向对矩阵A和矩阵B分别进行访问的过程中，使用mp次加载指令将矩阵A的元素存储到从V0开始依次至V(mp-1)的mp个向量寄存器中，并使用np次加载指令将矩阵B的元素存储到从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中。
ARMv8多核处理器在微内核打包过程中执行mp×np次向量乘加融合指令FMA，以产生矩阵C的N0×mp×np个中间结果元素，这些元素存储在从V(mp+np)依次至V31的向量寄存器中，同时将从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中的前N0个元素存储到内存空间Bc，Bc为临时开辟的用于存放矩阵B的部分元素的连续内存空间，该N0个元素之间在内存空间Bc上的物理距离为nr0(N0×np)个元素，从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中的相同位置的元素存储到内存空间Bc的相邻位置上。所述的从V(mp-1)开始依次至V(mp-2+np)的np个向量寄存器中的相同位置的元素，例如向量V7的4个元素分别编号为a0,a1,a2,a3；向量V8的4个元素分别编号为b0,b1,b2,b3；这里的a0与b0就属于相同位置的元素。所述的用于打包的微内核中实现小型矩阵乘计算使用内积公式，主微内核实现小型矩阵乘计算使用外积公式。
对于NN模式的小型矩阵乘，如果矩阵B的大小超过ARMv8多核处理器的一级缓存后，ARMv8多核处理器对矩阵B执行打包操作，将其打包到内存空间Bc中，并将打包操作和小型矩阵乘的计算操作同步进行,其流程如图2所示，该打包操作和计算操作均通过迭代来实现，否则直接执行步骤S3。每次迭代中，将矩阵A的元素存储到从V0开始依次至V(mp-1)的mp个向量寄存器中，将矩阵B的元素存储到从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中。从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中，每个向量寄存器可以装下N0个单精度浮点元素。从V0开始依次至V(mp-1)的mp个向量寄存器，在每次迭代中，每个向量寄存器都只有一个元素用来参与小型矩阵乘的计算操作，从而需要N0次迭代，才能使得每个向量寄存器中的元素均参与到小型矩阵乘的计算操作中。而从V(mp)开始依次至V(mp-1+np)的np个向量寄存器的所有元素在一次迭代中就使用完毕，在下一轮的迭代中又需要重新从矩阵B中取数据。上述的迭代指的是沿着矩阵A和矩阵B的K个维度进行循环的迭代。当所有迭代完成之后，矩阵C的1行至mr0行已经被更新，内存空间Bc已经被矩阵B的打包数据填充完整。
S3,ARMv8多核处理器使用内存空间Bc中的已经打包的数据来更新矩阵C的mr0行到M行的元素。步骤S2得到的内存空间Bc来提供计算矩阵C的mr0+1行至M行所需的数据，此步骤对于NN模式与NT模式，其实现流程是相同的。矩阵C的mr0+1行至M行的计算操作通过迭代来完成。ARMv8多核处理器从矩阵A中选择mr0行，对其每一行分别取N0个元素，存放到从V0开始依次至V(mr0-1)的mr0个向量寄存器中，该操作称为取数操作，该mr0行的地址是不连续的，相互之间的距离为N0×K个字节，在N0次迭代结束后再次执行该取数操作。在每次迭代中，内存空间Bc中的元素被取出并存储到从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中，将矩阵A中取出的每个元素作为标量，每次迭代中从内存空间Bc中取出nr0个元素作为向量，二者依次相乘，得到N0×mp×np个元素，存储在从V(mp+np)依次至V31的向量寄存器中。在矩阵C的mr0+1行至M行的计算过程中，有部分执行数据预取操作的指令被插入计算与访存指令的执行时间间隙中。数据预取操作是指在计算指令执行前，将数据从内存预先取出到缓存中。在迭代操作中，使用的循环展开因子为8，即要把8轮迭代用代码实现出来，然后嵌套在循环嵌套中，而不是只实现1次循环的迭代，然后就用循环嵌套。
S4,在经过前面的更新步骤后，矩阵C可能还会剩下一些边缘数据没有被更新，这是由于M％mr0或N％nr0不为0，即结果矩阵C的大小M×N不能被主微内核的大小mr×nr整除。当M％mr0不为0或N％nr0不为0时，调用ARMv8多核处理器中处理边缘数据的微内核，简称边缘微内核。对于边缘微内核的实现，其计算指令和访存指令需要按照计算需求进行排布，以确保计算指令所需的数据在执行计算指令前，已存放至向量寄存器中。计算指令FMA所需的数据提前通过取数据指令加载到向量寄存器中，如图3所示。这里的FMA是处理器的乘加融合指令，load是取数据指令。
S5,在沿着K维度进行的循环迭代结束后，将ARMv8多核处理器中存储矩阵C的相关向量寄存器中的数据存回到内存中的矩阵C。
对于基于ARMv8多核处理器的不规则形状矩阵乘的实现，
本发明公开了一种基于ARMv8多核处理器的不规则矩阵乘优化方法，利用ARMv8多核处理器实现,其步骤包括：
S1，对于参与相乘运算的两个矩阵，乘号左侧的矩阵称为矩阵A，乘号右侧的矩阵称为矩阵B，判断矩阵A和矩阵B是否被转置，在相乘运算中，若矩阵没被转置称为N模式，矩阵被转置称为T模式，若矩阵A为N模式，矩阵B为T模式，则将该小型矩阵乘法称为NT模式，若矩阵A为N模式，矩阵B为N模式，则将该小型矩阵乘法称为NN模式；建立矩阵存储空间，用于存放矩阵A与矩阵B相乘得到的结果矩阵C；
S2,确定对不规则形状矩阵乘的行方向和列方向的数据进行并行化处理所需要的并行线程数，为了保证并行线程之间的工作平衡，对矩阵C的行方向和列方向的数据分别进行并行化。将矩阵C划分为二维的子块网格，每个并行线程更新一个子块网格。当使用T个并行线程对矩阵C的计算负载进行划分时，每个并行线程将执行
Figure BDA0003336782700000111
次计算矩阵乘法的操作。每个并行线程所需的内存访问次数为
Figure BDA0003336782700000112
其中Tm和Tn分别是分配给M维和N维度数据的并行线程数，Tm×Tn＝T。对于一个子块网格的CMR，其计算公式为：
Figure BDA0003336782700000113
为了最大化CMR，使用算术-几何平均不等式，得到该CMR的取值范围为：
Figure BDA0003336782700000114
如果
Figure BDA0003336782700000115
则方程的两边将相等。换句话说，当
Figure BDA0003336782700000116
时，CMR取得最大值。考虑S2中执行微内核打包的开销，取Tn的上界值，即
Figure BDA0003336782700000117
Figure BDA0003336782700000118
以最大化CMR。我们注意到TmodTn＝0，以确保内核数量可以在并行线程之间平均分配。
Figure BDA0003336782700000119
表示向上取整。
S3，在NN模式的不规则形状矩阵乘中，如果矩阵B的大小超过最后一级缓存的容量，将t×nr0个元素打包到内存空间Bc中，其中t为自然数。打包微内核中的mr和nr的取值与主内核中对应值的取值相同(即，对于FP32，mr＝7和nr＝12)。小型和不规则形状的矩阵乘的t值分别设置为0和1。这意味着前者在每次迭代中只执行图4中的step①，不规则形状矩阵乘需要执行step①和②。
S4,在经过前面的更新步骤后，矩阵C可能还会剩下一些边缘数据没有被更新，这是由于M％mr0或N％nr0不为0，即结果矩阵C的大小M×N不能被主微内核的大小mr×nr整除。当M％mr0不为0或N％nr0不为0时，调用ARMv8多核处理器中处理边缘数据的微内核，简称边缘微内核。对于边缘微内核的实现，其计算指令和访存指令需要按照计算需求进行排布，以确保计算指令所需的数据在执行计算指令前，已存放至向量寄存器中。计算指令FMA所需的数据提前通过取数据指令加载到向量寄存器中，如图3所示。这里的FMA是处理器的乘加融合指令，load是取数据指令。
S5,在沿着K维度进行的循环迭代结束后，将ARMv8多核处理器中存储矩阵C的相关向量寄存器中的数据存回到内存中的矩阵C。
设计LibShalom的目的是为了克服在ARMv8处理器上优化小型和不规则GEMM性能不高的问题，该方法将计算与打包操作重叠看，对边缘微内核进行了高度优化，使用更加合理地并行化方法来开展优化。这里对LibShalom是否到达上述目的进行了有效性评估。图5和图6在Phyituym 2000+、KP920和ThunderX2这三个ARMv8处理器上评测单线程小型GEMM的性能。图5和图6中的(a)、(b)和(c)子图分别是在Phyituym 2000+、KP920和ThunderX2这三个ARMv8处理器上评测单线程小型GEMM的性能。图7和图8展示了多线程不规则GEMM的性能。图9展示了对于CP2K中使用到的矩阵的性能，其中(a)、(b)和(c)子图分别是在Phyituym2000+、KP920和ThunderX2这三个ARMv8处理器上得到的性能。图10显示了对于卷积神经网络VGG中使用到的矩阵的性能，其中(a)、(b)和(c)子图分别是在Phyituym 2000+、KP920和ThunderX2这三个ARMv8处理器上得到的性能。
图5和6给出了LibShalom(本发明方法)、OpenBLAS(传统的高性能BLAS库)、BLIS(在OpenBLAS之后发布的高性能BLAS库)、ARMPL(ARM公司官方的BLAS库)、LIBXSMM(Intel专门为小型GEMM设计的库)、BLASFEO(为多种处理器设计的针对小型GEMMBLAS库)的性能柱状图。实验中的数据规模M×N×K从8×8×8到120×120×120，运行20次计算平均值。图5和图6分别是热cache和冷cache的实验结果。基于图5到图6，可以得出下述结论：(1)矩阵越小的情况下，LibShalom的性能优势越大，这是因为打包开销在矩阵越小时开销占比越大。(2)在矩阵规模逐渐增大的情况下，LibShalom相比于第二好的BLASFEO也能取得加速，这得益于高度优化的微内核。
图7和图8展示了上述方法对于多线程不规则GEMM的性能，由于LIBXSMM和BLASFEO没有对不规则GEMM做优化，因此没有被操作参照对象。基于图5到图6，可以得出下述结论：(1)LibShalom采用的并行化方法能更好的使用处理器的各个核心，在Phyituym 2000+、KP920和ThunderX2，对比性能第二好的BLIS分别能取得1.8、1.6和1.3倍的加速；(2)与NN模式相比，LibShalom为在NT模式下运行的不规则形状的GEMM提供了更高的性能。这是因为，与NT模式不同，在N模式下,无法沿着K维度连续访问矩阵B的元素。总的来说，LibShalom是普遍适用的，可以在代表性的ARMv8处理器上提供可移植的性能。
图9展示了LibShalom对于CP2K程序中使用矩阵的性能，可以看出对于在实际应用中出现的矩阵，LibShalom依然能取得加速，这说明其可以被移植到ARM上的其他应用中。
图10展示了LibShalom对于卷积神经网络VGG中使用矩阵的性能，对于在新兴的深度学习计算负载，LibShalom依然能取得加速，这说明其可以用于加速深度学习应用。
以上所述仅为本申请的实施例而已，并不用于限制本申请。对于本领域技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本申请的权利要求范围之内。

Claims (2)

1.一种基于ARMv8多核处理器的小型矩阵乘优化方法，其特征在于，利用ARMv8多核处理器实现,其步骤包括：
S1，对于参与相乘运算的两个矩阵，乘号左侧的矩阵称为矩阵A，乘号右侧的矩阵称为矩阵B，判断矩阵A和矩阵B是否被转置，在相乘运算中，若矩阵没被转置称为N模式，矩阵被转置称为T模式，若矩阵A为N模式，矩阵B为T模式，则将该小型矩阵乘法称为NT模式，若矩阵A为N模式，矩阵B为N模式，则将该小型矩阵乘法称为NN模式；建立矩阵存储空间，用于存放矩阵A与矩阵B相乘得到的结果矩阵C；
S2,对于NT模式的小型矩阵乘，ARMv8多核处理器对矩阵B执行打包操作，并将打包操作和小型矩阵乘的计算操作同步进行；将ARMv8多核处理器中负责打包的微内核称为打包微内核，将ARMv8多核处理器中负责计算的微内核称为主微内核；为了将打包的访存开销隐藏在计算过程中，实现矩阵乘的微内核要在ARMv8处理器器中需要具有足够高的计算访存比CMR，其中CMR的计算公式和约束条件为：
Figure FDA0003336782690000011
其中mr和nr分别为ARMv8多核处理器上的向量寄存器所能存储的最大矩阵的行数和列数，同时也是主微内核在向量寄存器上所占用的存储空间矩阵的行数和列数；C0为计算访存比CMR的下限；ARMv8多核处理器提供了32个128bit的向量寄存器，该32个128bit的向量寄存器命名为寄存器V0、V1、V2、…、V31，因此mr和nr还需要满足的约束条件为：
Figure FDA0003336782690000012
其中j表示ARMv8多核处理器的一个向量寄存器可以加载数据元素的个数，％表示求模；以计算访存比的约束条件、mr和nr的约束条件为目标函数，使用拉格朗日乘子法，求解得到mr0和nr0，作为主微内核在向量寄存器上所占用的存储空间矩阵的行数和列数的取值；
在NT模式下实现小型矩阵乘时，用于对矩阵B进行打包的微内核，所占用向量寄存器的存储空间矩阵的维度为mp×np，其中mp＝mr0,np＝nr0/N0,经过N0次调用打包微内核，得到一个主微内核实现小型矩阵乘计算所需的数据；在对矩阵B进行打包的过程中，ARMv8多核处理器沿着矩阵A的行方向访问矩阵A，沿着矩阵B的列方向访问矩阵B，使用内积公式更新矩阵C；ARMv8多核处理器在沿上述方向对矩阵A和矩阵B分别进行访问的过程中，使用mp次加载指令将矩阵A的元素存储到从V0开始依次至V(mp-1)的mp个向量寄存器中，并使用np次加载指令将矩阵B的元素存储到从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中；
ARMv8多核处理器在微内核打包过程中执行mp×np次向量乘加融合指令FMA，以产生矩阵C的N0×mp×np个中间结果元素，这些元素存储在从V(mp+np)依次至V31的向量寄存器中，同时将从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中的前N0个元素存储到内存空间Bc，Bc为临时开辟的用于存放矩阵B的部分元素的连续内存空间，该N0个元素之间在内存空间Bc上的物理距离为nr0(N0×np)个元素，从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中的相同位置的元素存储到内存空间Bc的相邻位置上；所述的用于打包的微内核中实现小型矩阵乘计算使用内积公式，主微内核实现小型矩阵乘计算使用外积公式；
对于NN模式的小型矩阵乘，如果矩阵B的大小超过ARMv8多核处理器的一级缓存后，ARMv8多核处理器对矩阵B执行打包操作，将其打包到内存空间Bc中，并将打包操作和小型矩阵乘的计算操作同步进行，该打包操作和计算操作均通过迭代来实现，否则直接执行步骤S3；每次迭代中，将矩阵A的元素存储到从V0开始依次至V(mp-1)的mp个向量寄存器中，将矩阵B的元素存储到从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中；从V0开始依次至V(mp-1)的mp个向量寄存器，在每次迭代中，每个向量寄存器都只有一个元素用来参与小型矩阵乘的计算操作，从而需要N0次迭代，才能使得每个向量寄存器中的元素均参与到小型矩阵乘的计算操作中；而从V(mp)开始依次至V(mp-1+np)的np个向量寄存器的所有元素在一次迭代中就使用完毕，在下一轮的迭代中又需要重新从矩阵B中取数据；当所有迭代完成之后，矩阵C的1行至mr0行已经被更新，内存空间Bc已经被矩阵B的打包数据填充完整；
S3,ARMv8多核处理器使用内存空间Bc中的已经打包的数据来更新矩阵C的mr0行到M行的元素；步骤S2得到的内存空间Bc来提供计算矩阵C的mr0+1行至M行所需的数据，此步骤对于NN模式与NT模式，其实现流程是相同的；矩阵C的mr0+1行至M行的计算操作通过迭代来完成；ARMv8多核处理器从矩阵A中选择mr0行，对其每一行分别取N0个元素，存放到从V0开始依次至V(mr0-1)的mr0个向量寄存器中，该操作称为取数操作，该mr0行的地址是不连续的，相互之间的距离为N0×K个字节，在N0次迭代结束后再次执行该取数操作；在每次迭代中，内存空间Bc中的元素被取出并存储到从V(mp)开始依次至V(mp-1+np)的np个向量寄存器中，将矩阵A中取出的每个元素作为标量，每次迭代中从内存空间Bc中取出nr0个元素作为向量，二者依次相乘，得到N0×mp×np个元素，存储在从V(mp+np)依次至V31的向量寄存器中；在矩阵C的mr0+1行至M行的计算过程中，有部分执行数据预取操作的指令被插入计算与访存指令的执行时间间隙中；数据预取操作是指在计算指令执行前，将数据从内存预先取出到缓存中；
S4,当M％mr0不为0或N％nr0不为0时，调用ARMv8多核处理器中处理边缘数据的微内核，简称边缘微内核；对于边缘微内核的实现，其计算指令和访存指令需要按照计算需求进行排布，以确保计算指令所需的数据在执行计算指令前，已存放至向量寄存器中；计算指令FMA所需的数据提前通过取数据指令加载到向量寄存器中；
S5,在沿着K维度进行的循环迭代结束后，将ARMv8多核处理器中存储矩阵C的相关向量寄存器中的数据存回到内存中的矩阵C。
2.一种基于ARMv8多核处理器的不规则矩阵乘优化方法，其特征在于，利用ARMv8多核处理器实现,其步骤包括：
S1，对于参与相乘运算的两个矩阵，乘号左侧的矩阵称为矩阵A，乘号右侧的矩阵称为矩阵B，判断矩阵A和矩阵B是否被转置，在相乘运算中，若矩阵没被转置称为N模式，矩阵被转置称为T模式，若矩阵A为N模式，矩阵B为T模式，则将该不规则矩阵乘法称为NT模式，若矩阵A为N模式，矩阵B为N模式，则将该不规则矩阵乘法称为NN模式；建立矩阵存储空间，用于存放矩阵A与矩阵B相乘得到的结果矩阵C；
S2,确定对不规则形状矩阵乘的行方向和列方向的数据进行并行化处理所需要的并行线程数，为了保证并行线程之间的工作平衡，对矩阵C的行方向和列方向的数据分别进行并行化；将矩阵C划分为二维的子块网格，每个并行线程更新一个子块网格；当使用T个并行线程对矩阵C的计算负载进行划分时，每个并行线程将执行
Figure FDA0003336782690000041
次计算矩阵乘法的操作；每个并行线程所需的内存访问次数为
Figure FDA0003336782690000042
其中Tm和Tn分别是分配给M维和N维度数据的并行线程数，Tm×Tn＝T；对于一个子块网格的CMR，其计算公式为：
Figure FDA0003336782690000043
为了最大化CMR，使用算术-几何平均不等式，得到该CMR的取值范围为：
Figure FDA0003336782690000044
考虑S2中执行微内核打包的开销，取Tn的上界值，即
Figure FDA0003336782690000045
以最大化CMR，
Figure FDA0003336782690000046
表示向上取整；
S3，在NN模式的不规则形状矩阵乘中，如果矩阵B的大小超过最后一级缓存的容量，将t×nr0个元素打包到内存空间Bc中，其中t为自然数；打包微内核中的mr和nr的取值与主内核中对应值的取值相同；小型和不规则形状的矩阵乘的t值分别设置为0和1；
S4,当M％mr0不为0或N％nr0不为0时，调用ARMv8多核处理器中处理边缘数据的微内核，简称边缘微内核；对于边缘微内核的实现，其计算指令和访存指令需要按照计算需求进行排布，以确保计算指令所需的数据在执行计算指令前，已存放至向量寄存器中；计算指令FMA所需的数据提前通过取数据指令加载到向量寄存器中；
S5,在沿着K维度进行的循环迭代结束后，将ARMv8多核处理器中存储矩阵C的相关向量寄存器中的数据存回到内存中的矩阵C。
CN202111296475.6A 2021-11-03 2021-11-03 基于ARMv8多核处理器的小型和不规则矩阵乘优化方法 Active CN113987414B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111296475.6A CN113987414B (zh) 2021-11-03 2021-11-03 基于ARMv8多核处理器的小型和不规则矩阵乘优化方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111296475.6A CN113987414B (zh) 2021-11-03 2021-11-03 基于ARMv8多核处理器的小型和不规则矩阵乘优化方法

Publications (2)

Publication Number Publication Date
CN113987414A CN113987414A (zh) 2022-01-28
CN113987414B true CN113987414B (zh) 2022-09-09

Family

ID=79746261

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111296475.6A Active CN113987414B (zh) 2021-11-03 2021-11-03 基于ARMv8多核处理器的小型和不规则矩阵乘优化方法

Country Status (1)

Country Link
CN (1) CN113987414B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116881618A (zh) * 2023-08-25 2023-10-13 之江实验室 通用矩阵乘计算优化方法、装置及处理器
CN117556893B (zh) * 2024-01-12 2024-05-03 芯动微电子科技(武汉)有限公司 基于并行遗传算法的gpu运算gemm调优方法与装置

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102214160B (zh) * 2011-07-08 2013-04-17 中国科学技术大学 一种基于龙芯3a的单精度矩阵乘法优化方法
CN107168683B (zh) * 2017-05-05 2020-06-09 中国科学院软件研究所 申威26010众核cpu上gemm稠密矩阵乘高性能实现方法

Also Published As

Publication number Publication date
CN113987414A (zh) 2022-01-28

Similar Documents

Publication Publication Date Title
Frison et al. BLASFEO: Basic linear algebra subroutines for embedded optimization
KR102443546B1 (ko) 행렬 곱셈기
CN112445753B (zh) 从多维阵列预取多维元素块的硬件装置和方法
CN113987414B (zh) 基于ARMv8多核处理器的小型和不规则矩阵乘优化方法
EP3526665B1 (en) Sorting for data-parallel computing devices
CN111381939B (zh) 多线程处理器中的寄存器文件
Dang et al. CUDA-enabled Sparse Matrix–Vector Multiplication on GPUs using atomic operations
Krotkiewski et al. Efficient 3D stencil computations using CUDA
CN110516316B (zh) 一种间断伽辽金法求解欧拉方程的gpu加速方法
TW201737075A (zh) 複數乘法指令
Dufrechou et al. Solving sparse triangular linear systems in modern GPUs: a synchronization-free algorithm
CN114090954A (zh) 一种基于ft-2000+的整数矩阵乘法内核优化方法
Falch et al. Register caching for stencil computations on GPUs
Tang et al. Optimizing and auto-tuning iterative stencil loops for GPUs with the in-plane method
CN102722472A (zh) 一种复数矩阵的优化方法
CN101661457A (zh) 多处理器系统的三角线性方程组求解的方法和装置
Choi et al. A lightweight and efficient GPU for NDP utilizing data access pattern of image processing
Li et al. Autotsmm: An auto-tuning framework for building high-performance tall-and-skinny matrix-matrix multiplication on cpus
Cui et al. An implementation of tensor product patch smoothers on GPU
Singh et al. High-Performance Level-1 and Level-2 BLAS
CN112100099B (zh) 一种面向多核向量处理器的格子玻尔兹曼优化方法
Polok et al. Fast radix sort for sparse linear algebra on GPU
CN111857727B (zh) 一种多维循环自动向量化分块因子分块方法及装置
US20230289287A1 (en) Programmable Multi-Level Data Access Address Generator
Soliman et al. Exploiting ILP, DLP, TLP, and MPI to accelerate matrix multiplication on Xeon processors

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant