高性能计算例程: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的行数或列数了。如下图所示: