Technical Report 2014-001

Implementing Level-3 BLAS Routines in OpenCL
on Diﬀerent Processing Units

Kazuya Matsumoto, Naohito Nakasato, and Stanislav Sedukhin
October 22, 2014

Graduate School of Computer Science and Engineering
The University of Aizu
Tsuruga, Ikki-Machi, Aizu-Wakamatsu City
Fukushima, 965-8580 Japan

Technical Report 2014-001
Title:
Implementing Level-3 BLAS Routines in OpenCL on Diﬀerent Processing Units

Authors:
Kazuya Matsumoto, Naohito Nakasato, and Stanislav Sedukhin
Key Words and Phrases:
Level-3 BLAS, GPU, multi-core CPU, many-core processor, OpenCL, performance porting,auto-tuning

Abstract:
This paper presents an implementation of diﬀerent matrix-matrix multiplication routines in OpenCL.
We utilize the high-performance GEMM (GEneral Matrix-Matrix Multiply) implemen- tation fromour previous work for the present implementation of other matrix-matrix multiply routines in Level-3
BLAS (Basic Linear Algebra Subprograms). The other routines include SYMM (Symmetric MatrixMatrix Multiply), SYRK (Symmetric Rank-K Update), SYR2K (Symmetric Rank-2K Update), andTRMM (Triangular Matrix-Matrix Multiply). A key in our approach is to copy given matrix databy copying OpenCL kernels into a form such that a high-performance GEMM kernel can be utilized
for computation. We use a previously developed auto-tuning system for the highly optimized copying
kernels as well as for GEMM kernel. The performance evaluation of our implementation is conducted onfour diﬀerent GPUs (AMD Radeon R9 290X, FirePro W9100, Radeon HD 7970, and NVIDIA GeForceGTX Titan), a many-core processor (Intel Xeon Phi 5110P), and a multi-core processor (Core i7
3960X). The evaluation results show that the tuning on the copying kernels is eﬀective and contributes
to develop high-performance BLAS3 routines.

Report Date:
10/22/2014

Written Language:
EnglishAny Other Identifying Information of this Report:

Distribution Statement:
First Issue: 10 copies

Supplementary Notes:

Distributed Parallel Processing Laboratory
The University of Aizu
Aizu-Wakamatsu
Fukushima 965-8580
Japan

1

Introduction

As an interface of numerical software library, Basic Linear Algebra Subprograms (BLAS)
is widely used in scientiﬁc and engineering communities. The BLAS consists of threelevels. The Level-1 BLAS is a collection of vector-vector routines, the Level-2 BLAS is
dedicated to matrix-vector routines, and the Level-3 BLAS (BLAS3) includes matrixmatrix routines. Processor manufactures, universities, and research institutions havestruggled to develop highly-tuned BLAS libraries since the publication of BLAS standard [9, 4].
The General Matrix-Matrix Multiply (GEMM) routine is considered to be the mostfundamental BLAS3 routine. Utilizing a tuned GEMM routine as the core computation is generally accepted approach to implement other BLAS3 matrix-matrix routines by considering each of the other routines as a combination of GEMM and asmall amount of Level-1 and Level-2 BLAS computations [9]. Igual et al. [6, 5] presented high-performance BLAS3 based on cuBLAS library. They made experiments on
three diﬀerent GEMM-based algorithms for each BLAS3 routine. Their implementation
showed better performance without touching underlining CUDA kernels than simply
using BLAS3 routine of cuBLAS.
Auto-tuning is an important technique to resolve the problem of performance portability, and it is a well accepted solution for developing high-performance BLAS implementation. ATLAS (Automatically Tuned Linear Algebra Software) [15] is famous
projects for auto-tuned BLAS routines on CPUs. Also, several auto-tuning systems for
GEMM have been developed [2, 7, 12, 10, 3, 11]. In CUDA, an auto-tuning framework
for NVIDIA GPUs has been implemented [10]. In OpenCL, Du et al. [3] presented
auto-tuned GEMM routines for GPUs. In our previous work [11], we have also implemented an auto-tuning system for a high-performance GEMM kernel in OpenCL. The
GEMM routine is realized with an approach that copies all matrix data into buﬀers andthen runs the auto-tuned GEMM kernel. This present work extends the approach withthe matrix data copying for other BLAS3 routines including Symmetric Matrix-Matrix
multiply (SYMM), Symmetric Rank-K update (SYRK), Symmetric Rank-2K update(SYR2K), and Triangular Matrix-Matrix multiply (TRMM). We consider that this approach is eﬀective to develop highly eﬃcient BLAS3 routines especially on many-coreprocessors like GPUs.
Another contribution of this work is that we use the developed auto-tuning system
to tune copying kernels as well as GEMM kernel. Copying kernels are also an important component in our implementation. Particularly for small problem sizes, the timespending on the data copying is a performance bottleneck in our BLAS3 implementation. We show that the auto-tuning makes it possible to develop high-performancecopying kernels.
In this work, our BLAS3 is implemented in OpenCL. OpenCL is a standard framework for parallel programming [8, 1]. Programs written in OpenCL are functionallyportable across various processors that include CPUs, GPUs, and other computing devices such as Intel Xeon Phi co-processor and FPGA. This paper presents results of3

