高性能计算-gemm串行计算优化(3)

目标：Darknet 源码cpu矩阵乘法函数 gemm_nn 优化。参数说明：lda A的列数; ldb B的列数; ldc C的列数; M C的行数; K A的列数

测试方法：Darknet源码，makefile文件添加编译选项 -pg，编译后得到可执行程序 darknet，运行可执行程序：

./darknet detect cfg/yolov3-tiny.cfg yolov3-tiny.weights data/kite.jpg

得到文件 gmon.out，转换为文本，查看热点函数及其耗时。

gprof darknet gmon.out > out 

1. 矩阵乘法基础写法

for(int i=0; i<M; i++)//矩阵A: M*K
{
  for(int j=o; j<N; j++)//矩阵B: K*N
  {
		for(int k=0; k<K; k++)
		{
			C[i][j] += A[i][k] * B[k][j];
		}
  }
}

2. darknet 矩阵乘法源码 src/gemm.c

void gemm_nn(int M, int N, int K, float ALPHA, float *A, int lda, float*B, int ldb,loat *C, int ldc)
{
  int i,j,k;
	for (i = 0; i < M; ++i)//矩阵A: M*K
	{
		 for (j = 0; j < N; ++j)//矩阵B: K*N
		{
			register float c_sum = 0.0f;//使用局部变量减少对传入内存引用
			for (k = 0; k < K; ++k)//K 矩阵A列数；按行列顺序，逐一完成C 元素的计算
			{
				//lda 矩阵A的列数；lAb 矩阵B的列数
				c_sum += ALPHA * A[i * lda  + k] * B[k * ldb + j ];
			}
			//ldc 矩阵C的列数
			C[i * ldc +j] += c_sum;
		}
	}
}

3. 循环展开优化：2x和4x加速不明显

void gemm_nn(int M, int N, int K, float ALPHA, 
        float *A, int lda, 
        float *B, int ldb,
        float *C, int ldc)
{
    int i,j,k;
    for (i = 0; i < M; ++i)
    {
        for (j = 0; j < N; ++j)
        {
            register float c_sum = 0.0f;
			#if 0
			for (k = 0; k < K; ++k)//一个完整元素的完整计算过程 
            {
                c_sum += ALPHA * A[i * lda + k] * B[k * ldb + j];
            }
			#elif 1 //2x
			for (k = 0; k < K; k += 2)//一个完整元素的完整计算过程 
            {
                c_sum += ALPHA * A[i * lda + k] * B[k * ldb + j];
                c_sum += ALPHA * A[i * lda + k+1] * B[(k+1) * ldb + j];
            }
			#elif //4x
            for (k = 0; k < K; k += 4)//一个完整元素的完整计算过程 
            {
                c_sum += ALPHA * A[i * lda + k] * B[k * ldb + j];
                c_sum += ALPHA * A[i * lda + k+1] * B[(k+1) * ldb + j];
				c_sum += ALPHA * A[i * lda + k+2] * B[(k+2) * ldb + j];
                c_sum += ALPHA * A[i * lda + k+3] * B[(k+3) * ldb + j];
            }
			#elif 
			#endif
            C[i * ldc + j] += c_sum;
        }
    }
}

3. 矩阵乘法一维分块优化

void gemm_nn(int M, int N, int K, float ALPHA,
        float *A, int lda,
        float *B, int ldb,
        float *C, int ldc)
{
	memset(C,0,sizeof(float)*ldc*M);
	for(int m=0; m<M; m += 4)//矩阵A M*K 行方向遍历，定位块行号
	{
		for(int n=0; n<N; n += 4)//矩阵B K*N 列方向遍历，定位块列号
		{
			for(int k=0; k<K; k++)//K方向遍历，放在最内层也可
			{
				for(int mb=0; mb<4; mb++)//定位块内行号
				{
					for(int nb=0; nb<4; nb++)//定位快内列号
					{
						C[(m+mb)*ldc + (n+nb)] +=  ALPHA* A[(m+mb)*lda + (k)]*B[(k)*ldb + (n+nb)];
					}
				}
			}
		}
	}
}

4. gemm_nn 矩阵乘法二维分块优化：cannon 算法

void gemm_nn(int M, int N, int K, float ALPHA,
        float *A, int lda,
        float *B, int ldb,
        float *C, int ldc)
{
	memset(C,0,sizeof(float)*ldc*M);
	for(int m=0; m<M; m += 4)//矩阵A M*K M方向遍历，定位块行号
	{
		for(int n=0; n<N; n += 4)//矩阵B K*N N方向遍历，定位块列号
		{
			for(int k=0; k<K; k += 4)//K方向遍历，定位块号
			{
				for(int mb=0; mb<4; mb++)//定位块内行号
				{
					for(int nb=0; nb<4; nb++)//定位块内列号
					{
						for(int kb=0; kb<4; kb++)//快内按照 K 方向累加
						{
							C[(m+mb)*ldc + (n+nb)] = ALPHA * A [(m+mb)*lda + (k+kb)] * B[(k+kb)*ldb + (n+nb)];
						}
					}
				}
			}
		}
	}
}

5. 对B转置计算
   通常矩阵元素的计算公式：
   $$C[m][n]=\sum _0^{k-1}A[m][k]\ast B[k][n]$$
   优化后的矩阵元素计算公式：
   $$C[m][n]=\sum _0^{k-1}A[m][k]\ast B^T[n][k]$$
   代码实现：

void gemm_nn(int M, int N, int K, float ALPHA, float *A, int lda, float*B, int ldb,loat *C, int ldc)
{
	//对B转置再计算，可以增加缓存命中
        memset(C,0,sizeof(float)*ldc*M);
        float* Btemp=calloc(K*ldb,sizeof(float));
        for(int i=0;i<K;i++)
        {
          for(int j=0;j<ldb;j++)
            Btemp[j*K+i] = B[i*ldb+j];
        }
        for(int m=0;m<M;m++)
        {
          for(int n=0;n<ldb;n++)
          {
            for(int k=0;k<K;k++)
              C[m*ldc+n] += ALPHA * A[m*lda+k] * Btemp[n*K+k];
          }
        }
        free(Btemp);
}

6. 测试结果

gemm_nn 函数耗时：
darknet 源码	11.7s
loop 2x循环展开：11.5s
loop 4x循环展开：12.1s
一维分块4*K：9.2，比源码提升 27%
二维分块block 4*4 ：6.53s，比源码效率提升 79%
对B转置：5.17s，比源码提升 126%