yuenshome / yuenshome.github.io

https://yuenshome.github.io

MIT License

83 stars 15 forks source link

gemm optimize #31

Open ysh329 opened 5 years ago

ysh329 commented 5 years ago

深入理解CPU和异构计算芯片GPU/FPGA/ASIC （下篇） https://zhuanlan.zhihu.com/p/25996962

面向 GEMM 引入新封装的 API | 英特尔® 软件 https://software.intel.com/zh-cn/articles/introducing-the-new-packed-apis-for-gemm?utm_source=Weibo

ysh329 commented 5 years ago

内存、cache和寄存器之间的关系及区别

寄存器是中央处理器内的组成部份。寄存器是有限存贮容量的高速存贮部件，它们可用来暂存指令、数据和位址。在中央处理器的控制部件中，包含的寄存器有指令寄存器(IR)和程序计数器(PC)。在中央处理器的算术及逻辑部件中，包含的寄存器有累加器(ACC)。
内存包含的范围非常广，一般分为只读存储器（ROM）、随机存储器（RAM）和高速缓存存储器（cache）。
寄存器是CPU内部的元件，寄存器拥有非常高的读写速度，所以在寄存器之间的数据传送非常快。
Cache ：即高速缓冲存储器，是位于CPU与主内存间的一种容量较小但速度很高的存储器。由于CPU的速度远高于主内存，CPU直接从内存中存取数据要等待一定时间周期，Cache中保存着CPU刚用过或循环使用的一部分数据，当CPU再次使用该部分数据时可从Cache中直接调用,这样就减少了CPU的等待时间,提高了系统的效率。Cache又分为一级Cache(L1 Cache)和二级Cache(L2 Cache)，L1 Cache集成在CPU内部，L2 Cache早期一般是焊在主板上,现在也都集成在CPU内部，常见的容量有256KB或512KB L2 Cache。

总结：大致来说数据是通过内存-Cache-寄存器，Cache缓存则是为了弥补CPU与内存之间运算速度的差异而设置的的部件。

计算机的存储体系（Memory hierarchy）金字塔

计算机的存储体系

ysh329 commented 5 years ago

为什么寄存器（register）比内存快？

原因一：距离不同

距离不是主要因素，但是最好懂，所以放在最前面说。内存离CPU比较远，所以要耗费更长的时间读取。以3GHz的CPU为例，电流每秒钟可以振荡30亿次，每次耗时大约为0.33纳秒。光在1纳秒的时间内，可以前进30厘米。也就是说，在CPU的一个时钟周期内，光可以前进10厘米。因此，如果内存距离CPU超过5厘米，就不可能在一个时钟周期内完成数据的读取，这还没有考虑硬件的限制和电流实际上达不到光速。相比之下，寄存器在CPU内部，当然读起来会快一点。

距离对于桌面电脑影响很大，对于手机影响就要小得多。手机CPU的时钟频率比较慢（iPhone 5s为1.3GHz），而且手机的内存紧挨着CPU。

原因二：硬件设计不同

苹果公司新推出的iPhone 5s，CPU是A7，寄存器有6000多位（31个64位寄存器，加上32个128位寄存器）。而iPhone 5s的内存是1GB，约为80亿位（bit）。这意味着，高性能、高成本、高耗电的设计可以用在寄存器上，反正只有6000多位，而不能用在内存上。因为每个位的成本和能耗只要增加一点点，就会被放大80亿倍。

事实上确实如此，内存的设计相对简单，每个位就是一个电容和一个晶体管，而寄存器的设计则完全不同，多出好几个电子元件。并且通电以后，寄存器的晶体管一直有电，而内存的晶体管只有用到的才有电，没用到的就没电，这样有利于省电。这些设计上的因素，决定了寄存器比内存读取速度更快。

原因三：工作方式不同

寄存器的工作方式很简单，只有两步：（1）找到相关的位，（2）读取这些位。

内存的工作方式就要复杂得多：

找到数据的指针。（指针可能存放在寄存器内，所以这一步就已经包括寄存器的全部工作了。）
将指针送往内存管理单元（MMU），由MMU将虚拟的内存地址翻译成实际的物理地址。
将物理地址送往内存控制器（memory controller），由内存控制器找出该地址在哪一根内存插槽（bank）上。
确定数据在哪一个内存块（chunk）上，从该块读取数据。
数据先送回内存控制器，再送回CPU，然后开始使用。

内存的工作流程比寄存器多出许多步。每一步都会产生延迟，累积起来就使得内存比寄存器慢得多。

为了缓解寄存器与内存之间的巨大速度差异，硬件设计师做出了许多努力，包括在CPU内部设置缓存、优化CPU工作方式，尽量一次性从内存读取指令所要用到的全部数据等等。

ysh329 commented 5 years ago

gemm optimize

ysh329 commented 5 years ago

1. 1x1到1x4

1.1 naive 1x1

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; ++i)
{
    for(int j = 0; j < n; ++j)
    {
        for(int p = 0; p < k; ++p)
        {
            C(i, j) += A(i, p) * B(p, j);
        }
    }
}

1.2 1x1

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; ++i)
{
    for(int j = 0; j < n; ++j)
    {
        AddDot(&A(i, 0), &B(0, j), &C(i, j), k, n);
    }
}

#define B_COL(p) b_col[p*n]
AddDot(float *a_row, float *b_col, float *c, int k, int n)
{
    for(int p = 0; p < k; ++p)
    {
        *c += a_row[p] * B_COL(p);
    }
}

1.3 1x4

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; ++i)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot(&A(i, 0), &B(0, j+0), &C(i, j+0), k, n);
        AddDot(&A(i, 0), &B(0, j+1), &C(i, j+1), k, n);
        AddDot(&A(i, 0), &B(0, j+2), &C(i, j+2), k, n);
        AddDot(&A(i, 0), &B(0, j+3), &C(i, j+3), k, n);
    }
}

#define B_COL(p) b_col[p*n]
AddDot(float *a_row, float *b_col, float *c, int k, int n)
{
    for(int p = 0; p < k; ++p)
    {
        *c += a_row[p] * B_COL(p);
    }
}

1.4 1x4

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; ++i)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot1x4(&A(i, 0), &B(0, j), &C(i, j), k, n);
    }
}

AddDot1x4(float *a_row, float *b_col, float *c, int k, int n)
{
    for(int p = 0; p < k; ++p)
        C(0, 0) += a_row[p] * B(p, 0);
    for(int p = 0; p < k; ++p)
        C(0, 1) += a_row[p] * B(p, 1);
    for(int p = 0; p < k; ++p)
        C(0, 2) += a_row[p] * B(p, 2);
    for(int p = 0; p < k; ++p)
        C(0, 3) += a_row[p] * B(p, 3);
}

1.5 1x4的提升 L2 Cache

规模m = n = k > 500左右时会有2倍的提升，因为原本的4个for循环合并为一个，一次for循环同时算4个C(i, j)：

索引p每8个浮点操作（4个乘法+加法）只需要更新一次；
元素A( 0, p )只需从内存中取出1次（不像前面取出4次），之后都在L2 cache中，直到后面的数填满了L2 cache。

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; ++i)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot1x4(&A(i, 0), &B(0, j), &C(i, j), k, n);
    }
}

AddDot1x4(float *a_row, float *b_col, float *c, int k, int n)
{
    for(int p = 0; p < k; ++p)
    {
        C(0, 0) += A(0, p) * B(p, 0);
        C(0, 1) += A(0, p) * B(p, 1);
        C(0, 2) += A(0, p) * B(p, 2);
        C(0, 3) += A(0, p) * B(p, 3);
    }
}

1.6 1x4的提升 C和A放到Register

规模m = n = k < 500时会有2倍的提升，因为把C的一行四列的1x4的计算临时结果以及A(p, 0)放到了寄存器中，减少了cache与register之间的数据传输。

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; ++i)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot1x4(&A(i, 0), &B(0, j), &C(i, j), k, n);
    }
}

AddDot1x4(float *a_row, float *b_col, float *c, int k, int n)
{
    register float
       c_00_reg,   c_01_reg,   c_02_reg,   c_03_reg,  
       a_0p_reg;

    c_00_reg = 0.0; 
    c_01_reg = 0.0; 
    c_02_reg = 0.0; 
    c_03_reg = 0.0;

    for(int p = 0; p < k; ++p)
    {
        a_0p_reg = A(0, p);

        c_00_reg += a_0p_reg * B(p, 0);
        c_01_reg += a_0p_reg  * B(p, 1);
        c_02_reg += a_0p_reg  * B(p, 2);
        c_03_reg += a_0p_reg  * B(p, 3);
    }
    C(0, 0) = c_00_reg;
    C(0, 1) = c_01_reg;
    C(0, 2) = c_02_reg;
    C(0, 3) = c_03_reg;
}

1.7 1x4的提升用指针访问B

规模m = n = k < 500时会有1.5~2倍的提升，用指针访问B，bp0_pntr, bp1_pntr, bp2_pntr, and bp3_pntr 访问 B( p, 0 ), B( p, 1 ), B( p, 2 ), B( p, 3 )，减少数据索引找B中元素的开销。

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; ++i)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot1x4(&A(i, 0), &B(0, j), &C(i, j), k, n);
    }
}

