Developer Guide

Contents

Offloading oneMKL Computations onto the GPU

The Intel® oneAPI Math Kernel Library (oneMKL) improves performance with math routines for software applications that solve large computational problems. oneMKL provides BLAS, Sparse BLAS, and LAPACK linear algebra routines, fast Fourier transforms, vectorized math functions, random number generation functions, and other functionality.
The oneMKL distribution includes an examples directory which contains examples of various calls to oneMKL routines.
For more information about the Intel oneAPI Math Kernel Library, see:

OpenMP Directives to Offload oneMKL Computations

You can use OpenMP directives to offload oneMKL computations onto the GPU. There are two ways to do this.
One way involves using the Intel-specific OpenMP extension
target variant dispatch
directive. You would place the call to the oneMKL routine inside a
target variant dispatch
construct, as shown in the example below. In this example, arrays
A
,
B
, and
C
used in the multiplication are mapped to the device before the call to the oneMKL routine
cblas_dgemm
. The
use_device_ptr(A,B,C)
clause is used on the
target variant dispatch
directive to indicate that
A
,
B
, and
C
point to objects that have corresponding storage on the device. When
cblas_dgemm
is called, the corresponding device pointers for
A
,
B
, and
C
will be passed as arguments, and the device copies of
A
,
B
, and
C
will be used in the computation.
//============================================================== // Copyright © 2022 Intel Corporation // // SPDX-License-Identifier: MIT // ============================================================= // clang-format off #include <stdio.h> #include <stdlib.h> #include <math.h> #include <omp.h> #include "mkl.h" #include "mkl_omp_offload.h" #define min(x,y) (((x) < (y)) ? (x) : (y)) #define EPSILON 0.0001 int main() { double *A, *B, *C, *C_fl; int64_t m, n, k; double alpha, beta; double sum; int64_t i, j, q; int fail; printf ("\n This example computes real matrix C=alpha*A*B+beta*C using \n" " Intel oneMKL function dgemm, where A, B, and C are matrices and \n" " alpha and beta are double precision scalars\n\n"); m = 2000, k = 200, n = 1000; printf (" Initializing data for matrix multiplication C=A*B for matrix \n" " A(%li x %li) and matrix B(%li x %li)\n\n", m, k, k, n); alpha = 1.0; beta = 0.0; printf (" Allocating memory for matrices aligned on 64-byte boundary for better \n" " performance \n\n"); A = (double *)mkl_malloc( m * k * sizeof( double ), 64 ); B = (double *)mkl_malloc( k * n * sizeof( double ), 64 ); C = (double *)mkl_malloc( m * n * sizeof( double ), 64 ); C_fl = (double *)mkl_malloc( m*n*sizeof( double ), 64 ); if (A == NULL || B == NULL || C == NULL || C_fl == NULL) { printf( "\n ERROR: Cannot allocate memory for matrices. Exiting... \n\n"); return 1; } printf (" Intializing matrices \n\n"); for (i = 0; i < (m*k); i++) { A[i] = (double)(i+1); } for (i = 0; i < (k*n); i++) { B[i] = (double)(-i-1); } for (i = 0; i < (m*n); i++) { C[i] = 0.0; C_fl[i] = 0.0; } printf (" Computing matrix product using Intel oneMKL dgemm function via CBLAS interface \n\n"); #pragma omp target data map(to: A[0:m*k], B[0:k*n]) map(tofrom: C[0:m*n]) { #pragma omp target variant dispatch use_device_ptr(A, B, C) { cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, k, B, n, beta, C, n); } } printf ("\n Top left corner of matrix C: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C[j+i*n]); } printf ("\n"); } printf (" Computing matrix product using for-loops \n"); for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { sum = 0.0; for (q = 0; q < k; q++) { sum += A[k*i+q] * B[n*q+j]; } C_fl[n*i+j] = alpha * sum + beta * C_fl[n*i+j]; } } printf ("\n Top left corner of matrix C_fl: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C_fl[j+i*n]); } printf ("\n"); } printf (" Computations completed. Verifying... \n\n"); fail = 0; for (i = 0; i < (m*n); i++) { if (fabs(C[i] - C_fl[i]) > EPSILON) { fail = 1; break; } } if (fail) printf ("\n **** FAIL **** \n"); else printf ("\n **** PASS **** \n"); printf ("\n Deallocating memory \n\n"); mkl_free(A); mkl_free(B); mkl_free(C); return fail; }
Another way to inform the compiler that oneMKL computations should be offloaded onto the GPU is by using the OpenMP 5.1
dispatch
directive, as shown in the example below. In this example too, arrays
A
,
B
, and
C
are mapped to the device before the call to the oneMKL routine
cblas_dgemm
. When
cblas_dgemm
is called, the corresponding device pointers for
A
,
B
, and
C
will be passed as arguments, so the device copies of
A
,
B
, and
C
will be used in the computation.
The
use_device_ptr
clause is not needed on the
dispatch
directive. With OpenMP 5.1, the list of device pointers needed by the oneMKL routines is given in the oneMKL OpenMP offload header file,
mkl_omp_offload.h
, where the GPU variant function is declared. The user should carefully review the list of device pointers required in the oneMKL header file and make sure that the corresponding arrays are accessible from the device before calling the oneMKL routine.
Note that, depending on the version of the compiler you are using, you may need to add the compiler option
-fopenmp-version=51
in order for the
dispatch
directive to be accepted.
//============================================================== // Copyright © 2022 Intel Corporation // // SPDX-License-Identifier: MIT // ============================================================= // clang-format off #include <stdio.h> #include <stdlib.h> #include <math.h> #include <omp.h> #include "mkl.h" #include "mkl_omp_offload.h" #define min(x,y) (((x) < (y)) ? (x) : (y)) #define EPSILON 0.0001 int main() { double *A, *B, *C, *C_fl; int64_t m, n, k; double alpha, beta; double sum; int64_t i, j, q; int fail; printf ("\n This example computes real matrix C=alpha*A*B+beta*C using \n" " Intel oneMKL function dgemm, where A, B, and C are matrices and \n" " alpha and beta are double precision scalars\n\n"); m = 2000, k = 200, n = 1000; printf (" Initializing data for matrix multiplication C=A*B for matrix \n" " A(%li x %li) and matrix B(%li x %li)\n\n", m, k, k, n); alpha = 1.0; beta = 0.0; printf (" Allocating memory for matrices aligned on 64-byte boundary for better \n" " performance \n\n"); A = (double *)mkl_malloc( m * k * sizeof( double ), 64 ); B = (double *)mkl_malloc( k * n * sizeof( double ), 64 ); C = (double *)mkl_malloc( m * n * sizeof( double ), 64 ); C_fl = (double *)mkl_malloc( m*n*sizeof( double ), 64 ); if (A == NULL || B == NULL || C == NULL || C_fl == NULL) { printf( "\n ERROR: Cannot allocate memory for matrices. Exiting... \n\n"); return 1; } printf (" Intializing matrices \n\n"); for (i = 0; i < (m*k); i++) { A[i] = (double)(i+1); } for (i = 0; i < (k*n); i++) { B[i] = (double)(-i-1); } for (i = 0; i < (m*n); i++) { C[i] = 0.0; C_fl[i] = 0.0; } printf (" Computing matrix product using Intel oneMKL dgemm function via CBLAS interface \n\n"); #pragma omp target data map(to: A[0:m*k], B[0:k*n]) map(tofrom: C[0:m*n]) { #pragma omp target variant dispatch use_device_ptr(A, B, C) cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, k, B, n, beta, C, n); } printf ("\n Top left corner of matrix C: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C[j+i*n]); } printf ("\n"); } printf (" Computing matrix product using for-loops \n"); for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { sum = 0.0; for (q = 0; q < k; q++) { sum += A[k*i+q] * B[n*q+j]; } C_fl[n*i+j] = alpha * sum + beta * C_fl[n*i+j]; } } printf ("\n Top left corner of matrix C_fl: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C_fl[j+i*n]); } printf ("\n"); } printf (" Computations completed. Verifying... \n\n"); fail = 0; for (i = 0; i < (m*n); i++) { if (fabs(C[i] - C_fl[i]) > EPSILON) { fail = 1; break; } } if (fail) printf ("\n **** FAIL **** \n"); else printf ("\n **** PASS **** \n"); printf ("\n Deallocating memory \n\n"); mkl_free(A); mkl_free(B); mkl_free(C); return fail; }
Notes:
  • oneMKL routines expect the arrays/matrices to be on the device before the computation will be run. So the user has to map the data to the device, or allocate the data directly on the device, before calling a oneMKL routine.
  • If a oneMKL routine is not called from a
    target variant dispatch
    (or
    dispatch
    ) region, or if offload is disabled, then the oneMKL computations will be executed on the CPU.
  • Only one call to a oneMKL routine can be issued from an OpenMP
    target variant dispatch
    (or
    dispatch
    ) construct. If there are two consecutive calls to oneMKL routines, then the calls should be placed in separate
    target variant dispatch
    (or
    dispatch
    ) constructs.
Fortran
When calling oneMKL routines from Fortran code, be sure to add the following include statement:
include "mkl_omp_offload.f90
Also, if calling oneMKL Fortran API with
32-bit integers
, add the following module
use
statement:
use onemkl_blas_omp_offload_lp64
On the other hand, if calling oneMKL Fortran API with
64-bit integers
, add the following module
use
statement:
use onemkl_blas_omp_offload_ilp64
The following Fortran example illustrates how DGEMM is called from a Fortran program, and the
include
and
use
statements mentioned above.
!============================================================= ! Copyright © 2022 Intel Corporation ! ! SPDX-License-Identifier: MIT !============================================================= include "mkl_omp_offload.f90" program DGEMM_MAIN #if defined(MKL_ILP64) use onemkl_blas_omp_offload_ilp64 #else use onemkl_blas_omp_offload_lp64 #endif use omp_lib use iso_fortran_env implicit none integer, parameter :: m = 20 integer, parameter :: k = 5 integer, parameter :: n = 10 double precision a(m,k), b(k,n), c1(m,n), c2(m,n) double precision alpha, beta integer i, j print* print*,' D G E M M EXAMPLE PROGRAM' ! Initialize alpha = 1.025 beta = 0.75 do i = 1, m do j = 1, k a(i,j) = (i-1) - (0.25 * k) end do end do do i = 1, k do j = 1, n b(i,j) = -((i-1) + j) end do end do do i = 1, m do j = 1, n c1(i,j) = 0.2 + i - j c2(i,j) = 0.2 + i - j end do end do ! Execute DGEMM on host. call DGEMM('N','N',m,n,k,alpha,a,m,b,k,beta,c1,m) print * print *, 'c1 - After DGEMM host execution' do i=1,m print 110, (c1(i,j),j=1,n) end do print* ! Execute DGEMM on device !$omp target data map(to: a, b) map(tofrom: c2) !$omp dispatch call DGEMM('N','N',m,n,k,alpha,a,m,b,k,beta,c2,m) !$omp end target data print * print *, 'c2 - After DGEMM device execution' do i=1,m print 110, (c2(i,j),j=1,n) end do print * 101 format(7x,'M=',i5,' N=',i5,' K=',i5) 102 format(7x,'ALPHA=',f10.2,' BETA=',f10.2) 110 format(7x,10(f10.2,2x)) end
To compile and link the above Fortran example with
32-bit integers
:
ifx -fiopenmp -fopenmp-targets=spir64 -qmkl=parallel -fpp -free -c dgemm_example_f.f ifx -fiopenmp -fopenmp-targets=spir64 -qmkl=parallel -fsycl -L${MKLROOT}/lib/intel64 -liomp5 -lsycl -lOpenCL -lstdc++ -lpthread -lm -ldl -lmkl_sycl dgemm_example_f.o
To compile and link the above Fortran example with
64-bit integers
:
ifx -fiopenmp -fopenmp-targets=spir64 -qmkl=parallel -m64 -DMKL_ILP64 -i8 -fpp -free -c dgemm_example_f.f ifx -fiopenmp -fopenmp-targets=spir64 -qmkl=parallel -fsycl -L${MKLROOT}/lib/intel64 -liomp5 -lsycl -lOpenCL -lstdc++ -lpthread -lm -ldl -lmkl_sycl dgemm_example_f.o
After generating the executable (a.out), from a C/C++ or Fortran program, you can run the executable under
ze_tracer
and look for the heading “Device Timing Results” in the generated trace. Below that heading we should see the oneMKL kernels listed. This way we confirm that oneMKL computations have been offloaded onto the GPU.
Example run command:
OMP_TARGET_OFFLOAD=MANDATORY ZE_AFFINITY_MASK=0.0 ze_tracer -h -d ./a.out

Batching of oneMKL GEMM Calls

The oneMKL library includes “batch” routines that allow the user to batch several oneMKL calls into a single oneMKL call. At runtime, oneMKL will intelligently execute all of the matrix operations to optimize overall performance.
For example, the
cblas_dgemm
routine computes a matrix-matrix product of two general matrices
a
and
b
, returning the result in a matrix
c
. The
cblas_dgemm
interface is shown below.
void cblas_dgemm (const CBLAS_LAYOUT layout, const CBLAS_TRANSPOSE transa, const CBLAS_TRANSPOSE transb, const MKL_INT m, const MKL_INT n, const MKL_INT k, const double alpha, const double *a, const MKL_INT lda, const double *b, const MKL_INT ldb, const double beta, double *c, const MKL_INT ldc);
The
cblas_dgemm_batch
routine is similar to the
cblas_dgemm
routine, but the
cblas_dgemm_batch
routine performs matrix-matrix operations on groups of matrices, processing a number of groups at once.
The
cblas_dgemm_batch
interface is shown below. Note that the interface resembles the
cblas_dgemm
interface. However, it involves passing matrix arguments as arrays of pointers to matrices, and passing parameters as arrays of parameters.
void cblas_dgemm_batch (const CBLAS_LAYOUT layout, const CBLAS_TRANSPOSE* transa_array, const CBLAS_TRANSPOSE* transb_array, const MKL_INT* m_array, const MKL_INT* n_array, const MKL_INT* k_array, const double* alpha_array, const double **a_array, const MKL_INT* lda_array, const double **b_array, const MKL_INT* ldb_array, const double* beta_array, double **c_array, const MKL_INT* ldc_array, const MKL_INT group_count, const MKL_INT* group_size);
The batch operation is defined as follows:
idx = 0 for i = 0 .. group_count - 1 alpha and beta in alpha_array[i] and beta_array[i] for j = 0 .. group_size[i] - 1 a, b, and c matrices in a_array[idx], b_array[idx], and c_array[idx], respectively c := alpha*op(a)*op(b) + beta*c, idx = idx + 1 end for end for
where:
  • op(X) is one of op(X) = X, or op(X) = XT, or op(X) = XH,
  • alpha and beta are scalar elements of alpha_array and beta_array,
  • a, b and c are matrices such that for m, n, and k which are elements of m_array, n_array, and k_array:
  • op(a) is an m-by-k matrix,
  • op(b) is a k-by-n matrix,
  • C is an m-by-n matrix.
  • a, b, and c represent matrices stored at addresses pointed to by a_array, b_array, and c_array, respectively. The number of entries in a_array, b_array, and c_array is total_batch_count = the sum of all of the group_size entries.
It is possible to batch the multiplications of different shapes and parameters by packaging them into groups, where each group consists of multiplications of matrices of the same shapes (same m, n, and k) and the same parameters.
The basic assumption for the batch API are that all operations in a batch (whether in the same group or different groups) are independent of one another. So oneMKL does not guarantee any particular ordering between operations in a batch, and will try to execute multiple operations in parallel.
In general, the larger you can make the batch size, the better. This allows oneMKL to better parallelize the operations and distribute the work across the GPU.
We illustrate how two calls to
cblas_dgemm
can be replaced with one call to
cblas_dgemm_batch
. The following example includes two calls to
cblas_dgemm
.
//============================================================== // Copyright © 2022 Intel Corporation // // SPDX-License-Identifier: MIT // ============================================================= // clang-format off #include <stdio.h> #include <stdlib.h> #include <math.h> #include <omp.h> #include "mkl.h" #include "mkl_omp_offload.h" #define min(x,y) (((x) < (y)) ? (x) : (y)) #define epsilon 0.0000001f bool compare(double x, double y) { // returns true if x and y are the same return fabs(x - y) <= epsilon; } int main() { double *A1, *B1, *C1, *C1_fl; double *A2, *B2, *C2, *C2_fl; int m, n, k, i, j, q; double alpha, beta; double sum; int fail; double t_start, t_end; m = 2000, k = 200, n = 1000; alpha = 1.0; beta = 0.0; printf (" Allocating memory for matrices aligned on 64-byte boundary for better \n" " performance \n\n"); A1 = (double *)mkl_malloc (m*k*sizeof( double ), 64 ); B1 = (double *)mkl_malloc (k*n*sizeof( double ), 64 ); C1 = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); C1_fl = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); A2 = (double *)mkl_malloc (m*k*sizeof( double ), 64 ); B2 = (double *)mkl_malloc (k*n*sizeof( double ), 64 ); C2 = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); C2_fl = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); if (A1 == NULL || B1 == NULL || C1 == NULL || C1_fl == NULL || A2 == NULL || B2 == NULL || C2 == NULL || C2_fl == NULL) { printf( "\n ERROR: Can't allocate memory for matrices. Aborting... \n\n"); return 1; } printf (" Intializing matrix data \n\n"); for (i = 0; i < (m*k); i++) { A1[i] = A2[i] = (double)(i+1); } for (i = 0; i < (k*n); i++) { B1[i] = B2[i] = (double)(-i-1); } for (i = 0; i < (m*n); i++) { C1[i] = C2[i] = 0.0; C1_fl[i] = C2_fl[i] = 0.0; } printf (" \nComputing matrix product using Intel MKL cblas_dgemm function \n"); t_start = omp_get_wtime(); #pragma omp target data \ map(to: A1[0:m*k], B1[0:k*n], A2[0:m*k], B2[0:k*n]) \ map(tofrom: C1[0:m*n], C2[0:m*n]) { #pragma omp dispatch cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A1, k, B1, n, beta, C1, n); #pragma omp dispatch cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A2, k, B2, n, beta, C2, n); } t_end = omp_get_wtime(); printf ("\n Top left corner of matrix C1: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C1[j+i*n]); } printf ("\n"); } printf ("\n Top left corner of matrix C2: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C2[j+i*n]); } printf ("\n"); } printf (" \nComputing matrix product using for-loops \n"); for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { sum = 0.0; for (q = 0; q < k; q++) sum += A1[k*i+q] * B1[n*q+j]; C1_fl[n*i+j] = sum; } } for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { sum = 0.0; for (q = 0; q < k; q++) sum += A2[k*i+q] * B2[n*q+j]; C2_fl[n*i+j] = sum; } } printf ("\n Top left corner of matrix C1: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C1_fl[j+i*n]); } printf ("\n"); } printf ("\n Top left corner of matrix C2: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C2_fl[j+i*n]); } printf ("\n"); } printf ("\n Computations completed. Verifying... \n\n"); fail = 0; for (i = 0; i < (m*n); i++) { if (! compare(C1[i], C1_fl[i]) || ! compare(C2[i], C2_fl[i])) { fail = 1; break; } } if (fail) { printf (" **** FAIL **** \n"); } else { printf(" time = %lf seconds\n", t_end - t_start); printf (" **** PASS **** \n"); } mkl_free(A1); mkl_free(B1); mkl_free(C1); mkl_free(C1_fl); mkl_free(A2); mkl_free(B2); mkl_free(C2); mkl_free(C2_fl); return 0; }
The two calls to
cblas_dgemm
in the above example can be batched together, resulting in one call to
cblas_dgemm_batch
, as shown in the following example. Note that the batch is composed of one group of size 2, since we have two matrix multiplications with the same set of parameters (
layout
,
transa
,
transb
,
m
,
n
,
k
,
alpha
,
lda
,
ldb
,
beta
, and
ldc
). total_batch_size in this case is 2.
//============================================================== // Copyright © 2022 Intel Corporation // // SPDX-License-Identifier: MIT // ============================================================= // clang-format off #include <stdio.h> #include <stdlib.h> #include <math.h> #include <omp.h> #include "mkl.h" #include "mkl_omp_offload.h" #define min(x,y) (((x) < (y)) ? (x) : (y)) #define epsilon 0.0000001f bool compare(double x, double y) { // return true is x and y are the same return (fabs(x - y) <= epsilon); } int main() { double *A1, *B1, *C1, *C1_fl; double *A2, *B2, *C2, *C2_fl; int m, n, k, i, j, q; double alpha, beta; double sum; int fail; double t_start, t_end; m = 2000, k = 200, n = 1000; alpha = 1.0; beta = 0.0; printf (" Allocating memory for matrices aligned on 64-byte boundary for better \n" " performance \n\n"); A1 = (double *)mkl_malloc (m*k*sizeof( double ), 64 ); B1 = (double *)mkl_malloc (k*n*sizeof( double ), 64 ); C1 = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); C1_fl = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); A2 = (double *)mkl_malloc (m*k*sizeof( double ), 64 ); B2 = (double *)mkl_malloc (k*n*sizeof( double ), 64 ); C2 = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); C2_fl = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); if (A1 == NULL || B1 == NULL || C1 == NULL || C1_fl == NULL || A2 == NULL || B2 == NULL || C2 == NULL || C2_fl == NULL) { printf( "\n ERROR: Can't allocate memory for matrices. Aborting... \n\n"); return 1; } printf (" Intializing matrix data \n\n"); for (i = 0; i < (m*k); i++) { A1[i] = A2[i] = (double)(i+1); } for (i = 0; i < (k*n); i++) { B1[i] = B2[i] = (double)(-i-1); } for (i = 0; i < (m*n); i++) { C1[i] = C2[i] = 0.0; C1_fl[i] = C2_fl[i] = 0.0; } printf (" \nComputing matrix product using Intel MKL cblas_dgemm_batch function \n"); #define GRP_COUNT 1 // 1 group MKL_INT group_count = GRP_COUNT; MKL_INT group_sizes[GRP_COUNT] = {2}; // 8 matrix multiplications CBLAS_TRANSPOSE transa_array[GRP_COUNT] = {CblasNoTrans}; CBLAS_TRANSPOSE transb_array[GRP_COUNT] = {CblasNoTrans}; MKL_INT m_array[GRP_COUNT] = {m}; MKL_INT n_array[GRP_COUNT] = {n}; MKL_INT k_array[GRP_COUNT] = {k}; MKL_INT lda_array[GRP_COUNT] = {k}; MKL_INT ldb_array[GRP_COUNT] = {n}; MKL_INT ldc_array[GRP_COUNT] = {n}; double alpha_array[GRP_COUNT] = {alpha}; double beta_array[GRP_COUNT] = {beta}; // Number of matrix multiplications = 2 double **a_array, **b_array, **c_array; a_array = (double **)mkl_calloc(2, sizeof( double* ), 64); b_array = (double **)mkl_calloc(2, sizeof( double* ), 64); c_array = (double **)mkl_calloc(2, sizeof( double* ), 64); t_start = omp_get_wtime(); // Call cblas_dgemm_batch #pragma omp target enter data \ map(to: A1[0:m*k], B1[0:k*n], C1[0:m*n]) \ map(to: A2[0:m*k], B2[0:k*n], C2[0:m*n]) #pragma omp target data use_device_ptr(A1, B1, C1, A2, B2, C2) { a_array[0] = A1, a_array[1] = A2; b_array[0] = B1, b_array[1] = B2; c_array[0] = C1, c_array[1] = C2; } #pragma omp target data \ map(to:a_array[0:2], b_array[0:2], c_array[0:2]) { #pragma omp dispatch cblas_dgemm_batch ( CblasRowMajor, transa_array, transb_array, m_array, n_array, k_array, alpha_array, (const double **)a_array, lda_array, (const double **)b_array, ldb_array, beta_array, c_array, ldc_array, group_count, group_sizes); } // end target data map #pragma omp target exit data \ map(from: C1[0:m*n], C2[0:m*n]) t_end = omp_get_wtime(); printf ("\n Top left corner of matrix C1: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C1[j+i*n]); } printf ("\n"); } printf ("\n Top left corner of matrix C2: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C2[j+i*n]); } printf ("\n"); } printf (" \nComputing matrix product using for-loops \n"); for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { sum = 0.0; for (q = 0; q < k; q++) sum += A1[k*i+q] * B1[n*q+j]; C1_fl[n*i+j] = sum; } } for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { sum = 0.0; for (q = 0; q < k; q++) sum += A2[k*i+q] * B2[n*q+j]; C2_fl[n*i+j] = sum; } } printf ("\n Top left corner of matrix C1: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C1_fl[j+i*n]); } printf ("\n"); } printf ("\n Top left corner of matrix C2: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C2_fl[j+i*n]); } printf ("\n"); } printf ("\n Computations completed. Verifying... \n\n"); fail = 0; for (i = 0; i < (m*n); i++) { if (! compare(C1[i], C1_fl[i]) || ! compare(C2[i], C2_fl[i])) { fail = 1; break; } } if (fail) { printf (" **** FAIL **** \n"); } else { printf(" time = %lf seconds\n", t_end - t_start); printf (" **** PASS **** \n"); } mkl_free(A1); mkl_free(B1); mkl_free(C1); mkl_free(C1_fl); mkl_free(A2); mkl_free(B2); mkl_free(C2); mkl_free(C2_fl); return 0; }
The performance of the above two examples when running on the particular GPU used (1-tile only) was as follows:
dgemm_example_01_c.cpp (two calls to cblas_dgemm): 2.976183 seconds dgemm_batch_example_01_c.cpp (one call to cblas_dgemm_batch): 1.881641 seconds
A more complex example of batching is shown below. In this example, we have a batch composed of 3 groups (GROUP_COUNT=3). The size of each group is a randomly chosen number between 1 and 10. Several parameters (
layout
,
transA
,
transB
,
m
,
n
, and
k
) are chosen randomly, but in each group the parameters are the same for all the multiplications. The total_batch_size is equal to the sum of all the group sizes.
//============================================================== // Copyright © 2022 Intel Corporation // // SPDX-License-Identifier: MIT // ============================================================= // clang-format off /******************************************************************************* * Copyright 2019-2021 Intel Corporation. * * This software and the related documents are Intel copyrighted materials, and * your use of them is governed by the express license under which they were * provided to you (License). Unless the License provides otherwise, you may not * use, modify, copy, publish, distribute, disclose or transmit this software or * the related documents without Intel's prior written permission. * * This software and the related documents are provided as is, with no express * or implied warranties, other than those that are expressly stated in the * License. *******************************************************************************/ #include <stdio.h> #include <omp.h> #include "mkl.h" #include "mkl_omp_offload.h" #include "common.h" #define GROUP_COUNT 3 int dnum = 0; int main() { CBLAS_LAYOUT layout = (rand_int(0,1) == 0) ? CblasColMajor : CblasRowMajor; CBLAS_TRANSPOSE *transA, *transB; MKL_INT *m, *n, *k, *lda, *ldb, *ldc; double *alpha, *beta; MKL_INT *group_size, *sizea_array, *sizeb_array, *sizec_array, total_batch_size = 0, sizea, sizeb, sizec; double **a_array, **b_array, **c_array, **c_ref_array; double **a_array_dev, **b_array_dev, **c_array_dev; transA = (CBLAS_TRANSPOSE *)mkl_malloc(GROUP_COUNT * sizeof(CBLAS_TRANSPOSE), 64); transB = (CBLAS_TRANSPOSE *)mkl_malloc(GROUP_COUNT * sizeof(CBLAS_TRANSPOSE), 64); m = (MKL_INT *)mkl_malloc(GROUP_COUNT * sizeof(MKL_INT), 64); n = (MKL_INT *)mkl_malloc(GROUP_COUNT * sizeof(MKL_INT), 64); k = (MKL_INT *)mkl_malloc(GROUP_COUNT * sizeof(MKL_INT), 64); lda = (MKL_INT *)mkl_malloc(GROUP_COUNT * sizeof(MKL_INT), 64); ldb = (MKL_INT *)mkl_malloc(GROUP_COUNT * sizeof(MKL_INT), 64); ldc = (MKL_INT *)mkl_malloc(GROUP_COUNT * sizeof(MKL_INT), 64); group_size = (MKL_INT *)mkl_malloc(GROUP_COUNT * sizeof(MKL_INT), 64); alpha = (double *)mkl_malloc(GROUP_COUNT * sizeof(double), 64); beta = (double *)mkl_malloc(GROUP_COUNT * sizeof(double), 64); if ((m == NULL) || (n == NULL) || (k == NULL) || (lda == NULL) || (ldb == NULL) || (ldc == NULL) || (group_size == NULL) || (alpha == NULL) || (beta == NULL)) { printf("Cannot allocate input arrays\n"); return 1; } MKL_INT i, j, p, idx; for (i = 0; i < GROUP_COUNT; i++) { transA[i] = (rand_int(0,1) == 0) ? CblasNoTrans : CblasTrans; transB[i] = (rand_int(0,1) == 0) ? CblasNoTrans : CblasTrans; alpha[i] = rand_double_scalar(); beta[i] = rand_double_scalar(); m[i] = rand_int(1, 20); n[i] = rand_int(1, 20); k[i] = rand_int(1, 20); lda[i] = MAX(m[i], k[i]); ldb[i] = MAX(k[i], n[i]); ldc[i] = MAX(m[i], n[i]); group_size[i] = rand_int(1, 10); total_batch_size += group_size[i]; #ifdef MKL_ILP64 printf("Group %lld: layout = %s, transA = %s, transB = %s, m = %lld, n = %lld, k = %lld, lda = %lld, ldb = %lld, ldc = %lld, alpha = %lf, beta = %lf, group_size = %lld\n", i, (layout == CblasColMajor) ? "Column Major" : "Row Major", (transA[i] == CblasNoTrans) ? "Non Transpose" : "Transpose", (transB[i] == CblasNoTrans) ? "Non Transpose" : "Transpose", m[i], n[i], k[i], lda[i], ldb[i], ldc[i], alpha[i], beta[i], group_size[i]); #else printf("Group %d: layout = %s, transA = %s, transB = %s, m = %d, n = %d, k = %d, lda = %d, ldb = %d, ldc = %d, alpha = %lf, beta = %lf, group_size = %d\n", i, (layout == CblasColMajor) ? "Column Major" : "Row Major", (transA[i] == CblasNoTrans) ? "Non Transpose" : "Transpose", (transB[i] == CblasNoTrans) ? "Non Transpose" : "Transpose", m[i], n[i], k[i], lda[i], ldb[i], ldc[i], alpha[i], beta[i], group_size[i]); #endif } sizea_array = (MKL_INT *)mkl_malloc(sizeof(MKL_INT) * total_batch_size, 64); sizeb_array = (MKL_INT *)mkl_malloc(sizeof(MKL_INT) * total_batch_size, 64); sizec_array = (MKL_INT *)mkl_malloc(sizeof(MKL_INT) * total_batch_size, 64); a_array = (double **)mkl_malloc(sizeof(double *) * total_batch_size, 64); b_array = (double **)mkl_malloc(sizeof(double *) * total_batch_size, 64); c_array = (double **)mkl_malloc(sizeof(double *) * total_batch_size, 64); a_array_dev = (double **)mkl_malloc(sizeof(double *) * total_batch_size, 64); b_array_dev = (double **)mkl_malloc(sizeof(double *) * total_batch_size, 64); c_array_dev = (double **)mkl_malloc(sizeof(double *) * total_batch_size, 64); c_ref_array = (double **)mkl_malloc(sizeof(double *) * total_batch_size, 64); if ((sizea_array == NULL) || (sizeb_array == NULL) || (sizec_array == NULL) || (a_array == NULL) || (b_array == NULL) || (c_array == NULL) || (a_array_dev == NULL) || (b_array_dev == NULL) || (c_array_dev == NULL) || (c_ref_array == NULL)) { printf("Cannot allocate matrices and size arrays\n"); return 1; } idx = 0; for (i = 0; i < GROUP_COUNT; i++) { sizea = (((layout == CblasRowMajor) && (transA[i] == CblasTrans)) || ((layout == CblasColMajor) && (transA[i] == CblasNoTrans))) ? lda[i] * k[i] : m[i] * lda[i]; sizeb = (((layout == CblasRowMajor) && (transB[i] == CblasTrans)) || ((layout == CblasColMajor) && (transB[i] == CblasNoTrans))) ? ldb[i] * n[i] : k[i] * ldb[i]; sizec = (layout == CblasColMajor) ? ldc[i] * n[i] : ldc[i] * m[i]; for (j = 0; j < group_size[i]; j++) { a_array[idx] = (double *)mkl_malloc(sizeof(double) * sizea, 64); a_array_dev[idx] = a_array[idx]; sizea_array[idx] = sizea; if (a_array[idx] == NULL) { printf("cannot allocate a matrices\n"); return 1; } b_array[idx] = (double *)mkl_malloc(sizeof(double) * sizeb, 64); b_array_dev[idx] = b_array[idx]; sizeb_array[idx] = sizeb; if (b_array[idx] == NULL) { printf("cannot allocate b matrices\n"); return 1; } c_array[idx] = (double *)mkl_malloc(sizeof(double) * sizec, 64); c_array_dev[idx] = c_array[idx]; sizec_array[idx] = sizec; if (c_array[idx] == NULL) { printf("cannot allocate c matrices\n"); return 1; } c_ref_array[idx] = (double *)mkl_malloc(sizeof(double) * sizec, 64); if (c_ref_array[idx] == NULL) { printf("cannot allocate c_ref matrices\n"); return 1; } init_double_array(sizea, a_array[idx], 1); init_double_array(sizeb, b_array[idx], 1); init_double_array(sizec, c_array[idx], 1); for (p = 0; p < sizec_array[idx]; p++) c_ref_array[idx][p] = c_array[idx][p]; idx++; } } // run gemm_batch on host, use standard oneMKL interface cblas_dgemm_batch(layout, transA, transB, m, n, k, alpha, (const double **) a_array, lda, (const double **) b_array, ldb, beta, c_ref_array, ldc, GROUP_COUNT, group_size); double *a, *b, *c; for (i = 0; i < total_batch_size; i++) { a = a_array[i]; b = b_array[i]; c = c_array[i]; #pragma omp target enter data map(to:a[0:sizea_array[i]],b[0:sizeb_array[i]],c[0:sizec_array[i]]) #pragma omp target data use_device_ptr(a,b,c) { a_array_dev[i] = a; b_array_dev[i] = b; c_array_dev[i] = c; } } #pragma omp target data map(to:a_array_dev[0:total_batch_size], \ b_array_dev[0:total_batch_size], \ c_array_dev[0:total_batch_size]) device(dnum) { #pragma omp dispatch cblas_dgemm_batch(layout, transA, transB, m, n, k, alpha, (const double **) a_array_dev, lda, (const double **) b_array_dev, ldb, beta, c_array_dev, ldc, GROUP_COUNT, group_size); } for (i = 0; i < total_batch_size; i++) { a = a_array[i]; b = b_array[i]; c = c_array[i]; #pragma omp target exit data map(from:a[0:sizea_array[i]],b[0:sizeb_array[i]],c[0:sizec_array[i]]) } double computed, reference, diff; MKL_INT l; idx = 0; for (p = 0; p < GROUP_COUNT; p++) { for (l = 0; l < group_size[p]; l++) { for (i = 0; i < m[p]; i++) { for (j = 0; j < n[p]; j++) { if (layout == CblasColMajor) { computed = c_array[idx][i + j * ldc[p]]; reference = c_ref_array[idx][i + j * ldc[p]]; } else { computed = c_array[idx][j + i * ldc[p]]; reference = c_ref_array[idx][j + i * ldc[p]]; } diff = computed - reference; diff = (diff > 0) ? diff : -diff; if (diff > 0.0001) { #ifdef MKL_ILP64 printf("Error in matrix %lld (group = %lld, matrix index in group = %lld) at index [%lld][%lld], computed = %lf, reference = %lf, difference = %lf\n", idx, p, l, i, j, computed, reference, diff); #else printf("Error in matrix %d at index [%d][%d], computed = %lf, reference = %lf, difference = %lf\n", idx, i, j, computed, reference, diff); #endif free_double_matrices(a_array, total_batch_size); free_double_matrices(b_array, total_batch_size); free_double_matrices(c_array, total_batch_size); free_double_matrices(c_ref_array, total_batch_size); mkl_free(a_array); mkl_free(b_array); mkl_free(c_array); mkl_free(c_ref_array); mkl_free(a_array_dev); mkl_free(b_array_dev); mkl_free(c_array_dev); mkl_free(sizea_array); mkl_free(sizeb_array); mkl_free(sizec_array); mkl_free(transA); mkl_free(transB); mkl_free(m); mkl_free(n); mkl_free(k); mkl_free(lda); mkl_free(ldb); mkl_free(ldc); mkl_free(group_size); mkl_free(alpha); mkl_free(beta); return 1; } } } idx++; } } printf("Validation PASSED\n"); free_double_matrices(a_array, total_batch_size); free_double_matrices(b_array, total_batch_size); free_double_matrices(c_array, total_batch_size); free_double_matrices(c_ref_array, total_batch_size); mkl_free(a_array); mkl_free(b_array); mkl_free(c_array); mkl_free(c_ref_array); mkl_free(a_array_dev); mkl_free(b_array_dev); mkl_free(c_array_dev); mkl_free(sizea_array); mkl_free(sizeb_array); mkl_free(sizec_array); mkl_free(transA); mkl_free(transB); mkl_free(m); mkl_free(n); mkl_free(k); mkl_free(lda); mkl_free(ldb); mkl_free(ldc); mkl_free(group_size); mkl_free(alpha); mkl_free(beta); return 0; }

