General Matrix Multiply

From Wikipedia, the free encyclopedia

The General Matrix Multiply (GEMM) is a subroutine in the Basic Linear Algebra Subprograms (BLAS) which performs matrix multiplication, that is the multiplication of two matrices. This includes:

GEMM is often tuned by High Performance Computing vendors to run as fast as possible, because it is the building block for so many other routines. It is also the most important routine in the LINPACK benchmark. For this reason, implementations of fast BLAS library typically focus first on GEMM performance.

[edit] Operation

The xGEMM routine calculates the new value of matrix C based on the matrix-product of matrices A and B, and the old value of matrix C

C \leftarrow \alpha A B + \beta C

where α and β values are scalar coefficients.

The Fortran interface for these procedures are:

SUBROUTINE xGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC )

where TRANSA and TRANSB determines if the matrices A and B are to be transposed. M is the number of rows in matrix A and C. N is the number of columns in matrix B and C. K is the number of columns in matrix A and rows in matrix B. LDA, LDB and LDC specifies the size of the first dimension of the matrices, as laid out in memory; meaning the memory distance between the start of each row/column, depending on the memory structure.

The C interface for these procedures are similar:

void cblas_xgemm (const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,
                  const SCALAR alpha, const TYPE * A, const int lda, const TYPE * B, const int ldb, const SCALAR beta, TYPE * C, const int ldc)

where 'x' in name is s,d,c, or z meaning single precision real, double precision real, single precision complex, or double precision complex, and other types have similar meanings as in Fortran interface.

For C interface, higher API also exists. For example, in GNU Scientific Library, a corresponding interface is (take double precision real case for instance):

int gsl_blas_dgemm(CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB, double alpha, const gsl_matrix * A, const gsl_matrix * B, double beta, gsl_matrix * C)

where TransA is CblasNoTrans or CblasTrans for A or transposed A respectively, and TransB has similar choices.

[edit] Optimization

Not only is GEMM an important building block of other numeric software, it is often an important building block for calls to GEMM for larger matrices. By decomposing one or both of the input matrices into block matrices, GEMM can be used repeatedly on the smaller blocks to build up a result for the full matrix. This is one of the motivations for including the β parameter, so the results of previous blocks can be accumulated. Note that this requires the special case β = 1 which many implementations optimize for, thereby eliminating one multiplication for each value of C.

This decomposition allows for better locality of reference both in space and time of the data used in the product. This, in turn, takes advantage of the cache on the system. For systems with more than one level of cache, the blocking can be applied a second time to the order in which the blocks are used in the computation. Both of these levels of optimization are used in implementations such as ATLAS.