AddDot1x4(float *a_row, float *b_col, float *c, int k, int n)
{
    register float
       c_00_reg,   c_01_reg,   c_02_reg,   c_03_reg,  
       a_0p_reg;

    c_00_reg = 0.0; 
    c_01_reg = 0.0; 
    c_02_reg = 0.0; 
    c_03_reg = 0.0;

    float
    *bp0_pntr, *bp1_pntr, *bp2_pntr, *bp3_pntr; 

    bp0_pntr = &B( 0, 0 );
    bp1_pntr = &B( 0, 1 );
    bp2_pntr = &B( 0, 2 );
    bp3_pntr = &B( 0, 3 );

    for(int p = 0; p < k; ++p)
    {
        a_0p_reg = A(0, p);

        c_00_reg += a_0p_reg * *bp0_pntr++;
        c_01_reg += a_0p_reg * *bp1_pntr++;
        c_02_reg += a_0p_reg * *bp2_pntr++;
        c_03_reg += a_0p_reg * *bp3_pntr++;
    }
    C(0, 0) = c_00_reg;
    C(0, 1) = c_01_reg;
    C(0, 2) = c_02_reg;
    C(0, 3) = c_03_reg;
}

1.8 1x4的轻微下降对k方向展开p+=4，这个累加的数也可以自定

unroll这个for循环，从原本的+=1改为+=4混淆了编译器（unroll多少并非固定，可以根据不同的平台调节），反而有性能下降。

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; ++i)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot1x4(&A(i, 0), &B(0, j), &C(i, j), k, n);
    }
}

AddDot1x4(float *a_row, float *b_col, float *c, int k, int n)
{
    register float
       c_00_reg,   c_01_reg,   c_02_reg,   c_03_reg,  
       a_0p_reg;

    c_00_reg = 0.0; 
    c_01_reg = 0.0; 
    c_02_reg = 0.0; 
    c_03_reg = 0.0;

    float
    *bp0_pntr, *bp1_pntr, *bp2_pntr, *bp3_pntr; 

    bp0_pntr = &B( 0, 0 );
    bp1_pntr = &B( 0, 1 );
    bp2_pntr = &B( 0, 2 );
    bp3_pntr = &B( 0, 3 );

    for(int p = 0; p < k; p+=4)
    {
        a_0p_reg = A(0, p);
        c_00_reg += a_0p_reg * *bp0_pntr++;
        c_01_reg += a_0p_reg * *bp1_pntr++;
        c_02_reg += a_0p_reg * *bp2_pntr++;
        c_03_reg += a_0p_reg * *bp3_pntr++;

        a_0p_reg = A(0, p+1);
        c_00_reg += a_0p_reg * *bp0_pntr++;
        c_01_reg += a_0p_reg * *bp1_pntr++;
        c_02_reg += a_0p_reg * *bp2_pntr++;
        c_03_reg += a_0p_reg * *bp3_pntr++;

        a_0p_reg = A(0, p+2);
        c_00_reg += a_0p_reg * *bp0_pntr++;
        c_01_reg += a_0p_reg * *bp1_pntr++;
        c_02_reg += a_0p_reg * *bp2_pntr++;
        c_03_reg += a_0p_reg * *bp3_pntr++;

        a_0p_reg = A(0, p+3);
        c_00_reg += a_0p_reg * *bp0_pntr++;
        c_01_reg += a_0p_reg * *bp1_pntr++;
        c_02_reg += a_0p_reg * *bp2_pntr++;
        c_03_reg += a_0p_reg * *bp3_pntr++;
    }
    C(0, 0) = c_00_reg;
    C(0, 1) = c_01_reg;
    C(0, 2) = c_02_reg;
    C(0, 3) = c_03_reg;
}

1.9 1x4 间接寻址（indirect addressing）

间接寻址操作如下：

c_00_reg += a_0p_reg * *(bp0_pntr+1);

bp0_pntr指向元素B( p, 0 )。因此bp0_pntr+1指向元素B( p+1, 0 )，间接访问这种方法不需要更新指针的地址，而是直接通过+1访问下一个元素。这样，每算4次，b指针才会被更新，通过间接寻址减少了指针变量的更新次数。

注：没有性能提升，因为编译器已经做了这样的优化，我们只是显式地写了出来。

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; ++i)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot1x4(&A(i, 0), &B(0, j), &C(i, j), k, n);
    }
}

AddDot1x4(float *a_row, float *b_col, float *c, int k, int n)
{
    register float
       c_00_reg,   c_01_reg,   c_02_reg,   c_03_reg,  
       a_0p_reg;

    c_00_reg = 0.0; 
    c_01_reg = 0.0; 
    c_02_reg = 0.0; 
    c_03_reg = 0.0;

    float
    *bp0_pntr, *bp1_pntr, *bp2_pntr, *bp3_pntr; 

    bp0_pntr = &B( 0, 0 );
    bp1_pntr = &B( 0, 1 );
    bp2_pntr = &B( 0, 2 );
    bp3_pntr = &B( 0, 3 );

    for(int p = 0; p < k; p+=4)
    {
        a_0p_reg = A(0, p);
        c_00_reg += a_0p_reg * *bp0_pntr;
        c_01_reg += a_0p_reg * *bp1_pntr;
        c_02_reg += a_0p_reg * *bp2_pntr;
        c_03_reg += a_0p_reg * *bp3_pntr;

        a_0p_reg = A(0, p+1);
        c_00_reg += a_0p_reg * *(bp0_pntr+1);
        c_01_reg += a_0p_reg * *(bp1_pntr+1);
        c_02_reg += a_0p_reg * *(bp2_pntr+1);
        c_03_reg += a_0p_reg * *(bp3_pntr+1);

        a_0p_reg = A(0, p+2);
        c_00_reg += a_0p_reg * *(bp0_pntr+2);
        c_01_reg += a_0p_reg * *(bp1_pntr+2);
        c_02_reg += a_0p_reg * *(bp2_pntr+2);
        c_03_reg += a_0p_reg * *(bp3_pntr+2);

        a_0p_reg = A(0, p+3);
        c_00_reg += a_0p_reg * *(bp0_pntr+3);
        c_01_reg += a_0p_reg * *(bp1_pntr+3);
        c_02_reg += a_0p_reg * *(bp2_pntr+3);
        c_03_reg += a_0p_reg * *(bp3_pntr+3);

        bp0_pntr+=4;
        bp1_pntr+=4;
        bp2_pntr+=4;
        bp3_pntr+=4;
    }
    C(0, 0) = c_00_reg;
    C(0, 1) = c_01_reg;
    C(0, 2) = c_02_reg;
    C(0, 3) = c_03_reg;
}

ysh329 commented 5 years ago

2. 1x4到4x4

后面使用到向量指令以及向量寄存器，做到一次算4x4个结果。

有一部分特殊指令集如SSE3，一个周期内可执行两个乘加操作（两个乘法和两个加法操作），也就是每个周期可对4个浮点操作；
要把数据放到向量寄存器（vector register）才能实现如此操作；

因为有16个向量寄存器，每个都可以存放两个双精度double或四个单精度的浮点数。因而可以将32个双精度浮数点或64个单精度浮点数放到这16个向量寄存器中，若是双精度即我们4x4大小的矩阵C的元素，若是单精度即可以放8x8规模大小的矩阵C的元素。

2.1 4x4 A的行 B的列每次取4

扩展AddDot，每次取A的4行，B的4列，在AddDot4x4中计算16个C元素。

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; i+=4)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot4x4(&A(i, 0), &B(0, j), &C(i, j), n, k);
    }
}

void AddDot4x4(float *a, float *b, float *c, const int n, const int k)
{
  /* So, this routine computes a 4x4 block of matrix A
           C( 0, 0 ), C( 0, 1 ), C( 0, 2 ), C( 0, 3 ).  
           C( 1, 0 ), C( 1, 1 ), C( 1, 2 ), C( 1, 3 ).  
           C( 2, 0 ), C( 2, 1 ), C( 2, 2 ), C( 2, 3 ).  
           C( 3, 0 ), C( 3, 1 ), C( 3, 2 ), C( 3, 3 ).  
     Notice that this routine is called with c = C( i, j ) in the
     previous routine, so these are actually the elements 
           C( i  , j ), C( i  , j+1 ), C( i  , j+2 ), C( i  , j+3 ) 
           C( i+1, j ), C( i+1, j+1 ), C( i+1, j+2 ), C( i+1, j+3 ) 
           C( i+2, j ), C( i+2, j+1 ), C( i+2, j+2 ), C( i+2, j+3 ) 
           C( i+3, j ), C( i+3, j+1 ), C( i+3, j+2 ), C( i+3, j+3 ) 

     in the original matrix C 
     In this version, we "inline" AddDot */ 

    // 1st row of C
    AddDot(&A(0, 0), &B(0, 0), &C(0, 0), k, n);
    AddDot(&A(0, 0), &B(0, 1), &C(0, 1), k, n);
    AddDot(&A(0, 0), &B(0, 2), &C(0, 2), k, n);
    AddDot(&A(0, 0), &B(0, 3), &C(0, 3), k, n);

    // 2rd row of C
    AddDot(&A(1, 0), &B(0, 0), &C(1, 0), k, n);
    AddDot(&A(1, 0), &B(0, 1), &C(1, 1), k, n);
    AddDot(&A(1, 0), &B(0, 2), &C(1, 2), k, n);
    AddDot(&A(1, 0), &B(0, 3), &C(1, 3), k, n);

    // 3nd row of C
    AddDot(&A(2, 0), &B(0, 0), &C(2, 0), k, n);
    AddDot(&A(2, 0), &B(0, 1), &C(2, 1), k, n);
    AddDot(&A(2, 0), &B(0, 2), &C(2, 2), k, n);
    AddDot(&A(2, 0), &B(0, 3), &C(2, 3), k, n);

    // 4th row of C
    AddDot(&A(3, 0), &B(0, 0), &C(3, 0), k, n);
    AddDot(&A(3, 0), &B(0, 1), &C(3, 1), k, n);
    AddDot(&A(3, 0), &B(0, 2), &C(3, 2), k, n);
    AddDot(&A(3, 0), &B(0, 3), &C(3, 3), k, n);
}