Speeding Up Independent, Consecutive GEMM Calls

There are various ways to speed up the execution of consecutive GEMM calls that can be executed independently. One way is to batch the GEMM calls by calling the batch version of GEMM as shown above.
Another way is to enclose the calls to GEMM by an OpenMP
parallel
construct, so each OpenMP thread executing the parallel region dispatches one of the GEMM calls. This
parallel
approach is illustrated in the following example.
//============================================================== // Copyright © 2022 Intel Corporation // // SPDX-License-Identifier: MIT // ============================================================= // clang-format off #include <stdio.h> #include <stdlib.h> #include <math.h> #include <omp.h> #include "mkl.h" #include "mkl_omp_offload.h" #define min(x,y) (((x) < (y)) ? (x) : (y)) #define epsilon 0.0000001f bool compare(double x, double y) { // returns true if x and y are the same return fabs(x - y) <= epsilon; } int main() { double *A1, *B1, *C1, *C1_fl; double *A2, *B2, *C2, *C2_fl; int m, n, k, i, j, q; double alpha, beta; double sum; int fail; double t_start, t_end; m = 2000, k = 200, n = 1000; alpha = 1.0; beta = 0.0; printf (" Allocating memory for matrices aligned on 64-byte boundary for better \n" " performance \n\n"); A1 = (double *)mkl_malloc (m*k*sizeof( double ), 64 ); B1 = (double *)mkl_malloc (k*n*sizeof( double ), 64 ); C1 = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); C1_fl = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); A2 = (double *)mkl_malloc (m*k*sizeof( double ), 64 ); B2 = (double *)mkl_malloc (k*n*sizeof( double ), 64 ); C2 = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); C2_fl = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); if (A1 == NULL || B1 == NULL || C1 == NULL || C1_fl == NULL || A2 == NULL || B2 == NULL || C2 == NULL || C2_fl == NULL) { printf( "\n ERROR: Can't allocate memory for matrices. Aborting... \n\n"); return 1; } printf (" Intializing matrix data \n\n"); for (i = 0; i < (m*k); i++) { A1[i] = A2[i] = (double)(i+1); } for (i = 0; i < (k*n); i++) { B1[i] = B2[i] = (double)(-i-1); } for (i = 0; i < (m*n); i++) { C1[i] = C2[i] = 0.0; C1_fl[i] = C2_fl[i] = 0.0; } printf (" \nComputing matrix product using Intel MKL cblas_dgemm function \n"); t_start = omp_get_wtime(); #pragma omp target data \ map(to: A1[0:m*k], B1[0:k*n], A2[0:m*k], B2[0:k*n]) \ map(tofrom: C1[0:m*n], C2[0:m*n]) { #pragma omp parallel num_threads(2) { int id = omp_get_thread_num(); if (id == 0) { #pragma omp dispatch cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A1, k, B1, n, beta, C1, n); } else if (id == 1) { #pragma omp dispatch cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A2, k, B2, n, beta, C2, n); } } } t_end = omp_get_wtime(); printf ("\n Top left corner of matrix C1: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C1[j+i*n]); } printf ("\n"); } printf ("\n Top left corner of matrix C2: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C2[j+i*n]); } printf ("\n"); } printf (" \nComputing matrix product using for-loops \n"); for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { sum = 0.0; for (q = 0; q < k; q++) sum += A1[k*i+q] * B1[n*q+j]; C1_fl[n*i+j] = sum; } } for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { sum = 0.0; for (q = 0; q < k; q++) sum += A2[k*i+q] * B2[n*q+j]; C2_fl[n*i+j] = sum; } } printf ("\n Top left corner of matrix C1: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C1_fl[j+i*n]); } printf ("\n"); } printf ("\n Top left corner of matrix C2: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C2_fl[j+i*n]); } printf ("\n"); } printf ("\n Computations completed. Verifying... \n\n"); fail = 0; for (i = 0; i < (m*n); i++) { if (! compare(C1[i], C1_fl[i]) || ! compare(C2[i], C2_fl[i])) { fail = 1; break; } } if (fail) { printf (" **** FAIL **** \n"); } else { printf(" time = %lf seconds\n", t_end - t_start); printf (" **** PASS **** \n"); } mkl_free(A1); mkl_free(B1); mkl_free(C1); mkl_free(C1_fl); mkl_free(A2); mkl_free(B2); mkl_free(C2); mkl_free(C2_fl); return 0; }
Yet another way to speed up the execution of independent, consecutive GEMM calls is to use the
nowait
clause on the
dispatch
construct so the host thread does not have to wait for a dispatched GEMM call to complete before dispatching the next one. After the last GEMM call, we insert an OpenMP
taskwait
directive to guarantee that all the dispatched MKL calls complete before the host thread proceeds any further. This
nowait
approach is illustrated in the following example.
//============================================================== // Copyright © 2022 Intel Corporation // // SPDX-License-Identifier: MIT // ============================================================= // clang-format off #include <stdio.h> #include <stdlib.h> #include <math.h> #include <omp.h> #include "mkl.h" #include "mkl_omp_offload.h" #define min(x,y) (((x) < (y)) ? (x) : (y)) #define epsilon 0.0000001f bool compare(double x, double y) { // returns true if x and y are the same return fabs(x - y) <= epsilon; } int main() { double *A1, *B1, *C1, *C1_fl; double *A2, *B2, *C2, *C2_fl; int m, n, k, i, j, q; double alpha, beta; double sum; int fail; double t_start, t_end; m = 2000, k = 200, n = 1000; alpha = 1.0; beta = 0.0; printf (" Allocating memory for matrices aligned on 64-byte boundary for better \n" " performance \n\n"); A1 = (double *)mkl_malloc (m*k*sizeof( double ), 64 ); B1 = (double *)mkl_malloc (k*n*sizeof( double ), 64 ); C1 = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); C1_fl = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); A2 = (double *)mkl_malloc (m*k*sizeof( double ), 64 ); B2 = (double *)mkl_malloc (k*n*sizeof( double ), 64 ); C2 = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); C2_fl = (double *)mkl_malloc (m*n*sizeof( double ), 64 ); if (A1 == NULL || B1 == NULL || C1 == NULL || C1_fl == NULL || A2 == NULL || B2 == NULL || C2 == NULL || C2_fl == NULL) { printf( "\n ERROR: Can't allocate memory for matrices. Aborting... \n\n"); return 1; } printf (" Intializing matrix data \n\n"); for (i = 0; i < (m*k); i++) { A1[i] = A2[i] = (double)(i+1); } for (i = 0; i < (k*n); i++) { B1[i] = B2[i] = (double)(-i-1); } for (i = 0; i < (m*n); i++) { C1[i] = C2[i] = 0.0; C1_fl[i] = C2_fl[i] = 0.0; } printf (" \nComputing matrix product using Intel MKL cblas_dgemm function \n"); t_start = omp_get_wtime(); #pragma omp target data \ map(to: A1[0:m*k], B1[0:k*n], A2[0:m*k], B2[0:k*n]) \ map(tofrom: C1[0:m*n], C2[0:m*n]) { #pragma omp dispatch nowait cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A1, k, B1, n, beta, C1, n); #pragma omp dispatch nowait cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A2, k, B2, n, beta, C2, n); #pragma omp taskwait } t_end = omp_get_wtime(); printf ("\n Top left corner of matrix C1: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C1[j+i*n]); } printf ("\n"); } printf ("\n Top left corner of matrix C2: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C2[j+i*n]); } printf ("\n"); } printf (" \nComputing matrix product using for-loops \n"); for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { sum = 0.0; for (q = 0; q < k; q++) sum += A1[k*i+q] * B1[n*q+j]; C1_fl[n*i+j] = sum; } } for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { sum = 0.0; for (q = 0; q < k; q++) sum += A2[k*i+q] * B2[n*q+j]; C2_fl[n*i+j] = sum; } } printf ("\n Top left corner of matrix C1: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C1_fl[j+i*n]); } printf ("\n"); } printf ("\n Top left corner of matrix C2: \n"); for (i=0; i<min(m,6); i++) { for (j=0; j<min(n,6); j++) { printf ("%12.5G", C2_fl[j+i*n]); } printf ("\n"); } printf ("\n Computations completed. Verifying... \n\n"); fail = 0; for (i = 0; i < (m*n); i++) { if (! compare(C1[i], C1_fl[i]) || ! compare(C2[i], C2_fl[i])) { fail = 1; break; } } if (fail) { printf (" **** FAIL **** \n"); } else { printf(" time = %lf seconds\n", t_end - t_start); printf (" **** PASS **** \n"); } mkl_free(A1); mkl_free(B1); mkl_free(C1); mkl_free(C1_fl); mkl_free(A2); mkl_free(B2); mkl_free(C2); mkl_free(C2_fl); return 0; }

Product and Performance Information

1

Performance varies by use, configuration and other factors. Learn more at www.Intel.com/PerformanceIndex.