高性能计算例程：GEMM
一、背景：
- GEMM(General Matrix Multiplication）全称为“通用矩阵乘法”；
- GEMM是BLAS(Basic Linear Algebra Subprograms)中Level 3的例程，也是整个BLAS最重要的例程之一；
二、例程(以double类型为例)：
1.功能：
计算
α
?
A
′
?
B
′
+
β
?
C
\alpha*A'*B'+\beta*C
α?A′?B′+β?C的值，然后将结果存入C中,
即：
C
?
α
?
A
′
?
B
′
+
β
?
C
C\longleftarrow \alpha*A'*B'+\beta*C
C?α?A′?B′+β?C
其中：
（1）
α
和
β
为标量常数
\alpha \text{和} \beta \text{为标量常数}
α和β为标量常数
（2）
A
′
=
o
p
(
A
)
=
{
A
,
op?is?no?transpose
A
T
,
op?is?transpose
A'=op(A)=\left\{ \begin{aligned} &A & , & \text{op is no transpose} \\ &A^T& , & \text{op is transpose} \\ \end{aligned} \right.
A′=op(A)={?AAT?,,?op?is?no?transposeop?is?transpose?
（3）
B
′
=
o
p
(
B
)
=
{
B
,
op?is?no?transpose
B
T
,
op?is?transpose
B'=op(B)=\left\{ \begin{aligned} &B & , & \text{op is no transpose} \\ &B^T& , & \text{op is transpose} \\ \end{aligned} \right.
B′=op(B)={?BBT?,,?op?is?no?transposeop?is?transpose?
(4)
A
′
?is?
k
?by?
m
?:?
(
A
′
)
m
×
k
A' \text{ is }k\text{ by } m \text{ : }\qquad (A')_{m \times k}
A′?is?k?by?m?:?(A′)m×k?
(5)
B
′
?is?
n
?by?
k
?:?
(
B
′
)
k
×
n
B' \text{ is }n\text{ by } k \text{ : }\qquad (B')_{k \times n}
B′?is?n?by?k?:?(B′)k×n?
(6)
C
?is?
n
?by?
m
?:?
(
C
)
m
×
n
C \text{ is }n\text{ by } m \text{ : }\qquad (C)_{m \times n}
C?is?n?by?m?:?(C)m×n?
即：
C
=
α
?
o
p
(
A
)
?
o
p
(
B
)
+
β
?
C
C =\alpha*op(A)*op(B)+\beta*C
C=α?op(A)?op(B)+β?C
2.例程：
void cblas_dgemm
(
const CBLAS_LAYOUT Layout,
const CBLAS_TRANSPOSE transa,
const CBLAS_TRANSPOSE transb,
const CBLAS_INT m,
const CBLAS_INT n,
const CBLAS_INT k,
const double alpha,
const double *a,
const CBLAS_INT lda,
const double *b,
const CBLAS_INT ldb,
const double beta,
double *c,
const CBLAS_INT ldc
);
参数说明：
-
(1) const CBLAS_LAYOUT Layout:
矩阵的存储方式：列优先或者行优先。
列优先：Column-Major
行优先：Row-Major
作者（我 E2MCC）补充：
其实矩阵的有多种存储方式，如块行优先存储，块列优先存储，对于稀疏矩阵而言还有 COO(Triplet) CSC CSR ELL DIA 等存储方式。这里的 ColMaj 和RowMaj 是最常见的稠密矩阵的基本存储方式。 -
(2) const CBLAS_TRANSPOSE transa:
对
A
A
A 进行转置操作：转置或不转置（Transpose
A
A
A or not） -
(3) const CBLAS_TRANSPOSE transb:
对
B
B
B 进行转置操作：转置或不转置（Transpose
B
B
B or not） -
(4) const CBLAS_INT m:
m 表示实际运算时
A
A
A 矩阵的行维度(行数)，也就是 op
(
A
)
(A)
(A) 的行维度；
当然也是结果矩阵
C
C
C 的行维度(行数)。 -
(5) const CBLAS_INT k:
k 表示实际运算时
A
A
A 矩阵的列维度(列数)，也就是 op
(
A
)
(A)
(A) 的列维度；
同时也是实际运算时
B
B
B 矩阵的行维度(行数)，也就是 op
(
B
)
(B)
(B) 的行维度。 -
(6) const CBLAS_INT n:
n 表示实际运算时
B
B
B 矩阵的列维度(列数)，也就是 op
(
B
)
(B)
(B) 的列维度； -
(7) const double alpha:
常数项
α
\alpha
α。 -
(8) const double * a :
a 为一个指针，指向
A
A
A 矩阵数据在内存中的位置。 -
(9) const CBLAS_INT lda:
Leading Dimension of A：表示
A
A
A 矩阵数据在内存中存储时的步长。
注意：这里 lda 是指的 A 的 leading dimension 而不是 op
(
A
)
(A)
(A)的 leading dimension。 -
(10) const double * b:
b为一个指针，指向
B
B
B 矩阵数据在内存中的位置。 -
(11) const CBLAS_INT ldb:
Leading Dimension of B：表示
B
B
B 矩阵数据在内存中存储时的步长。
注意：这里 ldb 是指的 B 的 leading dimension 而不是 op
(
B
)
(B)
(B)的 leading dimension。 -
(12) const double beta:
常数项
β
\beta
β。 -
(13) double * c:
c 为一个指针，指向
C
C
C 矩阵数据在内存中的位置。 -
(14) const CBLAS_INT ldc:
Leading Dimension of C：表示
C
C
C 矩阵数据在内存中存储时的步长。
注意：这里 ldc 是指的 C 的 leading dimension 而不是 op
(
C
)
(C)
(C)的 leading dimension。
3.Leading dimension详细解释
（1）leading dimension 是什么?
leading dimension 是二维逻辑矩阵映射到一维内存空间时的映射方程的系数。
（2）为什么要使用leading dimension？leading dimension和m，n，k有什么区别？
当A矩阵就是一整个矩阵时，leading dimension确实为A矩阵的行数或者列数，但是若A矩阵仅仅只是一个更大的矩阵中的一部分时，那么leading dimension 必定不等于A的行数或列数了。如下图所示：
|