#define B_COL(p) b_col[p*n]
void AddDot(float *a_row, float *b_col, float *c, int k, int n)
{
    for(int p = 0; p < k; ++p)
    {
        *c += a_row[p] * B_COL(p);
    }
}

2.2 4x4 拆分AddDot到AddDot4x4中

无性能提升。

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; i+=4)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot4x4(&A(i, 0), &B(0, j), &C(i, j), n, k);
    }
}

void AddDot4x4(float *a, float *b, float *c, const int n, const int k)
{
    // 1st row of C
    for(int p = 0; p < k; ++p)
        C(0, 0) += A(0, p) * B(p, 0);
    for(int p = 0; p < k; ++p)
        C(0, 1) += A(0, p) * B(p, 1);
    for(int p = 0; p < k; ++p)
        C(0, 2) += A(0, p) * B(p, 2);
    for(int p = 0; p < k; ++p)
        C(0, 3) += A(0, p) * B(p, 3);

    // 2rd row of C
    for(int p = 0; p < k; ++p)
        C(1, 0) += A(1, p) * B(p, 0);
    for(int p = 0; p < k; ++p)
        C(1, 1) += A(1, p) * B(p, 1);
    for(int p = 0; p < k; ++p)
        C(1, 2) += A(1, p) * B(p, 2);
    for(int p = 0; p < k; ++p)
        C(1, 3) += A(1, p) * B(p, 3);

    // 3nd row of C
    for(int p = 0; p < k; ++p)
        C(2, 0) += A(2, p) * B(p, 0);
    for(int p = 0; p < k; ++p)
        C(2, 1) += A(2, p) * B(p, 1);
    for(int p = 0; p < k; ++p)
        C(2, 2) += A(2, p) * B(p, 2);
    for(int p = 0; p < k; ++p)
        C(2, 3) += A(2, p) * B(p, 3);

    // 4th row of C
    for(int p = 0; p < k; ++p)
        C(3, 0) += A(3, p) * B(p, 0);
    for(int p = 0; p < k; ++p)
        C(3, 1) += A(3, p) * B(p, 1);
    for(int p = 0; p < k; ++p)
        C(3, 2) += A(3, p) * B(p, 2);
    for(int p = 0; p < k; ++p)
        C(3, 3) += A(3, p) * B(p, 3);
}

2.3 4x4 提升合并AddDot4x4中的For循环

在m=n=k>500时，4~5倍的性能提升，合并16个For循环为1个。

矩阵变更大，不少数据在寄存器中可被复用地更多。

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; i+=4)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot4x4(&A(i, 0), &B(0, j), &C(i, j), n, k);
    }
}

void AddDot4x4(float *a, float *b, float *c, const int n, const int k)
{
    for(int p = 0; p < k; ++p)
    {
        // 1st row of C
        C(0, 0) += A(0, p) * B(p, 0);
        C(0, 1) += A(0, p) * B(p, 1);
        C(0, 2) += A(0, p) * B(p, 2);
        C(0, 3) += A(0, p) * B(p, 3);

        // 2rd row of C
        C(1, 0) += A(1, p) * B(p, 0);
        C(1, 1) += A(1, p) * B(p, 1);
        C(1, 2) += A(1, p) * B(p, 2);
        C(1, 3) += A(1, p) * B(p, 3);

        // 3nd row of C
        C(2, 0) += A(2, p) * B(p, 0);
        C(2, 1) += A(2, p) * B(p, 1);
        C(2, 2) += A(2, p) * B(p, 2);
        C(2, 3) += A(2, p) * B(p, 3);

        // 4th row of C
        C(3, 0) += A(3, p) * B(p, 0);
        C(3, 1) += A(3, p) * B(p, 1);
        C(3, 2) += A(3, p) * B(p, 2);
        C(3, 3) += A(3, p) * B(p, 3);
    }
}

2.4 4x4 提升 C和A放到Register中

m=n=k<500，性能有2倍提升；
m=n=k>500，性能有1.5倍提升。

得益于将A的4行1列和C的4行4列放到寄存器中。当使用超过寄存器容量的寄存器变量时，性能提升就固定了（即m=n=k>500，性能有1.5倍提升）。

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; i+=4)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot4x4(&A(i, 0), &B(0, j), &C(i, j), n, k);
    }
}

void AddDot4x4(float *a, float *b, float *c, const int n, const int k)
{
  register float
    /* hold contributions to
       C( 0, 0 ), C( 0, 1 ), C( 0, 2 ), C( 0, 3 ) 
       C( 1, 0 ), C( 1, 1 ), C( 1, 2 ), C( 1, 3 ) 
       C( 2, 0 ), C( 2, 1 ), C( 2, 2 ), C( 2, 3 ) 
       C( 3, 0 ), C( 3, 1 ), C( 3, 2 ), C( 3, 3 )   */

       c_00_reg,   c_01_reg,   c_02_reg,   c_03_reg,  
       c_10_reg,   c_11_reg,   c_12_reg,   c_13_reg,  
       c_20_reg,   c_21_reg,   c_22_reg,   c_23_reg,  
       c_30_reg,   c_31_reg,   c_32_reg,   c_33_reg,

    /* hold 
       A( 0, p ) 
       A( 1, p ) 
       A( 2, p ) 
       A( 3, p ) */

       a_0p_reg,
       a_1p_reg,
       a_2p_reg,
       a_3p_reg;

    c_00_reg = 0.0;   c_01_reg = 0.0;   c_02_reg = 0.0;   c_03_reg = 0.0;
    c_10_reg = 0.0;   c_11_reg = 0.0;   c_12_reg = 0.0;   c_13_reg = 0.0;
    c_20_reg = 0.0;   c_21_reg = 0.0;   c_22_reg = 0.0;   c_23_reg = 0.0;
    c_30_reg = 0.0;   c_31_reg = 0.0;   c_32_reg = 0.0;   c_33_reg = 0.0;

    for(int p = 0; p < k; ++p)
    {
        a_0p_reg = A(0, p);
        a_1p_reg = A(1, p);
        a_2p_reg = A(2, p);
        a_3p_reg = A(3, p);

        // 1st row of C
        c_00_reg += a_0p_reg * B( p, 0 );
        c_01_reg += a_0p_reg * B( p, 1 );
        c_02_reg += a_0p_reg * B( p, 2 );
        c_03_reg += a_0p_reg * B( p, 3 );

        // 2rd row of C
        c_10_reg += a_1p_reg * B( p, 0 );
        c_11_reg += a_1p_reg * B( p, 1 );
        c_12_reg += a_1p_reg * B( p, 2 );
        c_13_reg += a_1p_reg * B( p, 3 );

        // 3nd row of C
        c_20_reg += a_2p_reg * B( p, 0 );
        c_21_reg += a_2p_reg * B( p, 1 );
        c_22_reg += a_2p_reg * B( p, 2 );
        c_23_reg += a_2p_reg * B( p, 3 );

        // 4th row of C
        c_30_reg += a_3p_reg * B( p, 0 );
        c_31_reg += a_3p_reg * B( p, 1 );
        c_32_reg += a_3p_reg * B( p, 2 );
        c_33_reg += a_3p_reg * B( p, 3 );
    }
    C( 0, 0 ) += c_00_reg;   C( 0, 1 ) += c_01_reg;   C( 0, 2 ) += c_02_reg;   C( 0, 3 ) += c_03_reg;
    C( 1, 0 ) += c_10_reg;   C( 1, 1 ) += c_11_reg;   C( 1, 2 ) += c_12_reg;   C( 1, 3 ) += c_13_reg;
    C( 2, 0 ) += c_20_reg;   C( 2, 1 ) += c_21_reg;   C( 2, 2 ) += c_22_reg;   C( 2, 3 ) += c_23_reg;
    C( 3, 0 ) += c_30_reg;   C( 3, 1 ) += c_31_reg;   C( 3, 2 ) += c_32_reg;   C( 3, 3 ) += c_33_reg;
}

2.5 4x4 提升用指针代替循环中的B

轻微性能提升，用指针代替访问B，减少索引开销。

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; i+=4)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot4x4(&A(i, 0), &B(0, j), &C(i, j), n, k);
    }
}