performance evaluation on the following diﬀerent processing units:
1. AMD Radeon R9 290X (Volcanic Islands);
2. AMD FirePro W9100 (Volcanic Islands);
3. AMD Radeon HD 7970 (Southern Islands);
4. NVIDIA GeForce GTX Titan (Kepler);
5. Intel Xeon Phi 5110P (MIC);
6. Intel Core i7 3960X (Sandy Bridge).
The rest of this paper is organized as follows: Section 2 presents the basic information on OpenCL. Section 3 describes our approach for high-performance BLAS3
implementation. Section 4 shows the performance evaluation of the BLAS3 implementation and provides comparison with the existing vendor provided BLAS libraries.
Finally, Section 6 gives the concluding remarks of this study.

2

OpenCL Basics

OpenCL (Open Computing Language) is an open standard for general purpose parallel
programming on heterogeneous computing platforms. The OpenCL framework includes
a C99-based programming language for writing parallel functions called kernels, andruntime APIs for controlling OpenCL platforms and devices.
An OpenCL platform consists of one or more OpenCL devices. An OpenCL device
comprises multiple compute units (CUs), each of which has multiple processing elements(PEs). When an OpenCL kernel is executed on the device, an N-dimensional index
space, which is called NDRange, is deﬁned. In this work, we consider a two-dimensional
index space only, which is suitable for dealing with matrix data. Each instance in an
NDRange, is called a work-item. Every work-item has a unique ID. Several work-items
organize a work-group. A work-item runs on one or more PEs. A task on a work-groupis processed by all the PEs of a CU.
In an OpenCL kernel, four distinct memory regions are accessible to work-items.
1. Global memory (gm) is a memory region in which data can be read/written by all
work-items. It is impossible to synchronize work-items during kernel execution inthis memory.
2. Constant memory (cm) is a read-only region of global memory. Data in thisregion are not changed during execution.
3. Local memory (lm) is a speciﬁc memory region to a work-group. Work-items ina work-group can share data in local memory.

44. Private memory (pm) is a speciﬁc memory region to a work-item. Data in private
memory of a work-item is not visible to other work-items. On most OpenCL
devices, private memory is in the register ﬁle.

3

Implementation

We use a common strategy in our implementation of diﬀerent matrix-matrix multiply
routines. The strategy is to copy matrix data into a storage layout such that a highperformance GEMM kernel can be utilized for computation.

3.1

GEMM