void AddDot4x4(float *a, float *b, float *c, const int n, const int k)
{
    register float 
    c_00_reg,   c_01_reg,   c_02_reg,   c_03_reg,  
    c_10_reg,   c_11_reg,   c_12_reg,   c_13_reg,  
    c_20_reg,   c_21_reg,   c_22_reg,   c_23_reg,  
    c_30_reg,   c_31_reg,   c_32_reg,   c_33_reg,

    a_0p_reg,
    a_1p_reg,
    a_2p_reg,
    a_3p_reg;

    float
    *b_p0_pntr, *b_p1_pntr, *b_p2_pntr, *b_p3_pntr; 

    c_00_reg = 0.0;   c_01_reg = 0.0;   c_02_reg = 0.0;   c_03_reg = 0.0;
    c_10_reg = 0.0;   c_11_reg = 0.0;   c_12_reg = 0.0;   c_13_reg = 0.0;
    c_20_reg = 0.0;   c_21_reg = 0.0;   c_22_reg = 0.0;   c_23_reg = 0.0;
    c_30_reg = 0.0;   c_31_reg = 0.0;   c_32_reg = 0.0;   c_33_reg = 0.0;

    for(int p = 0; p < k; ++p)
    {
        a_0p_reg = A(0, p);
        a_1p_reg = A(1, p);
        a_2p_reg = A(2, p);
        a_3p_reg = A(3, p);

        b_p0_pntr = &B( p, 0 );
        b_p1_pntr = &B( p, 1 );
        b_p2_pntr = &B( p, 2 );
        b_p3_pntr = &B( p, 3 );

        // 1st row of C
        c_00_reg += a_0p_reg * *b_p0_pntr;
        c_01_reg += a_0p_reg * *b_p1_pntr;
        c_02_reg += a_0p_reg * *b_p2_pntr;
        c_03_reg += a_0p_reg * *b_p3_pntr;

        // 2rd row of C
        c_10_reg += a_1p_reg * *b_p0_pntr;
        c_11_reg += a_1p_reg * *b_p1_pntr;
        c_12_reg += a_1p_reg * *b_p2_pntr;
        c_13_reg += a_1p_reg * *b_p3_pntr;

        // 3nd row of C
        c_20_reg += a_2p_reg * *b_p0_pntr;
        c_21_reg += a_2p_reg * *b_p1_pntr;
        c_22_reg += a_2p_reg * *b_p2_pntr;
        c_23_reg += a_2p_reg * *b_p3_pntr;

        // 4th row of C
        c_30_reg += a_3p_reg * *b_p0_pntr++;
        c_31_reg += a_3p_reg * *b_p1_pntr++;
        c_32_reg += a_3p_reg * *b_p2_pntr++;
        c_33_reg += a_3p_reg * *b_p3_pntr++;
    }

    C( 0, 0 ) += c_00_reg;   C( 0, 1 ) += c_01_reg;   C( 0, 2 ) += c_02_reg;   C( 0, 3 ) += c_03_reg;
    C( 1, 0 ) += c_10_reg;   C( 1, 1 ) += c_11_reg;   C( 1, 2 ) += c_12_reg;   C( 1, 3 ) += c_13_reg;
    C( 2, 0 ) += c_20_reg;   C( 2, 1 ) += c_21_reg;   C( 2, 2 ) += c_22_reg;   C( 2, 3 ) += c_23_reg;
    C( 3, 0 ) += c_30_reg;   C( 3, 1 ) += c_31_reg;   C( 3, 2 ) += c_32_reg;   C( 3, 3 ) += c_33_reg;
}

ysh329 commented 5 years ago

3. 进一步优化4x4（Further optimizing）

3.1 4x4 下降将循环中的B用指针索引后，放到寄存器

将从B的第j行的指针中，取出的变量放到寄存器中。这样做性能些许下降，但是为了后面能作进一步的优化。

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; i+=4)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot4x4(&A(i, 0), &B(0, j), &C(i, j), n, k);
    }
}

void AddDot4x4(float *a, float *b, float *c, const int n, const int k)
{
    register float
    c_00_reg,   c_01_reg,   c_02_reg,   c_03_reg,  
    c_10_reg,   c_11_reg,   c_12_reg,   c_13_reg,  
    c_20_reg,   c_21_reg,   c_22_reg,   c_23_reg,  
    c_30_reg,   c_31_reg,   c_32_reg,   c_33_reg,

    a_0p_reg,
    a_1p_reg,
    a_2p_reg,
    a_3p_reg,

    b_p0_reg,
    b_p1_reg,
    b_p2_reg,
    b_p3_reg;

    float
    *b_p0_pntr, *b_p1_pntr, *b_p2_pntr, *b_p3_pntr; 
    b_p0_pntr = &B( 0, 0 );
    b_p1_pntr = &B( 0, 1 );
    b_p2_pntr = &B( 0, 2 );
    b_p3_pntr = &B( 0, 3 );

    c_00_reg = 0.0;   c_01_reg = 0.0;   c_02_reg = 0.0;   c_03_reg = 0.0;
    c_10_reg = 0.0;   c_11_reg = 0.0;   c_12_reg = 0.0;   c_13_reg = 0.0;
    c_20_reg = 0.0;   c_21_reg = 0.0;   c_22_reg = 0.0;   c_23_reg = 0.0;
    c_30_reg = 0.0;   c_31_reg = 0.0;   c_32_reg = 0.0;   c_33_reg = 0.0;

    for(int p = 0; p < k; ++p)
    {
        a_0p_reg = A(0, p);
        a_1p_reg = A(1, p);
        a_2p_reg = A(2, p);
        a_3p_reg = A(3, p);

        b_p0_reg = *b_p0_pntr++;
        b_p1_reg = *b_p1_pntr++;
        b_p2_reg = *b_p2_pntr++;
        b_p3_reg = *b_p3_pntr++;

        // 1st row of C
        c_00_reg += a_0p_reg * b_p0_reg;
        c_01_reg += a_0p_reg * b_p1_reg;
        c_02_reg += a_0p_reg * b_p2_reg;
        c_03_reg += a_0p_reg * b_p3_reg;

        // 2nd row of C
        c_10_reg += a_1p_reg * b_p0_reg;
        c_11_reg += a_1p_reg * b_p1_reg;
        c_12_reg += a_1p_reg * b_p2_reg;
        c_13_reg += a_1p_reg * b_p3_reg;

        // 3rd row of C
        c_20_reg += a_2p_reg * b_p0_reg;
        c_21_reg += a_2p_reg * b_p1_reg;
        c_22_reg += a_2p_reg * b_p2_reg;
        c_23_reg += a_2p_reg * b_p3_reg;

        // 4th row of C
        c_30_reg += a_3p_reg * b_p0_reg;
        c_31_reg += a_3p_reg * b_p1_reg;
        c_32_reg += a_3p_reg * b_p2_reg;
        c_33_reg += a_3p_reg * b_p3_reg;
    }

    C( 0, 0 ) += c_00_reg;   C( 0, 1 ) += c_01_reg;   C( 0, 2 ) += c_02_reg;   C( 0, 3 ) += c_03_reg;
    C( 1, 0 ) += c_10_reg;   C( 1, 1 ) += c_11_reg;   C( 1, 2 ) += c_12_reg;   C( 1, 3 ) += c_13_reg;
    C( 2, 0 ) += c_20_reg;   C( 2, 1 ) += c_21_reg;   C( 2, 2 ) += c_22_reg;   C( 2, 3 ) += c_23_reg;
    C( 3, 0 ) += c_30_reg;   C( 3, 1 ) += c_31_reg;   C( 3, 2 ) += c_32_reg;   C( 3, 3 ) += c_33_reg;
}

3.2 4x4 提升变为vector operation进一步转为intrinsic

规模m=n=k<500下，基本都有2~3倍的性能提升（但规模在m=n=k>500下，则几乎没有提升，后面将讲到如何优化使用L2 Cache来做更大规模的提升），提升的原因是用上了向量寄存器和向量操作（或者说是SIMD）。这个实现过程，先考虑后续intrinsic的写法，即先写成对应的不用intrinsic的代码，然后再转成intrinsic函数。

实现过程中参考这个double的4x4的做法：Optimization_4x4_9 · flame/how-to-optimize-gemm Wiki，但是这种情况与我们的不同：1. 咱们是float；2. 咱们是行优先。

因此，对B按行取4个的前提下，对于矩阵A在取数时有如下三个思路：

对B按行取4个的前提下，有三个对A处理的思路：

每次取A的第1列的4个，与B的第1行的4个，即两个float4。因行优先，对B来说是连续的，对A则需要逐个元素竖着排布拼成float4；
取A行的第1个元素，重复4次，即float4，这样取A的四行，得到4个float4；
对B转置为B_trans，与A操作一样，按照行取4行（实际为B的4列）；或相反，对A转置为A_trans，则与B操作一样。

我实现的是第二种，即对A中的元素进行重复，如下是实现代码：

#define A(i, j) a[k*(i) + (j)]
#define B(i, j) a[n*(i) + (j)]
#define C(i, j) a[n*(i) + (j)]

for(int i = 0; i < m; i+=4)
{
    for(int j = 0; j < n; j+=4)
    {
        AddDot4x4(&A(i, 0), &B(0, j), &C(i, j), n, k);
    }
}

#include <mmintrin.h>
#include <xmmintrin.h>  // SSE
#include <pmmintrin.h>  // SSE2
#include <emmintrin.h>  // SSE3

typedef union v4f
{
  __m128 v;
  float s[4];
} v4f_t;

void AddDot4x4(float *a, float *b, float *c, const int n, const int k)
{   
    v4f_t
      c_00_01_02_03_vreg,
      c_10_11_12_13_vreg,
      c_20_21_22_23_vreg,
      c_30_31_32_33_vreg,

      b_p0_p1_p2_p3_vreg,

      a_0p_x4_vreg,
      a_1p_x4_vreg,
      a_2p_x4_vreg,
      a_3p_x4_vreg;

    float 
      *b_p0_pntr = &B( 0, 0 );

    c_00_01_02_03_vreg.v = _mm_setzero_ps();
    c_10_11_12_13_vreg.v = _mm_setzero_ps();
    c_20_21_22_23_vreg.v = _mm_setzero_ps();
    c_30_31_32_33_vreg.v = _mm_setzero_ps();

    for(int p = 0; p < k; ++p)
    {
        // duplicate each element in A as float4
        a_0p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 0, p ) );
        a_1p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 1, p ) );
        a_2p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 2, p ) );
        a_3p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 3, p ) );

        b_p0_pntr += p * n;
        b_p0_p1_p2_p3_vreg.v = __mm_load_ps( (float *) b_p0_pntr );

        c_00_01_02_03_vreg.v += a_0p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
        c_10_11_12_13_vreg.v += a_1p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
        c_20_21_22_23_vreg.v += a_2p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
        c_30_31_32_33_vreg.v += a_3p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
    }

    // 1st row of C
    C( 0, 0 ) += c_00_01_02_03_vreg.s[0];
    C( 0, 1 ) += c_00_01_02_03_vreg.s[1];
    C( 0, 2 ) += c_00_01_02_03_vreg.s[2];
    C( 0, 3 ) += c_00_01_02_03_vreg.s[3];
    // 2nd row of C
    C( 1, 0 ) += c_10_11_12_13_vreg.s[0];
    C( 1, 1 ) += c_10_11_12_13_vreg.s[1];
    C( 1, 2 ) += c_10_11_12_13_vreg.s[2];
    C( 1, 3 ) += c_10_11_12_13_vreg.s[3];
    // 3rd row of C
    C( 2, 0 ) += c_20_11_12_13_vreg.s[0];
    C( 2, 1 ) += c_20_11_12_13_vreg.s[1];
    C( 2, 2 ) += c_20_11_12_13_vreg.s[2];
    C( 2, 3 ) += c_20_11_12_13_vreg.s[3];
    // 4th row of C
    C( 3, 0 ) += c_30_31_32_33_vreg.s[0];
    C( 3, 1 ) += c_30_31_32_33_vreg.s[1];
    C( 3, 2 ) += c_30_31_32_33_vreg.s[2];
    C( 3, 3 ) += c_30_31_32_33_vreg.s[3];
}

ysh329 commented 5 years ago

英特尔® 高级矢量扩展指令集简介

Optimization Notice

虽然针对非英特尔处理器也有优化，但不保证有效，具体参考用户和参考指南，其中包括了本注意事项涵盖的特定指令集的更多相关信息。

英特尔® 高级矢量扩展指令集（英特尔® AVX）是在英特尔® 架构 CPU 上执行单指令多数据 (SIMD) 运算的指令集。这些指令添加了以下特性，对之前的 SIMD 产品——MMX™ 指令和英特尔® 数据流单指令多数据扩展指令集（英特尔® SSE）进行了扩展：

将 128 位 SIMD 寄存器扩展至 256 位。英特尔® AVX 的目标是在未来可以支持 512 或 1024 位。
添加了 3 操作数非破坏性运算。之前在 A = A + B 类运算中执行的是 2 操作数指令，它将覆盖源操作数，而新的操作数可以执行 A = B + C 类运算，且保持原始源操作数不变。.
少数几个指令采用 4 寄存器操作数，通过移除不必要的指令，支持更小、更快的代码。
对于操作数的内存对齐要求有所放宽。
新的扩展编码方案 (VEX) 旨在使得以后添加更加容易以及所执行的指令编码更小、速度更快。

与这些改进密切相关的是新的*融合乘加 (FMA) 指令，它可以更快、更精确地执行专门运算，例如单指令 A = A B + C**。第二代英特尔® 酷睿™ CPU 中将可提供 FMA 指令。其他特性还包括处理高级加密标准 (AES) 加密和解密的指令、用于某些加密基元的紧缩 carry-less 乘法运算 (PCLMULQDQ) 以及某些用于未来指令的预留槽，如硬件随机数生成器。

指令集概述

新指令使用英特尔的 VEX prefix 进行编码，VEX prefix 是一个 2个或 3 个字节的前缀，旨在降低当前和未来 x86/x64 指令编码的复杂性。两个新的 VEX prefix 形成于两个过时的 32 位指令——使用 DS 的 Load Pointer（LDS-0xC4，3 字节形式）和使用 ES 的 Load Pointer（LES-0xC5， 2 字节形式），它们以 32 位模式加载 DS 和 ES 段寄存器。在 64 位模式中，操作码 LDS 和 LES 生成一个无效操作码异常，但是在英特尔® AVX 下，这些操作码可另外用作编码新指令前缀。最终，VEX 指令只能在 64 位模式中运行时使用。前缀可比之前的 x86 指令编码更多的寄存器，可用于访问新的 256 位 SIMD 寄存器或者使用 3 和 4 操作数语法。作为用户，您无需担心这个问题（除非您正在编写汇编器或反汇编器）。

注：下文假定运算都是在 64 位模式中进行。

SIMD 指令可以在一个单一步骤中处理多个片段的数据，加速许多任务的吞吐量，包括视频编解码、图像处理、数据分析和物理模拟。在 32 位长度（称之为单精度）和64 位长度（称之为双精度）中，英特尔® AVX 指令在IEEE-754浮点值上运行。IEEE-754 是一个标准定义的、可复制的强大浮点运算，是大多数主流数字计算的标准。

较早的相关英特尔® SSE 指令还支持各种带符号和无符号整数大小，包括带符号和无符号字节（B，8 位）、字（W，16 位）、双字（DW，32 位）、四字（QW，64 位）和双四字（DQ，128 位）长度。并非所有指令都适用于所有大小组合，更多详细信息，敬请参阅“更多信息”中提供的链接。请参阅本文后文中的图 2，了解数据类型的图形表示法。

支持英特尔® AVX（和 FMA）的硬件包括 16 个 256 位 YMM 寄存器 YMM0-YMM15 和一个名为 MXCSR的 32 位控制/状态寄存器。YMM 寄存器替代了英特尔® SSE 使用的较早的 128 位 XMM 寄存器，它将 XMM 寄存器视作相应 YMM 寄存器的下层部分，如图 1 所示。

图 1. XMM 寄存器覆盖 YMM 寄存器。

MXCSR 的 0-5 位显示设置“粘滞”位后 SIMD 浮点异常，除非使用位LDMXCSR 或 FXRSTOR清除，否则它们将一直保持设置。设置时 7-12 位屏蔽个体异常，可通过启动进行初始设置或重置。0-5 位分别显示为无效运算、非法、被零除、溢出和精度。如欲获取详细信息，请参阅“更多信息”部分提供的链接。

图 2 显示了英特尔® SSE 和英特尔® AVX 指令中使用的数据类型。对于英特尔® AVX，允许使用增加至 128 或 256 位的 32 位或 64 位浮点类型的倍数以及增加至 128 位的任意整数类型的倍数。

图 2.英特尔® AVX 和英特尔® SSE 数据类型

指令通常分为标量版本和矢量版本，如图 3 所示。矢量版本通过将数据以并行“SIMD”模式在寄存器中处理进行运算；而标量版本则只在每个寄存器的一个条目中进行运算。这一区别减少了某些算法中的数据移动，提供了更加出色的整体吞吐量。

图 3. SIMD 与标量运算

内存对齐

当数据以 n 字节内存界限上存储的 n 字节块进行运算时，数据为内存对齐数据。例如，将 256 位数据加载至 YMM 寄存器中时，如果数据源为 256 位对齐，则数据被称为“对齐”.

对于英特尔® SSE 运算，除非明确规定，否则都需要内存对齐。例如，在英特尔® SSE 下，有针对内存对齐和内存未对齐运算的特定指令，如MOVAPD（移动对齐的紧缩双精度值）和 MOVUPD（移动非对齐的紧缩双精度值）指令。没有像这样分为两个的指令需要执行对齐访问。

英特尔® AVX 已经放宽了某些内存对齐要求，因此默认情况下，英特尔® AVX 允许未对齐的访问；但是，这样的访问将导致性能下降，因此旨在要求数据保持内存对齐的原规则仍然是不错的选择（面向 128 位访问的 16 位对齐和面向 256 位访问的 32 位对齐）。主要例外是明确指出需要内存对齐数据的 SSE 指令的VEX 扩展版本：这些指令仍然要求使用对齐的数据。其他要求对齐访问的特定指令请参阅英特尔® 高级矢量扩展指令集编程指南中的表2.4（请参阅“更多信息”中提供的链接）。

除了未对齐数据问题外，另外一个性能问题是混合使用旧的仅 XMM 的指令和较新的英特尔® AVX 指令会导致延迟，因此需要最大限度地控制 VEX 编码指令和旧的英特尔® SSE 代码之间的转换。也就是说，不要将 VEX 前缀的指令和非 VEX 前缀的指令混合使用，以实现最佳吞吐量。如果您非要这么做，请对同一 VEX/非 VEX 级别的指令进行分组，最大限度地控制二者之间的转换。此外，如果上层 YMM 位通过 VZEROUPPER 或 VZEROALL 设置为零（编译器应自动插入），则无转换损失。该插入要求另外一个指令，因此推荐使用分析 (profiling)。