Let us review the GEMM implementation presented in our previous paper [11]. We
focus on the matrix multiply-add C ← C + AB operation. The current implementation
supposes column-major as the initial matrix storage order for consistency with BLAS
standard1 .
Our GEMM routine is implemented with a GEMM (AB T + C) kernel and a data
copying kernel in OpenCL. The copying kernel is not simply replicating matrix data into
additionally allocated memory space (buﬀer). For computing AB + C, transposing thematrix B is necessary to utilize the auto-tuned AB T + C kernel. Moreover, the tuned
GEMM kernel assumes that matrix data are aligned in a block-major order for highperformance, and zeros are padded if a matrix dimension is not in multiples of blocking
factor. In summary, the copying kernel conducts matrix transposition, storage layoutchanging and zero padding if necessary. For simplicity, the following descriptions assume
cases where matrices are square, all matrix dimensions are multiples of any blocking
factors, and matrix data are aligned in column-major order if not speciﬁed.
The copying kernel as well as GEMM kernel uses two levels of blocking for eﬃciently
utilizing caches of processing units. In this paper, we describe values on the ﬁrst levelof blocking for work-item processing as Mwi , Nwi , Kwi and values on the other level ofblocking for work-group processing as Mwg , Nwg , Kwg . Fig. 1 delineates a kernel which
copies matrix A data in column-major to a buﬀer in a form of column-block columnmajor (CBC) order for our GEMM routine. The matrix-matrix multiply operation partis generally represented as GEMM BODY in Fig. 2. After copying matrices A and B tobuﬀer, our GEMM routine calls the GEMM kernel (see Fig. 3).

3.2

SYMM

For SYMM routine, we focus on C ← AB + C, where A is a symmetric matrix and only
the lower triangular part stores the matrix data. The diﬀerence between GEMM andSYMM is only in an operation on the matrix A: our SYMM routine are implemented byonly changing an OpenCL kernel which copies the matrix A from the GEMM routine.
1

In our previous work, the GEMM OpenCL kernel fundamentally supposes row-major order.

5

k
Kwg
Kwi
0

8 16 24 32 40 48 56

04

8 12 32 36 40 44

1

9 17 25 33 41 49 57

1

5

9 13 33 37 41 45

2 10 18 26 34 42 50 58

2

6 10 14 34 38 42 46

3

7 11 15 35 39 43 47

3 11 19 27 35 43 51 59
4 12 20 28 36 44 52 60

Copy

16 20 24 28 48 52 56 60

5 13 21 29 37 45 53 61

17 21 25 29 49 53 57 61

6 14 22 30 38 46 54 62

18 22 26 30 50 54 58 62

7 15 23 31 39 47 55 63

19 23 27 31 51 55 59 63

Mwi
Mwg
m