ysh329 commented 5 years ago

4. 4x4 提升利用L2 Cache对矩阵blocking

4.1 参考代码分析

参考：Optimization_4x4_11 · flame/how-to-optimize-gemm Wiki

对于参考代码的分析，参考代码更新后，规模在m=n=k>500有2倍率性能提升，对矩阵C做blocking（每次内层循环InnerKernel计算一块规模为mc x n大小的矩阵C，但并没有完全算完，只算了pb/k * ib/m，其中pb = min( k-p, kc ), ib = min( m-i, mc )）。对矩阵的分块充分利用L2 Cache。


/* Block sizes */
#define mc 256
#define kc 128

  /* This time, we compute a mc x n block of C by a call to the InnerKernel */
  for ( p=0; p<k; p+=kc ){
    pb = min( k-p, kc );
    for ( i=0; i<m; i+=mc ){
      ib = min( m-i, mc );
      InnerKernel( ib, n, pb, &A( i,p ), lda, &B(p, 0 ), ldb, &C( i,0 ), ldc );
    }
  }

外层的两个循环改为对k（矩阵A的列、B的行）和m（矩阵A的行）的遍历，同时分别设定这两个遍历的自增值kc、mc，每次调用完一个InnerKernel，即得到一块大小为mc x n的C矩阵元素。但该过程并没有将这一块大小为mc x n的C矩阵元素的结果算完（因为对k遍历的for循环在最外层）。

对由于分块带来的变化的分析（kc和mc）：

把原本的AddDot4x4作为micro kernel放到InnerKernel中，InnerKernel作为kernel；
每次内层循环InnerKernel计算一块规模为mc x n大小的矩阵C，但并没有完全算完，只算了pb/k，其中pb = min( k-p, kc )。pb作为最内层新的k；
pb = min( k-p, kc )和ib = min( m-i, mc )是用来处理最后剩下的边界问题；
外层两个for循环遍历的是k和m，每次自增大小分别为kc和mc（这也是矩阵分块blocking的大小），进一步有pb = min( k-p, kc )和ib = min( m-i, mc )；
InnerKernel中两个for循环遍历的是m和n，每次自增均为4（前文中提到，这个值自己可以调）；
AddDot4x4中的循环，每次自增1，和上一次的代码没有变化，但需要注意传进来的k值不是矩阵A的列数（B的行数），而是pb，即一个AddDot4x4，对矩阵C中的元素只算了pb/k，其中pb = min( k-p, kc )，此外AddDot4x4本身被InnerKernel（行列自增值都为4）调用，每次算C矩阵具体到每个元素的完成度只有：4x4xpb，即A的4行的前pb列，B的4列的前pb行，其中pb = min( k-p, kc )。

4.2 对现有代码加入blocking

考虑到参考代码是列优先，我们是行优先，那么需要对矩阵C做blocking时，做以下变化，同时修改变量名，如inc_p，inc_j，更清晰：

/* Block sizes */
#define inc_j 256  // inc_j along n
#define inc_p 128  // inc_p along k
#define min( i, j ) ( (i)<(j) ? (i): (j) )

/* This time, we compute a m x block_n block of C by a call to the InnerKernel */
for ( p=0; p<k; p+=inc_p )
{
    block_k = min( k-p, inc_p );
    for ( j=0; j<n; j+=inc_j )
    {
        block_n = min( n-j, inc_j );
        InnerKernel( m, block_n, block_k, &A( i,p ), lda, &B(p, 0 ), ldb, &C( i,0 ), ldc );
    }
}

#define A(i, j) a[lda*(i) + (j)]
#define B(i, j) b[ldb*(i) + (j)]
#define C(i, j) c[ldc*(i) + (j)]

/* Block sizes */
#define inc_j 256  // inc_j along n
#define inc_p 128  // inc_p along k
#define min( i, j ) ( (i)<(j) ? (i): (j) )

/* This time, we compute a m x block_n block of C by a call to the InnerKernel */
for ( p=0; p<k; p+=inc_p )
{
    block_k = min( k-p, inc_p );
    for ( j=0; j<n; j+=inc_j )
    {
        block_n = min( n-j, inc_j );
        InnerKernel( m, block_n, block_k, &A( i,p ), lda, &B(p, 0 ), ldb, &C( i,0 ), ldc );
    }
}

void InnerKernel( const int m, const int n, const int k,
                  float *a, const int lda, 
                  float *b, const int ldb,
                  float *c, const int ldc )
{
    for( int i=0; i<m; i+=4 )
    {
        for( int j=0; j<n; j+=4 )
        {
            AddDot4x4( k, &A( i,0 ), lda, &B( 0,j ), ldb, &C( i,j ), ldc );
        }
    }
}

#include <mmintrin.h>
#include <xmmintrin.h>  // SSE
#include <pmmintrin.h>  // SSE2
#include <emmintrin.h>  // SSE3

typedef union v4f
{
  __m128 v;
  float s[4];
} v4f_t;

void AddDot4x4( const int k, 
                float *a, const int lda,
                float *b, const int ldb,
                float *c, const int ldc)
{
    v4f_t
      c_00_01_02_03_vreg,
      c_10_11_12_13_vreg,
      c_20_21_22_23_vreg,
      c_30_31_32_33_vreg,

      b_p0_p1_p2_p3_vreg,

      a_0p_x4_vreg,
      a_1p_x4_vreg,
      a_2p_x4_vreg,
      a_3p_x4_vreg;

    float 
      *b_p0_pntr = &B( 0, 0 );

    c_00_01_02_03_vreg.v = _mm_setzero_ps();
    c_10_11_12_13_vreg.v = _mm_setzero_ps();
    c_20_21_22_23_vreg.v = _mm_setzero_ps();
    c_30_31_32_33_vreg.v = _mm_setzero_ps();

    for(int p = 0; p < k; ++p)
    {
        // duplicate each element in A as float4
        a_0p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 0, p ) );
        a_1p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 1, p ) );
        a_2p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 2, p ) );
        a_3p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 3, p ) );

        b_p0_pntr += p * n;
        b_p0_p1_p2_p3_vreg.v = __mm_load_ps( (float *) b_p0_pntr );

        c_00_01_02_03_vreg.v += a_0p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
        c_10_11_12_13_vreg.v += a_1p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
        c_20_21_22_23_vreg.v += a_2p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
        c_30_31_32_33_vreg.v += a_3p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
    }

    // 1st row of C
    C( 0, 0 ) += c_00_01_02_03_vreg.s[0];
    C( 0, 1 ) += c_00_01_02_03_vreg.s[1];
    C( 0, 2 ) += c_00_01_02_03_vreg.s[2];
    C( 0, 3 ) += c_00_01_02_03_vreg.s[3];
    // 2nd row of C
    C( 1, 0 ) += c_10_11_12_13_vreg.s[0];
    C( 1, 1 ) += c_10_11_12_13_vreg.s[1];
    C( 1, 2 ) += c_10_11_12_13_vreg.s[2];
    C( 1, 3 ) += c_10_11_12_13_vreg.s[3];
    // 3rd row of C
    C( 2, 0 ) += c_20_11_12_13_vreg.s[0];
    C( 2, 1 ) += c_20_11_12_13_vreg.s[1];
    C( 2, 2 ) += c_20_11_12_13_vreg.s[2];
    C( 2, 3 ) += c_20_11_12_13_vreg.s[3];
    // 4th row of C
    C( 3, 0 ) += c_30_31_32_33_vreg.s[0];
    C( 3, 1 ) += c_30_31_32_33_vreg.s[1];
    C( 3, 2 ) += c_30_31_32_33_vreg.s[2];
    C( 3, 3 ) += c_30_31_32_33_vreg.s[3];
}

ysh329 commented 5 years ago

5. 4x4 提升对矩阵A打包（pack）

5.1 参考代码分析

参考：Optimization_4x4_12 · flame/how-to-optimize-gemm Wiki

packA的操作将性能提升2倍，但是该过程分两步，第一步会将性能下降10%。因为InnerKernel中每次两个外层循环计算AddDot4x4前，都会对A的 4 行 k 列做PackMatrixA（这里k = block_k = min( k-p, inc_p );），即从A中取出一个4 x k规模并转存到packedA变量中，但这个过程在两层for循环中，导致第一次取了后，后面重复计算了（下一节我们会讲如何跳过这个重复的步骤）。列优先的部分参考代码如下：

// 参考代码：packA，列优先
void InnerKernel( int m, int n, int k, double *a, int lda, 
                                       double *b, int ldb,
                                       double *c, int ldc )
{
  int i, j;
  double 
    packedA[ m * k ];

  for ( j=0; j<n; j+=4 ){        /* Loop over the columns of C, unrolled by 4 */
    for ( i=0; i<m; i+=4 ){        /* Loop over the rows of C */
      /* Update C( i,j ), C( i,j+1 ), C( i,j+2 ), and C( i,j+3 ) in
     one routine (four inner products) */
      PackMatrixA( k, &A( i, 0 ), lda, &packedA[ i*k ] );
      AddDot4x4( k, &packedA[ i*k ], 4, &B( 0,j ), ldb, &C( i,j ), ldc );
    }
  }
}

5.2 实现对矩阵B的部分打包（pack）

因参考代码是列优先，我们代码是行优先，因而是对矩阵B做打包（pack）操作。

// packB，行优先
void InnerKernel( int m, int n, int k, float*a, int lda, 
                                       float *b, int ldb,
                                       float*c, int ldc )
{
    int i, j;
    float packedB[ k * n ];

    for ( i=0; i<m; i+=4 ) /* Loop over the rows of C */
    {        
        for ( j=0; j<n; j+=4 ) /* Loop over the columns of C, unrolled by 4 */
        {        /* Update C( i,j ), C( i,j+1 ), C( i,j+2 ), and C( i,j+3 ) in
                     one routine (four inner products) */
            PackMatrixB( k, &B( 0, j ), ldb, &packedB[ k*j ] );
            AddDot4x4( k, &packedB[ k*j ], 4, &B( 0,j ), ldb, &C( i,j ), ldc );
        }
    }
}

以上代码操作过程是，PackMatrixB每次取B矩阵的4列k行做pack进packedB变量中，下面是包括了PackMatrixB实现的整体代码：

#define A(i, j) a[lda*(i) + (j)]
#define B(i, j) b[ldb*(i) + (j)]
#define C(i, j) c[ldc*(i) + (j)]

/* Block sizes */
#define inc_j 256  // inc_j along n
#define inc_p 128  // inc_p along k
#define min( i, j ) ( (i)<(j) ? (i): (j) )

/* This time, we compute a m x block_n block of C by a call to the InnerKernel */
for ( p=0; p<k; p+=inc_p )
{
    block_k = min( k-p, inc_p );
    for ( j=0; j<n; j+=inc_j )
    {
        block_n = min( n-j, inc_j );
        InnerKernel( m, block_n, block_k, &A( i,p ), lda, &B(p, 0 ), ldb, &C( i,0 ), ldc );
    }
}

void InnerKernel( const int m, const int n, const int k,
                  float *a, const int lda, 
                  float *b, const int ldb,
                  float *c, const int ldc )
{
    int i, j;
    float packedB[ k * n ];
    for ( i=0; i<m; i+=4 ) /* Loop over the rows of C, unrolled by 4 */
    {        
        for ( j=0; j<n; j+=4 ) /* Loop over the columns of C, unrolled by 4 */
        {
          /* Update C( i,j ), C( i,j+1 ), C( i,j+2 ), and C( i,j+3 ) in
         one routine (four inner products) */
            PackMatrixB( k, &B( 0, j ), ldb, &packedB[ k*j ] );
            AddDot4x4( k, &packedB[ k*j ], 4, &B( 0,j ), ldb, &C( i,j ), ldc );
        }
    }
}

// B是逐行取，和intrinsic对应，第p行的4个，
void PackMatrixB( const int k, float *b, const int ldb, float * b_to)
{
    for(int p = 0; p < k; ++p)
    {
        float *b_floatx4_row = &B( p, 0 );

        *b_to     = *b_floatx4_row;
        *(b_to+1) = *(b_floatx4_row+1);
        *(b_to+2) = *(b_floatx4_row+2);
        *(b_to+3) = *(b_floatx4_row+3);

        b_to += 4;
    }
}

#include <mmintrin.h>
#include <xmmintrin.h>  // SSE
#include <pmmintrin.h>  // SSE2
#include <emmintrin.h>  // SSE3

typedef union v4f
{
  __m128 v;
  float s[4];
} v4f_t;

void AddDot4x4( const int k, 
                float *a, const int lda,
                float *b, const int ldb,
                float *c, const int ldc)
{
    v4f_t
      c_00_01_02_03_vreg,
      c_10_11_12_13_vreg,
      c_20_21_22_23_vreg,
      c_30_31_32_33_vreg,

      b_p0_p1_p2_p3_vreg,

      a_0p_x4_vreg,
      a_1p_x4_vreg,
      a_2p_x4_vreg,
      a_3p_x4_vreg;

    float 
      *b_p0_pntr = &B( 0, 0 );

    c_00_01_02_03_vreg.v = _mm_setzero_ps();
    c_10_11_12_13_vreg.v = _mm_setzero_ps();
    c_20_21_22_23_vreg.v = _mm_setzero_ps();
    c_30_31_32_33_vreg.v = _mm_setzero_ps();

    for(int p = 0; p < k; ++p)
    {
        // duplicate each element in A as float4
        a_0p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 0, p ) );
        a_1p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 1, p ) );
        a_2p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 2, p ) );
        a_3p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 3, p ) );

        b_p0_pntr += p * n;
        b_p0_p1_p2_p3_vreg.v = __mm_load_ps( (float *) b_p0_pntr );

        c_00_01_02_03_vreg.v += a_0p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
        c_10_11_12_13_vreg.v += a_1p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
        c_20_21_22_23_vreg.v += a_2p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
        c_30_31_32_33_vreg.v += a_3p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
    }

    // 1st row of C
    C( 0, 0 ) += c_00_01_02_03_vreg.s[0];
    C( 0, 1 ) += c_00_01_02_03_vreg.s[1];
    C( 0, 2 ) += c_00_01_02_03_vreg.s[2];
    C( 0, 3 ) += c_00_01_02_03_vreg.s[3];
    // 2nd row of C
    C( 1, 0 ) += c_10_11_12_13_vreg.s[0];
    C( 1, 1 ) += c_10_11_12_13_vreg.s[1];
    C( 1, 2 ) += c_10_11_12_13_vreg.s[2];
    C( 1, 3 ) += c_10_11_12_13_vreg.s[3];
    // 3rd row of C
    C( 2, 0 ) += c_20_11_12_13_vreg.s[0];
    C( 2, 1 ) += c_20_11_12_13_vreg.s[1];
    C( 2, 2 ) += c_20_11_12_13_vreg.s[2];
    C( 2, 3 ) += c_20_11_12_13_vreg.s[3];
    // 4th row of C
    C( 3, 0 ) += c_30_31_32_33_vreg.s[0];
    C( 3, 1 ) += c_30_31_32_33_vreg.s[1];
    C( 3, 2 ) += c_30_31_32_33_vreg.s[2];
    C( 3, 3 ) += c_30_31_32_33_vreg.s[3];
}

5.3 B部分打包后的加速

这一步操作就可将性能提升2倍，InnerKernel中有两层for循环，这步操作因为外层循环会复用内层循环的packedB，因而不需要重复做，这步只是加了一个if条件判断：

void InnerKernel( const int m, const int n, const int k,
                  float *a, const int lda, 
                  float *b, const int ldb,
                  float *c, const int ldc )
{
    int i, j;
    float packedB[ k * n ];
    for ( i=0; i<m; i+=4 ) /* Loop over the rows of C, unrolled by 4 */
    {        
        for ( j=0; j<n; j+=4 ) /* Loop over the columns of C, unrolled by 4 */
        {
            /* Update C( i,j ), C( i,j+1 ), C( i,j+2 ), and C( i,j+3 ) in
           one routine (four inner products) */
            if( i == 0 ) PackMatrixB( k, &B( 0, j ), ldb, &packedB[ k*j ] );
            AddDot4x4( k, &packedB[ k*j ], 4, &B( 0,j ), ldb, &C( i,j ), ldc );
        }
    }
}

ysh329 commented 5 years ago

6. 4x4 提升对矩阵B打包（pack）

6.1 参考代码分析

参考：Optimization_4x4_14 · flame/how-to-optimize-gemm Wiki

打包（pack）矩阵B的kx4个元素，因为在for循环里重复打包所以会有些许性能下降，同样，再加入if判断，让其不会重复打包。下面是列优先、double类型的参考代码：

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, p, pb, ib;

  /* This time, we compute a mc x n block of C by a call to the InnerKernel */

  for ( p=0; p<k; p+=kc ){
    pb = min( k-p, kc );
    for ( i=0; i<m; i+=mc ){
      ib = min( m-i, mc );
      InnerKernel( ib, n, pb, &A( i,p ), lda, &B(p, 0 ), ldb, &C( i,0 ), ldc, i==0 );
    }
  }
}

void InnerKernel( int m, int n, int k, double *a, int lda, 
                                       double *b, int ldb,
                                       double *c, int ldc, int first_time )
{
  int i, j;
  double 
    packedA[ m * k ];
  static double 
    packedB[ kc*nb ];    /* Note: using a static buffer is not thread safe... */

  for ( j=0; j<n; j+=4 ){        /* Loop over the columns of C, unrolled by 4 */
    if ( first_time )
      PackMatrixB( k, &B( 0, j ), ldb, &packedB[ j*k ] );
    for ( i=0; i<m; i+=4 ){        /* Loop over the rows of C */
      /* Update C( i,j ), C( i,j+1 ), C( i,j+2 ), and C( i,j+3 ) in
     one routine (four inner products) */
      if ( j == 0 ) 
    PackMatrixA( k, &A( i, 0 ), lda, &packedA[ i*k ] );
      AddDot4x4( k, &packedA[ i*k ], 4, &packedB[ j*k ], k, &C( i,j ), ldc );
    }
  }
}

void PackMatrixA( int k, double *a, int lda, double *a_to )
{
  int j;

  for( j=0; j<k; j++){  /* loop over columns of A */
    double 
      *a_ij_pntr = &A( 0, j );

    *a_to     = *a_ij_pntr;
    *(a_to+1) = *(a_ij_pntr+1);
    *(a_to+2) = *(a_ij_pntr+2);
    *(a_to+3) = *(a_ij_pntr+3);

    a_to += 4;
  }
}