Figure 1: Copying matrix data from column-major order into a column-block columnmajor order with blocking factors Mwg , Kwg , Mwi , Kwi
GEMM BODY(A, B, Start, End)
1: Cpm := 0
2: for p = Start to End − Kwi step Kwi do
3:Load Mwi · Kwi elements of A from global memory into private memory (Apm )4:
Load Kwi · Nwi elements of B from global memory into private memory (Bpm )5:
Cpm + = Apm × Bpm
6: end for
7: Return Cpm
Figure 2: Body of GEMM kernel
GEMM KERNEL(m, n, k, A, B, C)
1: Cpm := GEMM BODY(A, B, 0, k)
2: // Accumulating the product Cpm to the corresponding Mwi · Nwi data in global
memory (Cgm )3: Cgm := Cpm + Cgm
Figure 3: GEMM kernel algorithm
The newly prepared kernel copies data of A into a buﬀer A0 in a form of symmetricmatrix. This means that the lower triangular and the diagonal parts of A0 is ﬁlled withthe lower triangular matrix A straightforwardly and the upper triangular part of A0 is
ﬁlled with transposed data of the lower triangular matrix A:{
a(i, p) if i ≥ p;0
a (i, p) =a(p, i) if i < p,6

where a(i, p) and a0 (i, p) represents the element of A and A0 in (i, p) position, respectively. Since we can treat SYMM as GEMM after the copying, the SYMM routine
simply calls the GEMM KERNEL.

3.3

SYRK

For SYRK routine, we select C ← C + AAT , where only the lower triangular part of
C is stored and updated. Because of the lower triangular property, we are required toprepare a diﬀerent matrix multiplication kernel for SYRK based on the GEMM kernel.
Fig. 4 gives the SYRK kernel algorithm. The auto-tuned GEMM kernel (GEMM BODY
in Fig. 2) is called inside the SYRK algorithm.
The newly prepared SYRK kernel only updates the lower triangular part. To realize the operation, an work-group containing global work-item ID (y, x) updates the
responsible Mwg · Nwg elements of a block if the block is in lower triangular part ofC, otherwise, it returns immediately. All the work-items of a work-group assigned for
updating a diagonal block are responsible to compute the block in a parallel fashion
redundantly to avoid conditional statements which yield large latencies.
SYRK KERNEL(n, k, A, B, C)
1: Get a two-dimensional global ID (y, x)
2: if y · (Mwi /Mwg + 1) · Mwi < x · Nwi /Nwg then3:
Return // Assigned block for (y, x) is not in the lower triangular part
4: end if
5: Cpm := GEMM BODY(A, B, 0, k)
6: Merge Cpm with Mwi · Nwi elements only of the lower triangular part of C
Figure 4: SYRK (C := C + AAT ) kernel algorithm

3.4

SYR2K

For SYR2K routine, we focus on C := C + AB T + BAT , where only the lower triangularpart of C is stored and updated. Our SYR2K implementation simply calls the aboveSYRK routine twice, i.e., C := C+AB T and C := C+BAT computations are conductedin a consecutive order.

3.5

TRMM

For TRMM routine, we focus on B := AB, where A is a upper triangular matrix. Like
the SYRK case, our TRMM OpenCL kernel is a slightly diﬀerent from the GEMM
kernel. Since a multiplication of the upper triangular part of A and the general matrix
B is necessary, the TRMM kernel shifts the starting indexes of the multiplication ina work-group in units of the work-group blocking factor Kwg . Fig. 5 schematically7

Figure 5: TRMM (B := AB) scheme
illustrates our TRMM implementation. The following statement shows the operations:
BIJ :=K−1∑

AIP BP J(0 ≤ I, J ≤ m/Mwg − 1 and K = m/Kwg ).

P =I

In addition to the TRMM kernel, we have prepared a copying kernel to make the
TRMM computation possible. Like in the SYMM case, the prepared kernel copies dataof A into a buﬀer A0 . This indicates that data elements of the upper triangular anddiagonal parts of A0 are ﬁlled with the corresponding data of A, and all elements of thelower triangular part of A0 are initialized with zero:{a(i, p) if i ≤ p;0
a (i, p) =0
if i > p.

3.6

Performance Tuning of Copying Kernels

As mentioned above, several copying kernels are required in addition to the matrixmatrix multiplication kernels. The amount of matrix data copying is O(n2 ) for generalsquare (n × n) matrices while the amount of the computations (multiply-adds) is O(n3 )in case m = n = k. The larger the matrix size is, the copying time becomes less criticalin the total computation time. Nevertheless, we try to tune copying kernels as well to
avoid using these kernels whose performance is low.
Among BLAS3 routines implemented in this work, GEMM, SYMM, and TRMM
routines require respectively diﬀerent copying kernels. We have newly implemented
code generators which produce copying kernels for the auto-tuning. Each code generator
takes a set of six parameters which are two work-group blocking factors Mwg , Nwg , twowork-item blocking factors Mwi , Nwi , and two loop unrolling degrees Mud , Nud . Fig. 6
represents a sample code algorithm of copying kernel (for GEMM without any data
converting operations like matrix transposition). Note that the work-group blocking
factors do not appear in the code; these appear in a routine interface code. We utilize
the developed auto-tuning system for tuning the copying kernels in the same manner
as for the GEMM kernel.
8

COPY KERNEL(m, k, Asrc , Adst )1: Get a two-dimensional global ID (y, x)
2: #pragma unroll (Mud )3: for i = y to y + Mwi − 1 do
4:
#pragma unroll (Nud )5:
for j = x to x + Nwi − 1 do
6:
Adst (i, j) := Asrc (i, j)
7:
end for
8: end for
Figure 6: Sample algorithm of copying kernel
An advantage of the present approach is that the necessary time for tuning copying
kernels are much smaller than the necessary time for GEMM kernel. A search spacefor each copying kernel is very small compared to a search space for the GEMM kernel.
The number of tunable parameters on a copying kernel is six; hence, to ﬁnd a nearoptimal parameter setting, we are required to measure the performance of hundredsof kernels only. On the GEMM kernel, the number of tunable parameters is about 20
and we are required to test tens of thousands of kernels, which takes us much time
for the auto-tuning. A disadvantage of the approach is that every routine requires
additional memory space as destination buﬀers of data copying. In case of a BLAS3
computation on square (n × n) matrices, additional memory space for temporarily
storing two matrices (2n2 data elements) is needed.

4

Experimental Results

We evaluate our implementation of BLAS3 routines on four diﬀerent GPUs (AMDRadeon R9 290X, FirePro W9100, HD 7970, and NVIDIA GeForce GTX Titan), a
many-core co-processor (Intel Xeon Phi 5110P), and a multi-core processor (Intel Xeon
Core i7 3960X). Speciﬁcations of the processors are shown in Table 1. This paper
describes the performance of each BLAS3 routine, and does not take into account data
communication time between the host and OpenCL device. Notable characteristics ofeach device are as follows:• The AMD Radeon R9 290X is one of the AMD Volcanic Islands series that adopts
the Graphic Core Next (GCN) architecture. The GCN is a RISC (Reduced Instruction Set Computer) SIMD (Single-Instruction Multiple-Data) architecture.
The R9 290X contains 44 compute units (CUs). All CUs share a 1 MB L2 cache.
Each CU has four vector units, a scalar unit, a 16 KB read/write L1 cache, anda 64 KB local data share (32 KB is available to a single wavefront). Each vector unit consists of 16 processing elements (PEs); thus, the R9 290X is equippedwith 2816 (= 44 · 4 · 16) PEs. The performance in double-precision is limited toone-eighth that in single precision.
9• The AMD FirePro W9100 is one of the AMD Volcanic Islands series for usage on
workstations. The W9100 and R9 290X have the same GPU codename (“Hawaii”)
and characteristics of W9100 are mostly same as R9 290X. An important diﬀerence
between the two GPUs is that the performance in double precision is half that insingle precision (no limitation).
• The AMD Radeon HD 7970 is one of the AMD Southern Islands series that also
adopts the GCN architecture. The HD 7970 contains 32 compute units (CUs).
All CUs share a 768 KB L2 cache. Each CU has four vector units, a scalar unit,
a 16 KB read/write L1 cache, and a 64 KB local data share (32 KB is available
to a single wavefront). Each vector unit consists of 16 processing elements (PEs);
thus, the HD 7970 is equipped with 2048 (= 32 · 4 · 16) PEs.
• The GeForce GTX Titan is a NVIDIA Kepler GF110 GPU. The GTX Titan
contains 14 Streaming Multiprocessor eXtremes (SMXs). All SMXs share a 1.5
MB L2 cache. Each SMX has a 64 KB L1 cache, a 64 KB shared memory, 192
CUDA cores; thus, the GTX Titan is equipped with 2688 (= 14cdot192). The
performance in double precision is one-third of the single precision.
• The Xeon Phi 5110P is one of the Intel’s ﬁrst commercial products in MIC (ManyIntegrated Core) architecture. The Xeon Phi features 60 in-order cores (59 coresavailable for computation) connected by a high-performance bidirectional ringinterconnect. Each core supports four-way hyper-threading to hide memory andinstruction latencies. The computation unit of the core is a vector unit which
has a 512-bit wide register ﬁle (32 registers per thread) and can execute 8-way
double-precision or 16-way single-precision ﬂoating-point SIMD instructions in asingle clock. Additionally, each core has a 32 KB L1 data cache, a 32 KB L1
instruction cache, and a 512 KB L2 cache.
• The Core i7 3960X is a six-core CPU of Intel’s Sandy Bridge architecture. The
CPU has 15 MB L3 cache. The processing core has a 32 KB L1 cache and 256
KB L2 cache. Also, the core has a 256-bit wide register ﬁle and can execute fourway double-precision or eight-way single-precision ﬂoating-point SIMD (AVX)
instructions in a single clock.

10

11

Southern
Islands
Tahiti
GK110

GeForce
GTX
Titan
NVIDIA
Kepler

1000
930
925
876
5632
5632
4096
5376
704
2816
1024
1792
5632
5238
3789
4709
704
2619
947
1570
4
16
3
6
320
320
264
288
32
32
32
48
Scratchpad Scratchpad
Scratchpad Scratchpad
AMD APP SDK 2.9
CUDA 5.5aba
14.4
13.352.1014
14.4
319.21ca: Catalyst driver versionb: FirePro Uniﬁed Driver versionc: CUDA driver version

Hawaii

CodenameCore clock speed [MHz]Maximum SP operations per clock
Maximum DP operations per clock
Peak SP performance [GFlop/s]
Peak DP performance [GFlop/s]
Global memory size [GBytes]
Global memory bandwidth [GByte/s]
Local memory size [KBytes]
Local memory type
OpenCL SDK
Display driver version
SP: Single-precision
DP: Double-precision

AMD
Volcanic Islands

VendorArchitecture

Device name

Table 1: OpenCL device speciﬁcations
RadeonFireProRadeonR9 290X
W9100
HD 7970

Sandy
Bridge-Core i7
3960X

Knights
Corner
1053
3300
1888
96
944
48
1988
316.8
994
158.4
8
320
32
32
GlobalGlobalIntel SDK 2013
-

MIC

Intel

Xeon Phi
5110P

4.1

Performance on Radeon R9 290X

This section presents performance measurement results on the Radeon R9 290X GPU.
We compare our implementation with clMathLibraries clBLAS 2.2.0. The clBLAS
is an open-source BLAS library in OpenCL and is being developed mainly by AMD
employees. In the clBLAS, a program for tuning BLAS routines is provided, and theperformance results of clBLAS are measured after applying the tuning program.
Fig. 7 shows the performance of SSYMM routine. Our implementation demonstrates higher performance for most of matrix sizes than clBLAS, and the performanceis 2.8 times higher on average. Relative performances of BLAS3 implementation to
clBLAS are shown in Fig. 8 (higher is better). The STRMM and SGEMM routines
achieves better performance for large matrix sizes. We do not see a large diﬀerence inperformances of the SSYRK and SSYR2K routines.
5000

Our work
clBLAS

Performance [GFlop/s]

4000

3000

2000

1000

0
0

1024

2048

3072 4096 5120
Matrix size [m=n]

6144

7168

8192

Figure 7: SSYMM performance on Radeon R9 290X
Fig. 9 shows a measured bandwidth of copying kernel for the SSYMM routine, and it
compares the bandwidth of auto-tuned copying kernel with that of non-tuned copying
kernel2 . The auto-tuning enables to develop a copying kernel whose performance ismuch higher than the non-tuned kernel. The bandwidth of tuned kernel is up to 110
GB/s while the bandwidth of non-tuned kernel is up to 2 GB/s. To show anotherevidence why tuning of a copying kernel is necessary, we present the ratio of copying
time to total computation time for the SSYMM routine in Fig. 10. When we do nottune the copying kernel, the copying time dominates over 65% of total computation
time even for matrix sizes of around 8000. By the auto-tuning, the ratio decreases
to less than 5% for these sizes. The auto-tuning of the copying kernel contributes todevelop high-performance routines.
2

All parameters (Mwg , Nwg , Mwi , Nwi , Mud , Nud ) of the non-tuned copying kernel are set as 1.

12

8

SSYMM
STRMM
SGEMM
SSYR2K
SSYRK

Relative performance

7
6
5
4
3
2
1
0
0

1024

2048

3072
4096
5120
Matrix size [m=n=k]

6144

7168

8192

Figure 8: Relative performance of our BLAS3 implementation to clBLAS
120

Tuned
Non-tuned

Bandwidth [GB/s]

100

80

60

40

20

0
0

1024

2048

3072 4096 5120
Matrix size [m=n]

6144

7168

8192

Figure 9: Bandwidth of copying kernel for SSYMM routine
Our SSYMM routine demonstrates a similar performance to our DGEMM routine
since we use the identical GEMM kernel usage in both routines. However, the otherthree routines (SSYRK, SSYR2K, and STRMM) requires small modiﬁcations in GEMM
kernel usage to implement as described in Section 3.3-3.5. Fig. 11 depicts the relative
performance of each BLAS3 routine to the performance of the SGEMM routine. The
performance of the three routines is lower than that of SSYMM. The SSYRK and
SSYR2K routines demonstrate an almost identical performance since the two routines

13

100

Ratio [%]

80

60
Non-tuned
Tuned

40

20

0
0

1024

2048

3072 4096 5120
Matrix size [m=n]

6144

7168

8192

Figure 10: Retio of copying time to total computation time for routine
share utilized OpenCL kernels. The performance of STRMM routine is the lowest
among all our BLAS3 routines on the GPU.
3

SSYMM
SSYR2K
SSYRK
STRMM

Relative performance

2.5

2

1.5

1

0.5

0
0

1024

2048

3072 4096 5120
Matrix size [m=n=k]

6144

7168

8192

Figure 11: Relative performance of each BLAS3 routine to SGEMM routine

4.2

Performance on Diﬀerent Processing Units

This section presents experimental results on diﬀerent processing units (see Table 1).
Table 2 describes the maximum measured performance of our implementation as well14

as a vendor BLAS library on the processing units. The R9 290X and W9100 show alittle diﬀerent performance’s behaviour though the GPU architecture is same (Volcanic
Islands). The STRMM routine works more eﬃciently than SSYRK and SSYR2K on
the W9100 while it is opposite on the R9 290X. Among the present devices, W9100 is
the fastest for the double-precision routines while R9 290X is the fastest for the singleprecision routines. For the GTX Titan, our implementation in OpenCL works slower
than cuBLAS. We suppose that the OpenCL runtime in CUDA toolkit is less tuned
than the CUDA runtime. Our implementation on the Xeon Phi and Core i7 is worse
than MKL. MKL is highly tuned in an assembly language level, and this result shows
another evidence that it is hard to outperform the program written in assembly level
using a high-level language.
Table 2: Maximum performance in GFlop/s on diﬀerent OpenCL devices
DeviceR9 290X
W9100
HD 7970
GTX Titan
Xeon Phi
Impl.
Ours#1 Ours#1 Ours#1 Ours#2 Ours#3
SGEMM
4650 3617 4223 3341 2913 2472 2300 3374
545 1701
SSYMM
4584 2013 4113 1698 2871 1047 2292 2620
544 1188
SSYRK
3853 3619 2295 3325 2096 2448 2183 3248
91 768SSYR2K
3822 3623 2223 3307 2774 2448 2211 3164
89 205
STRMM
3623 2649 3157 2463 2758 1779 2226 2908
13 510DGEMM
670 654 1778 1019
823 662958 1349
269 841DSYMM
665 488 1708 584
818 407955 1229
266 649
DSYRK
657 664 1138 1008
734 739929 1332
34 646DSYR2K655 665 1096 991
796 740910 1307
32 277
DTRMM
643 576 1576 732
789 538
820 1189
10 281#1: clBLAS (clMathLibraries clBLAS) 2.2.0#2: cuBLAS (CUDA BLAS in CUDA) 5.5 Toolkit#3: MKL (Math Kernel Library) 11.0 update 1

5Core i7
Ours #3
97 277
97 27395 25695 25491 26260 13659 13155 12555 11754 128Related Work

Igual et al. [6, 5] designed and evaluated a high-performance implementation of different Level-3 BLAS routines. Their implementation takes advantage of existing the
NVIDIA cuBLAS library and attains huge speedups over the cuBLAS. The cuBLAS
at that time does not show an equal performance for all Level-3 BLAS routines (theGEMM routine runs much more eﬃciently than the other routines). For their BLAS3
implementation, they presented matrix-panel (MP), panel-matrix (PM), panel-panel(PP) product based variants, where both matrix dimensions is large in the “matrix”
and one of the dimensions is small in the “panel”. In this point of view, our approach,presented in this paper, is considered to be a matrix-matrix (MM) product variant.

15

As for SYRK routine implementation on a GPU, Nath et al.[14] have used a threadblock reordering technique to limit the computation only to data blocks that are on the
diagonal or in the in the lower/upper triangular part of the matrix. All the threads
in diagonal blocks compute half of the block in a parallel fashion redundantly to avoid
conditional statements that have been required otherwise. Their SYRK implementation
showed higher performance than cuBLAS at that time [13]. Our approach for SYRK
routine also limits the computation in the lower triangular part and implements theredundant computation for the diagonal blocks.

6

Conclusion

In this paper, we have presented the implementation of diﬀerent Level-3 BLAS (BLAS3)
routines in OpenCL. A key in the implementation is to utilize data copying kernels to
realize each BLAS3 routine. In the present work, the copying kernels as well as GEMM
kernel are highly tuned by an auto-tuning technique. We have parameterized eachcopying kernel and implemented code generators which produce copying kernels forthe auto-tuning. Experimental results show that tuning the copying kernels greatly
contributes to developing of high-performance BLAS3 routines.
The proposed approach utilizes a highly-tuned GEMM kernel for other existing
BLAS3 routines. We found that the performance of the SYMM routine is almost
equivalent to that of GEMM routine while the performance of the SYRK, SYR2K, andTRMM routines is a little lower than that of GEMM. A future work is to investigate
the reason of such low performance for a better BLAS implementation.

References[1] AMD Inc.: AMD Accelerated Parallel Processing OpenCL Programming Guide,
rev2.4 (Dec 2012)[2] Bilmes, J., Asanovic, K., Chin, C.W., Demmel, J.: Optimizing matrix multiply
using PHiPAC : a portable, high-performance, ANSI C Coding methodology. Tech.
rep., Computer Science Department, University of Tennessee (May 1996), http://www.netlib.org/lapack/lawnspdf/lawn111.pdf[3] Du, P., Weber, R., Luszczek, P., Tomov, S., Peterson, G., Dongarra, J.: From
CUDA to OpenCL: Towards a performance-portable solution for multi-platform
GPU programming. Parallel Computing 38(8), 391–407 (Oct 2011)[4] Goto, K., van De Geijn, R.: High-performance implementation of the level-3 BLAS.
ACM Transactions on Mathematical Software (TOMS) 35(1), 4:1–4:14 (Jul 2008)[5] Igual, F.D.: Matrix computations on graphics processors and clusters of GPUs.
Ph.D. thesis, Universitat Jaume I. Castellón (2011)16

[6] Igual, F.D., Quintana-Qrtı́, G., van De Geijn, R.A.: Level-3 BLAS on a GPU:
Picking the low hanging fruit. Tech. Rep. DICC 2009-04-01, Depto. de Ingenieraı́yCiencia de Computadores, Universidad Jaime I de Castellón (2009)[7] Jiang, C., Snir, M.: Automatic tuning matrix multiplication performance on graphics hardware. pp. 185–194. IEEE (Sep 2005)[8] Khronos Group: OpenCL - The open standard for parallel programming of heterogeneous systems (Accessed Sep 18, 2014), http://www.khronos.org/opencl[9] Kågström, B., Ling, P., Van Loan, C.: GEMM-based level 3 BLAS: Highperformance model implementations and performance evaluation benchmark. ACM
Transactions on Mathematical Software 24(3), 268–302 (1998)[10] Kurzak, J., Tomov, S., Dongarra, J.: Autotuning GEMM kernels for the Fermi
GPU. IEEE Transactions on Parallel and Distributed Systems 23(11), 2045–2057(Nov 2012)[11] Matsumoto, K., Nakasato, N., Sedukhin, S.G.: Performance tuning of matrix multiplication in OpenCL on diﬀerent GPUs and CPUs. In: Proceedings of the 2012
SC Companion: High Performance Computing, Networking, Storage and Analysis(SCC 2012). pp. 396–405. IEEE CS’s Conference Publishing Service, Salt Lake
City, Utah, USA (Nov 2012)[12] Nath, R., Tomov, S., Dongarra, J.: An improved MAGMA GEMM for Fermigraphics processing units. International Journal of High Performance Computing
Applications 24(4), 511–515 (2010)[13] Nath, R., Tomov, S., Dongarra, J.: BLAS for GPUs. In: Bailey, D.H., Lucas, R.F.,
Williams, S. (eds.) Performance Tuning of Scientiﬁc Applications, chap. 4. CRC
Press (2010)[14] Nath, R., Tomov, S., Dongarra, J.: Accelerating GPU kernels for dense linear
algebra. In: Proceedings of the 9th International Meeting High Performance Computing for Computational Science (VECPAR 2010). LNCS, vol. 6449, pp. 83–92.
Springer (2011)[15] Whaley, R.C., Petitet, A., Dongarra, J.J.: Automated empirical optimizations ofsoftware and the ATLAS project. Parallel Computing 27(1-2), 3–35 (Jan 2001)17