void PackMatrixB( int k, double *b, int ldb, double *b_to )
{
  int i;
  double 
    *b_i0_pntr = &B( 0, 0 ), *b_i1_pntr = &B( 0, 1 ),
    *b_i2_pntr = &B( 0, 2 ), *b_i3_pntr = &B( 0, 3 );

  for( i=0; i<k; i++){  /* loop over rows of B */
    *b_to++ = *b_i0_pntr++;
    *b_to++ = *b_i1_pntr++;
    *b_to++ = *b_i2_pntr++;
    *b_to++ = *b_i3_pntr++;
  }
}

6.2 实现对矩阵A的部分打包（pack）

现在，跟着上一步，考虑到我们的代码是行优先、float类型。

#define A(i, j) a[lda*(i) + (j)]
#define B(i, j) b[ldb*(i) + (j)]
#define C(i, j) c[ldc*(i) + (j)]

/* Block sizes */
#define inc_j 256  // inc_j along n
#define inc_p 128  // inc_p along k
#define min( i, j ) ( (i)<(j) ? (i): (j) )

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
    /* This time, we compute a m x block_n block of C by a call to the InnerKernel */
    for ( p=0; p<k; p+=inc_p )
    {
        block_k = min( k-p, inc_p );
        for ( j=0; j<n; j+=inc_j )
        {
            block_n = min( n-j, inc_j );
            InnerKernel( m, block_n, block_k, &A( i,p ), lda, &B(p, 0 ), ldb, &C( i,0 ), ldc, j==0 );
        }
    }
}

void InnerKernel( const int m, const int n, const int k,
                  float *a, const int lda, 
                  float *b, const int ldb,
                  float *c, const int ldc,
                  const int first_time)
{
    int i, j;
    float packedB[ k * n ];
    for ( i=0; i<m; i+=4 ) /* Loop over the rows of C, unrolled by 4 */
    {
        if (first_time) PackMatrixA( k, &A( i, 0 ), lda, &packedA[ i*k ] );
        for ( j=0; j<n; j+=4 ) /* Loop over the columns of C, unrolled by 4 */
        {
          /* Update C( i,j ), C( i,j+1 ), C( i,j+2 ), and C( i,j+3 ) in
         one routine (four inner products) */
            if( i == 0 ) PackMatrixB( k, &B( 0, j ), ldb, &packedB[ k*j ] );
            AddDot4x4( k, &packedB[ k*j ], 4, &B( 0,j ), ldb, &C( i,j ), ldc );
        }
    }
}

// B是逐行取，和intrinsic对应，第p行的4个，
void PackMatrixB( const int k, float *b, const int ldb, float * b_to)
{
    for(int p = 0; p < k; ++p)
    {
        float *b_floatx4_row = &B( p, 0 );

        *b_to     = *b_floatx4_row;
        *(b_to+1) = *(b_floatx4_row+1);
        *(b_to+2) = *(b_floatx4_row+2);
        *(b_to+3) = *(b_floatx4_row+3);

        b_to += 4;
    }
}

// A是重复取，和intrinsic对应，第p列的第一个重复4次
void PackMatrixA( const int k, float *a, const int lda, float * a_to)
{
    // 【】【】【】【】【】【】】【】【】
  double 
    *b_i0_pntr = &B( 0, 0 ), *b_i1_pntr = &B( 0, 1 ),
    *b_i2_pntr = &B( 0, 2 ), *b_i3_pntr = &B( 0, 3 );

  for( i=0; i<k; i++){  /* loop over rows of B */
    *b_to++ = *b_i0_pntr++;
    *b_to++ = *b_i1_pntr++;
    *b_to++ = *b_i2_pntr++;
    *b_to++ = *b_i3_pntr++;
  }
   // 【】【】【】【】【】【】【】

    float *a_
    for(int p = 0; p < k; ++p)
    {

    }
}

#include <mmintrin.h>
#include <xmmintrin.h>  // SSE
#include <pmmintrin.h>  // SSE2
#include <emmintrin.h>  // SSE3

typedef union v4f
{
  __m128 v;
  float s[4];
} v4f_t;

void AddDot4x4( const int k, 
                float *a, const int lda,
                float *b, const int ldb,
                float *c, const int ldc)
{
    v4f_t
      c_00_01_02_03_vreg,
      c_10_11_12_13_vreg,
      c_20_21_22_23_vreg,
      c_30_31_32_33_vreg,

      b_p0_p1_p2_p3_vreg,

      a_0p_x4_vreg,
      a_1p_x4_vreg,
      a_2p_x4_vreg,
      a_3p_x4_vreg;

    float 
      *b_p0_pntr = &B( 0, 0 );

    c_00_01_02_03_vreg.v = _mm_setzero_ps();
    c_10_11_12_13_vreg.v = _mm_setzero_ps();
    c_20_21_22_23_vreg.v = _mm_setzero_ps();
    c_30_31_32_33_vreg.v = _mm_setzero_ps();

    for(int p = 0; p < k; ++p)
    {
        // duplicate each element in A as float4
        a_0p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 0, p ) );
        a_1p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 1, p ) );
        a_2p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 2, p ) );
        a_3p_x4_vreg.v = __mm_loaddup_ps( (float *) &A( 3, p ) );

        b_p0_pntr += p * n;
        b_p0_p1_p2_p3_vreg.v = __mm_load_ps( (float *) b_p0_pntr );

        c_00_01_02_03_vreg.v += a_0p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
        c_10_11_12_13_vreg.v += a_1p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
        c_20_21_22_23_vreg.v += a_2p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
        c_30_31_32_33_vreg.v += a_3p_x4_vreg.v * b_p0_p1_p2_p3_vreg.v;
    }

    // 1st row of C
    C( 0, 0 ) += c_00_01_02_03_vreg.s[0];
    C( 0, 1 ) += c_00_01_02_03_vreg.s[1];
    C( 0, 2 ) += c_00_01_02_03_vreg.s[2];
    C( 0, 3 ) += c_00_01_02_03_vreg.s[3];
    // 2nd row of C
    C( 1, 0 ) += c_10_11_12_13_vreg.s[0];
    C( 1, 1 ) += c_10_11_12_13_vreg.s[1];
    C( 1, 2 ) += c_10_11_12_13_vreg.s[2];
    C( 1, 3 ) += c_10_11_12_13_vreg.s[3];
    // 3rd row of C
    C( 2, 0 ) += c_20_11_12_13_vreg.s[0];
    C( 2, 1 ) += c_20_11_12_13_vreg.s[1];
    C( 2, 2 ) += c_20_11_12_13_vreg.s[2];
    C( 2, 3 ) += c_20_11_12_13_vreg.s[3];
    // 4th row of C
    C( 3, 0 ) += c_30_31_32_33_vreg.s[0];
    C( 3, 1 ) += c_30_31_32_33_vreg.s[1];
    C( 3, 2 ) += c_30_31_32_33_vreg.s[2];
    C( 3, 3 ) += c_30_31_32_33_vreg.s[3];
}

yuenshome / yuenshome.github.io

gemm optimize #31

内存、cache和寄存器之间的关系及区别

计算机的存储体系（Memory hierarchy）金字塔

计算机的存储体系

为什么寄存器（register）比内存快？

原因一：距离不同

原因二：硬件设计不同

原因三：工作方式不同

gemm optimize

1. 1x1到1x4

1.1 naive 1x1

1.2 1x1

1.3 1x4

1.4 1x4

1.5 1x4的提升 L2 Cache

1.6 1x4的提升 C和A放到Register

1.7 1x4的提升 用指针访问B

1.8 1x4的轻微下降 对k方向展开p+=4，这个累加的数也可以自定

1.9 1x4 间接寻址（indirect addressing）

2. 1x4到4x4

2.1 4x4 A的行 B的列 每次取4

2.2 4x4 拆分AddDot到AddDot4x4中

2.3 4x4 提升 合并AddDot4x4中的For循环

2.4 4x4 提升 C和A放到Register中

2.5 4x4 提升 用指针代替循环中的B

3. 进一步优化4x4（Further optimizing）

3.1 4x4 下降 将循环中的B用指针索引后，放到寄存器

3.2 4x4 提升 变为vector operation进一步转为intrinsic

英特尔® 高级矢量扩展指令集简介

Optimization Notice

指令集概述

内存对齐

4. 4x4 提升 利用L2 Cache对矩阵blocking

4.1 参考代码分析

4.2 对现有代码加入blocking

5. 4x4 提升 对矩阵A打包（pack）

5.1 参考代码分析

5.2 实现对矩阵B的部分打包（pack）

5.3 B部分打包后的加速

6. 4x4 提升 对矩阵B打包（pack）

6.1 参考代码分析

6.2 实现对矩阵A的部分打包（pack）

6.3 A部分打包后的加速

1.7 1x4的提升用指针访问B

1.8 1x4的轻微下降对k方向展开p+=4，这个累加的数也可以自定

2.1 4x4 A的行 B的列每次取4

2.3 4x4 提升合并AddDot4x4中的For循环

2.5 4x4 提升用指针代替循环中的B

3.1 4x4 下降将循环中的B用指针索引后，放到寄存器

3.2 4x4 提升变为vector operation进一步转为intrinsic

4. 4x4 提升利用L2 Cache对矩阵blocking

5. 4x4 提升对矩阵A打包（pack）

6. 4x4 提升对矩阵B打包（pack）