Technical Report 2012-002

Implementing a Code Generator
for Fast Matrix Multiplication in OpenCL
on the GPU

Kazuya Matsumoto, Naohito Nakasato, Stanislav G. Sedukhin
July 2, 2012

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

Technical Report 2012-002
Title:
Implementing a Code Generator for Fast Matrix Multiplication in OpenCL on the GPU

Authors:
Kazuya Matsumoto, Naohito Nakasato, Stanislav G. Sedukhin
Key Words and Phrases:
matrix multiplication, OpenCL, GPU, auto-tuning

Abstract:
This paper presents results of an implementation of code generator for fast general matrix
multiply (GEMM) kernels. When a set of parameters is given, the code generator produces
the corresponding GEMM kernel written in OpenCL. The produced kernels are optimized
for high-performance implementation on GPUs from AMD. Access latencies to GPU global
memory is the main drawback for high performance. This study shows that storing matrix
data in a block-major layout increases the performance and stability of GEMM kernels. On the
Tahiti GPU (Radeon HD 7970), our DGEMM (double-precision GEMM) and SGEMM (singleprecision GEMM) kernels achieve the performance up to 848 GFlop/s (90% of the peak) and
2646 GFlop/s (70%), respectively.

Report Date:
7/2/2012

Written Language:
EnglishAny Other Identifying Information of this Report:
Manuscript has been accepted for Special Session: Auto-Tuning for Multicore and GPU (ATMG) and is toappear in the proceedings of the IEEE 6th International Symposium on Embedded Multicore SoCs (MCSoC12), Aizu-Wakamatsu City, Fukushima, Japan, September 20-22, 2012.

Distribution Statement:
First Issue: 10 copies

Supplementary Notes:

Distributed Parallel Processing Laboratory
The University of Aizu
Aizu-Wakamatsu
Fukushima 965-8580
JapanImplementing a Code Generator for Fast Matrix Multiplication
in OpenCL on the GPU
Kazuya Matsumoto, Naohito Nakasato, and Stanislav G. Sedukhin
Graduate School of Computer Science and Engineering
The University of Aizu
Aizu-Wakamatsu City, Fukushima, 965-8580 Japan
Email: {d8121101, nakasato, sedukhin}@u-aizu.ac.jp
Abstract—This paper presents results of an implementation of code generator for fast general matrix multiply(GEMM) kernels. When a set of parameters is given, thecode generator produces the corresponding GEMM kernel
written in OpenCL. The produced kernels are optimized
for high-performance implementation on GPUs from AMD.
Access latencies to GPU global memory is the main drawback
for high performance. This study shows that storing matrix
data in a block-major layout increases the performance andstability of GEMM kernels. On the Tahiti GPU (Radeon HD7970), our DGEMM (double-precision GEMM) and SGEMM(single-precision GEMM) kernels achieve the performance up
to 848 GFlop/s (90% of the peak) and 2646 GFlop/s (70%),
respectively.
Keywords-matrix multiplication; OpenCL; GPU; autotuning

main difference from the similar works [9], [13] is that
our code generator supports generating GEMM kernels
where matrix data are stored in memory not only in astandard row-/column-major layout, but also in a blockmajor layout.
The rest of this paper is organized as follows. Section
II describes our GEMM code generator. This section also
explains relations between parameters of the generator andthe OpenCL execution model. Section III shows results ofperformance evaluation on the AMD Tahiti GPU (Radeon
HD 7970). Finally, Section IV concludes the paper.
II. C ODE G ENERATORIn BLAS [1], GEMM is defined as
C ← αop(A)op(B) + βC,

I. I NTRODUCTION
Matrix-matrix multiply-add is a fundamental routine inlinear algebra, which is referred to as GEMM (GEneral
Matrix Multiply) in BLAS (Basic Linear Algebra Subroutines) standard [1]. GEMM appears in a lot of importantnumerical algorithms and is a building block of LAPACK(Linear Algebra PACKage) [2] and other Level-3 BLAS
routines [3]. The computational intensity and the regularity
of GEMM algorithms make them good candidates for
performance acceleration.
There are a lot of work on GEMM performance acceleration on CPUs [4]–[6] and GPUs [7]–[10]. Automatic performance tuning (or auto-tuning in short) is
an important technique to generate near-optimal GEMM
implementations. ATLAS (Automatically Tuned Linear
Algebra Software) [4] is a prominent auto-tuning software
library for BLAS routines running on CPUs. An autotuning system was also developed for GEMM kernels
in CUDA (Compute Unified Device Architecture) for
NVIDIA GPUs [7], [11], [12]. In addition, there is areport on auto-tuning methods for developing near-optimal
GEMM kernels in OpenCL (Open Compute Language) for
AMD GPUs [13].
Auto-tuning system needs two core components: code
generator and heuristic search engine. Code generator
produces codes parameterized from a pre-defined code
template. Heuristic search engine runs the generated codesand searches the best set of parameters for increasing
performance using a feedback loop. This study proposes acode generator for fast GEMM kernels written in OpenCL
as a preliminary work for an auto-tuning system. The

where both α and β are scalar values, and op(A), op(B)and C are m × k, k × n and m × n matrices, respectively.
The op(X) takes X (non-transposed matrix) or X T (transposed matrix); thus, there are four multiplication types:(a) C ← αAB + βC,(b) C ← αAB T + βC,(c) C ← αAT B + βC,(d) C ← αAT B T + βC,
in each real (double or single) precision. In the following
description, if not specified, we assume that matrix data
are stored in a row-major layout, which is standard in C
Programming Language.
Our code generator takes a set of parameters as the
input. When the input is given, the code generator produces the corresponding GEMM kernel code written in
OpenCL as the output. We can set different parameters
to the generator for each GEMM type (14 parametersfor (c) type and 12 parameters for the other (a), (b),(d) types). Six parameters define the blocking factors andthe other parameters are related to optimization of a way
for accessing matrix data. The followings describe themeaning of every parameter.
A. Blocking
A straightforward implementation of GEMM is a threenested-loop program. Blocking a GEMM algorithm isnecessary for getting high performance on a processor witha multi-level memory hierarchy. It is because the blocking
increases the data reuse ratio which is required to toleratememory access latencies.

We can consider that the GPU has three levels of memory spaces: off-chip memory (global memory), on-chip
memory (cache or local memory), and private memory(register file). Roughly speaking, the off-chip memoryis the largest memory space with the lowest bandwidth,
the private memory is the smallest memory space withthe highest bandwidth, and the on-chip memory is in the
middle of them.
OpenCL is an open standard for general-purpose parallel programming on heterogeneous platforms [14], [15].
The OpenCL execution model covers the three memory
spaces. OpenCL specifies a C99-based language that allows to write parallel functions called kernels. When akernel is submitted for execution on the GPU, an index
space, which is called NDRange, is defined. Each instancein NDRange is named work-item, and several work-items
are organized into a work-group. Every work-item canaccess its own private memory which is not visible fromother work-items. Work-items in a single work-groupconcurrently run on the processing elements of a computeunit, and share their own on-chip memory. Global memory
can be accessed from all work-items, though there is noway for synchronization of work-items in different workgroups during a kernel execution.
We use a two-level (larger and smaller) blocking in
our GEMM kernel template since it is effective for highperformance implementation to utilize the three levelsof memory and execution model. Let ml , nl , kl be thelarger blocking factors. The blocking divides the threematrices A, B, C into (m/ml ) × (k/kl ), (k/kl ) × (n/nl )and (m/ml )×(n/nl ) grids of ml ×kl , kl ×nl and ml ×nl
blocks, respectively. Fig. 1 depicts a matrix multiplyadd partitioned with the blocking factors. Workloads fora single ml × nl block of C are allocated to a workgroup. The work-group is associated with multiplicationof a ml × k stripe of A with a k × nl stripe of B andaddition of the product to an ml × nl block of C for the
final result.
An algorithm for the stripe-by-stripe multiplication iterates k/kl times in the outermost loop of our GEMM
algorithm (let us consider matrix sizes m, n, k as multiples of ml , nl , kl , respectively). In every iteration of
the outermost loop, the work-group produces the partial
result of ml × nl block by multiplying ml × kl blockon kl × nl block and adding the product to the ml × nl
block. Fig. 2 shows the block-by-block multiplication.
Each block is further divided with the smaller blocking
factor (ms , ns , ks ) such that a work-item in the workgroup is in charge of multiplication of ms × ks sub-blockof A by ks × ns sub-block of B and accumulation of the
product on an ms × ns sub-block of C.
B. Way for accessing matrix data
For storing matrix data in OpenCL, we can choose
either of two types of memory objects: buffer objects andimage objects. Buffer objects are like one-dimensional
arrays in C Language and buffer data are sequentiallystored in global memory. Image objects are intended to

n

n

k

nl

kl
ml

m +=

C

kl
ml

nl

m ×

A

k

B

Figure 1. Matrix-matrix multiply-add with the larger blocking factor(ml , nl , kl ) in an OpenCL kernel on a GPU

nl
ms

kl

ns

ks
ml

C

nl
ks

ms
ml

+=

× ns
kl

A
B

Figure 2. Block-block multiply-add with the smaller blocking factor(ms , ns , ks ) in a work-group on a compute unit

contain pixel data for image processing, and image data
are stored in special global memory called texture memory.
Using images may take advantage of spatial memory
locality for increasing performance. Our code generator
has a parameter to designate using either buffer or image
for reading data from matrices A and B. For data related
to matrix C, buffer is always used.
As another parameter, we can designate a width ofvector variables to the generator. Kernels using different
widths show different performance. Note that, when akernel uses image objects, the width is limited to two in
double-precision (double2) and four in single-precision(float4).
GPUs have a local memory1 to share data between
work-items in a single work-group. The bandwidth of thelocal memory is higher than that of the L1 cache. For
instance, in case of the Tahiti GPU (Radeon HD 7970),
the peak bandwidth of the local memory is 3789 GB/s
while that of the L1 cache is 1894 GB/s [15]. A drawbackof using the local memory is that it requires a barrier
synchronization between the work-items, which costs a
certain amount of time. The drawback leads to the factthat using the local memory is not always good for higher
performance; hence, we have added to the code generator
a parameter for either using the local memory for eachmatrix of A and B or not.
It is known that AT B + C kernel is the fastest among
all four GEMM types in row-major on the Cypress GPU,because of its efficient memory access patterns [10], [16].
For example, the maximum performance of AT B + C
double-precision kernel is 472 GFlop/s while that of
AB + C kernel is 359 GFlop/s on the Cypress GPU(Radeon HD 5870) [10]. AT B + C kernel is the fastest
also on the Tahiti GPU (experimental results are shownin Section III). This means that using the AT B + C
1 In GPUs from AMD, the local memory is named local data store(LDS).

Table ITAHITI GPU SPECIFICATIONS
Codename
Product nameCore clock speed [MHz]Number of compute units (CUs)
Number of double-precision (DP) units
Number of single-precision (SP) units
Peak DP performance [GFlop/s]
Peak SP performance [GFlop/s]
Memory clock speed [MHz]Global memory size / GPU [GB]L2 cache size / GPU [kB]
L1 cache size / CU [kB]
Local memory size / CU [kB]Global memory peak bandwidth [GB/s]
L2 cache peak bandwidth [GB/s]
L1 cache peak bandwidth [GB/s]
Local memory peak bandwidth [GB/s]

Tahiti
Radeon HD 7970
925
32
512
2048
947
3789
1375
3
768
16
64
264
710
1894
3789

kernel also for the other three GEMM types will contribute
performance improvement. A solution for the kernel utilization is to copy matrix data in a transposed form
if necessary. Such optimization was used previously by
Du et al. [9]. They showed that the optimization works
effectively for large matrices in which the O(n2 ) time of
the copying is negligible compared with O(n3 ) time of
execution of GEMM kernel.
We extend the solution with matrix transposition suchthat matrix data are stored in global memory in a blockmajor layout to increase a spatial locality. Currently, as an
option, the code generator supports producing AT B + C
kernels where matrix data are supposed to be aligned ineither of two block-major layouts shown in Fig. 3 (on an
m×k transposed matrix A with the blocking factor ml , kl ).
Fig. 3(a) shows a column-block major layout and data ineach k × ml column-block are sequentially stored in a
row-major order. We name this layout column block layout(CBL). In CBL, matrix data needed for a multiplication ofk × ml stripe and k × nl stripe are aligned in contiguous
memory locations. Fig. 3(b) depicts a row-block major
layout and data in each kl × ml sub-block of a kl × m
row-block are stored in a row-major order. Let us namethis layout row block layout (RBL). In RBL, matrix datarequired for a multiplication of kl × ml block and kl × nl
block are aligned in memory. There are past studies on thebenefits of using block major layouts to improve memory
utilization in CPUs [17], [18].
III. P ERFORMANCE E VALUATION
We evaluate the performance of DGEMM (doubleprecision GEMM) and SGEMM (single-precision GEMM)
kernels generated by our code generator on the AMD
Tahiti GPU (Radeon HD 7970). The GPU specifications
are shown in Table I. The GPU system runs on Ubuntu
10.04. The installed display driver is AMD Catalyst 12.3.
We use AMD Accelerated Parallel Processing (APP)SDK v2.6 for OpenCL software development. The present
performance evaluation does not take into account datatransfer time between host (CPU) and GPU.
The implemented code generator has been used tosearch the best set of parameters which produces thefastest GEMM kernel. We tested at least ten thousand
kernel variants per GEMM type. Such large number ofvariants were heuristically chosen. Our heuristic searchengine measures the performance of generated kernels inorder according to a feedback of the measured performance. Note that the search engine is not matured. It takes
around three hours to find the best set of parameters foreach GEMM type and takes 24 hours to find all the fastest
DGEMM and SGEMM kernels shown in this paper. The
following shows the procedure for selecting the fastest
kernel:1) Measuring the performance in GFlop/s of every
generated GEMM kernel for two problem sizes
n = 1536 and 4096.
2) Further measuring the performance of the fastest 50
kernels for problems sizes n (256 ≤ n ≤ 8192 in
multiples of 256) among a large number of kernels
tested in 1).
3) Selecting the fastest kernel among the 50 kernels
tested in 2).
The performance of the fastest DGEMM and SGEMM
kernels is shown in Fig. 4. Table II shows sets of parameters and the observed maximum performance of the kernels. In row-major, the maximum DGEMM performance(790 GFlop/s) of AT B + C kernel is higher than that(689 GFlop/s) of AB + C kernel, though the AT B + C
performance fluctuates depending on matrix sizes probably
due to memory bank conflicts. It was found that there
exist several AT B + C kernels demonstrating a stable(non-fluctuated) performance tendency with around 620
GFlop/s. The AT B + C kernel in CBL is most superior
to the other shown kernels in terms of the performanceand the stability; the achieved performance for this kernel
is up to 848 GFlop/s (90% of the peak performance)
in double-precision and 2646 GFlop/s (70%) in singleprecision. The present code generator limits the size of all
blocking factors to a power of two. Therefore, there may
be some room for improving the performance.
To make use of the AT B + C kernel in CBL, matrix
data have to be copied into CBL before executing thekernel. We need to prepare two kinds of copying kernels: akernel copying data into CBL without matrix transposition
and a kernel copying data into CBL with transposition.
For instance, in case of AB + C, the matrix A has tobe copied with transposition and the matrix B has to becopied without transposition. Every copying kernel for asquare matrix-matrix multiplication (i.e., m = n = k)reads n×n matrix from a memory space in global memory
and then, writes the n×n matrix to another memory spacein global memory. The memory bandwidth of every n × nmatrix-to-matrix copying kernel can be calculated as
BW (n) [Bytes / sec] =2n2 · S [Bytes],
Time for copying [sec]where S = 8 for the double precision and S = 4 for the
single precision. The measured memory bandwidth of thecopying kernels is shown in Fig. 5. The copying kernels

ml

ml

kl

k

k

m

m

(a) Column block layout (CBL)(b) Row block layout (RBL)

Figure 3.

Block-major layouts on an m × k transposed matrix with the blocking factor ml , kl
Table IIPARAMETERS AND PERFORMANCE OF THE FASTEST GEMM KERNELS

Kernel type
AT B + C in CBL

ml , nl , kl

ms , ns , ks

VectoraSharedb

ImagecPerformance [Gflop/s]

64,16,16
64,16,16
256,8,16
32,128,128
128,128,256
128,32,8
128,64,4
128,64,256

4,4,2
4,4,2
4,4,2
4,4,4
16,8,4
8,8,8
8,4,4
8,4,4

2
2
2
4
4
2
4
4

B
B
A
B
-

A

848
812
790
689
2646
2577
2488
2382

AT B + C in RBL
AT B + C in row-major
AB + C in row-major
AT B + C in CBL
AT B + C in RBL
SGEMM
AT B + C in row-major
AB + C in row-major
a. Width of vector variables
b. Matrix whose data are shared in local memoryc. Matrix whose data are loaded from image
DGEMM

900
2500

800

Performance [GFlop/s]

Performance [GFlop/s]

700
600
500
400
300
T

200
100

A B + C in CBLT
A B + C in RBLT
A B + C in row-major
AB + C in row-major

2000

1500

1000
ATB + C in CBL
ATB + C in RBLT
A B + C in row-major
AB + C in row-major

500

0

0
1024

2048

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

6144

7168

8192

DGEMM

1024

2048

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

6144

7168

8192

SGEMM

Figure 4. Performance of the fastest GEMM kernels produced by our code generator on the Radeon HD 7970. Matrix sizes are in multiples of 256.

are not deeply tuned; hence, the measured bandwidthof 20-160 GB/s is not high compared with the peakbandwidth (264 GB/s).
Fig. 6 shows the performance of our GEMM implementation with the fastest AT B + C kernel in CBL supportedby the coying kernels. In the figure, the performance iscompared with that of the AB + C kernel in row-major
and GEMM routines from AMD Accelerated Parallel Processing Math Libraries (APPML) clBLAS 1.8.269 [19].
The time for the copying is amortized when matrix sizes
become larger, and our implementation shows higher
performance than the AB + C kernel in row-major whenn ≥ 1536 in DGEMM and n ≥ 4096 in SGEMM. Also,
the performance of our implementation is higher than thatof APPML. The maximum performance of our DGEMM
and SGEMM AB + C implementations is 823 GFlop/s(87% of the peak) and 2541 GFlop/s (67%), respectively.
The performance difference among all four GEMM types
is not large and within 3% for DGEMM and 5% for
SGEMM.
When a matrix size is not in multiples of a blocking
factor, we use a zero padding technique. In the padding
technique, the values of padded portion in matrices A, B
are initialized as zero such that our GEMM implementa-

160

140

140

120

120

Bandwidth [GB/s]

Bandwidth [GB/s]

160

100
80
60
40

100
80
60
40

Copying A w/o transposition
Copying B w/o transposition
Copying A with transposition
Copying B with transposition

20
0
1024

2048

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

Copying A w/o transposition
Copying A with transposition

0
6144

7168

8192

1024

Double-precision

2048

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

6144

7168

8192

Single-precision

Measured memory bandwidth of copying kernels for utilizing the AT B +C kernel in CBL. In single-precision case, since the corresponding

Figure 5.
bandwidth for A and B is identical, the bandwidth for A is only shown.
900

2500

800

Performance [GFlop/s]

Performance [GFlop/s]

700
600
500
400
300

ATB + C with ATB+C kernel in CBL
AB + C with ATB+C kernel in CBL
AB + C kernel in row-major
ATB + C of APPML
AB + C of APPML

200
100

2000

1500

1000

ATB + C with ATB+C kernel in CBL
AB + C with ATB+C kernel in CBL
AB + C kernel in row-major
ATB + C of APPML
AB + C of APPML

500

0

0
1024

2048

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

6144

7168

8192

DGEMM

1024

2048

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

6144

7168

8192

SGEMM

Figure 6. Performance of our GEMM implementation using the fastest AT B + C kernel in CBL supported by the copying kernels, and theperformance comparison with the AB + C kernel in row-major and GEMM routines of AMD APPML clBLAS 1.8.269 [19].

tion can utilize the fastest GEMM kernel in CBL, which
is optimized for matrix sizes in multiples of the blocking
factor. The performance of the GEMM implementation
with the padding is shown in Fig. 7. One disadvantageof the padding is that it causes a performance degradation due to the extra computation for padded portions.
The performance degradation is small for relatively large
problem sizes. Another side effect of the padding is thatit requires three extra buffers for matrices, and as a result,the DGEMM computation becomes impossible in the caseof m = n = k > 6800 on the 3 GB memory.

In the SGEMM case of Fig. 7, two performance results
are shown. The green line indicates the performance of
the implementation with the already mentioned AT B + C
kernel in CBL which uses 256 as one of the blocking
factors. In the padding, the larger a blocking factor is,the bigger the performance deterioration becomes. To see
effects of using different blocking factors on the performance, we test another SGEMM implementation with an
AT B +C kernel in CBL having 64 as the largest blocking

factor2 . The gray line of Fig. 7 represents the performance.
As shown, using the smaller blocking factor refrains adrastic performance deterioration although the maximum
performance (2430 GFlop/s) of the implementation is alittle lower.
IV. C ONCLUSION
This paper has presented our code generator for searching fast GEMM kernels. The performance of GEMM
kernels running on GPUs is determined by many factors
and it is extremely difficult and time-consuming to develop
highly optimized kernels by hand tuning. We think thatsearching a number of kernel pattens with an auto-tuning
system is a good solution for quick development of fast
GEMM kernels. Currently, the heuristic search engine in
our implemented system is not matured and our system
is not a perfect auto-tuning system; however, results ofthis study show that searching exhaustive kernel patterns
by the code generator still contributes development of fast
2 Parameters of the SGEMM kernel: {m , n , k } = {64, 32, 32},ll l
{ms , ns , ks } = {4, 8, 4}, the width of vector variables is 4, andcommon data of the matrix B is shared in local memory.

900
2500

800

Performance [GFlop/s]

Performance [GFlop/s]

700
600
500
400
300
200

2000

1500

1000

500

100

{128,128,256}{64,32,32}

0

0
1024

2048

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

6144

7168

8192

DGEMM

1024

2048

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

6144

7168

8192

SGEMM

Figure 7. Performance of our AB + C GEMM implementation with a padding technique. In SGEMM, two performance results are shown in
different blocking factors {ml , nl , kl }.

GEMM kernels with modest development time. The code
generator supports production of AT B+C GEMM kernels
where matrix data are sequentially stored in a block-major
layout. It was shown that using the best kernel in theblock-major layout results in achieving higher maximum
performance and performance stability than using a kernel
in a row-major layout on the Tahiti GPU (Radeon HD7970).
Since programs written in OpenCL are portable across
OpenCL supported devices including multi-core CPUs andGPUs from different vendors, a possible future work is toevaluate GEMM kernels produced by the code generator
on several different architectures of CPUs and GPUs.
Implementing a code generator for other BLAS routines
is another future work.
R EFERENCES[1] Basic Linear Algebra Subprograms Technical Forum,“Basic Linear Algebra Subprograms Technical (BLAST)
Forum Standard,” Aug. 2001. [Online]. Available: http://www.netlib.org/blas/blast-forum/blas-report.pdf[2] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel,
J. Dongarra, J. D. Croz, A. Greenbaum, S. Hammarling,
A. McKenney, and D. Sorensen, LAPACK User’s Guide,
3rd ed. Philadelphia, PA, USA: Society for Industrial andApplied Mathematics, 1999.
[3] B. Kågström, P. Ling, and C. van Loan, “GEMM-based
level 3 BLAS: High-performance model implementations
and performance evaluation benchmark,” ACM Transactions on Mathematical Software, vol. 24, no. 3, pp. 268–
302, 1998.
[4] R. C. Whaley, A. Petitet, and J. J. Dongarra, “Automated empirical optimizations of software and the ATLAS
project,” Parallel Computing, vol. 27, no. 1-2, pp. 3–35,
Jan. 2001.
[5] K. Goto and R. A. van de Geijn, “Anatomy of highperformance matrix multiplication,” ACM Transactions onMathematical Software (TOMS), vol. 34, no. 3, pp. 1–25,
May 2008.

[6] J. Kurzak, W. Alvaro, and J. Dongarra, “Optimizing matrix
multiplication for a short-vector SIMD architecture CELL
processor,” Parallel Computing, vol. 35, no. 3, pp. 138–
150, Mar. 2009.
[7] J. Kurzak, S. Tomov, and J. Dongarra, “Autotuning
GEMM kernels for the Fermi GPU,” IEEE Transactions on
Parallel and Distributed Systems, 2012 (in press). [Online].
Available: http://dx.doi.org/10.1109/TPDS.2011.311[8] G. Tan, L. Li, S. Triechle, E. Phillips, Y. Bao, and N. Sun,
“Fast implementation of DGEMM on Fermi GPU,” inProceedings of 2011 International Conference for High
Performance Computing, Networking, Storage and Analysis(SC ’11). Seattle, WA, USA: ACM, Nov. 2011.
[9] P. Du, R. Weber, P. Luszczek, S. Tomov, G. Peterson,
and J. Dongarra, “From CUDA to OpenCL: Towards a
performance-portable solution for multi-platform GPU programming,” Parallel Computing, vol. 38, no. 8, pp. 391–
407, Oct. 2011.
[10] N. Nakasato, “A fast GEMM implementation on the Cypress GPU,” ACM SIGMETRICS Performance Evaluation
Review, vol. 38, no. 4, pp. 50–55, Mar. 2011.
[11] J. Kurzak, P. Luszczek, S. Tomov, and J. Dongarra,“Preliminary results of autotuning GEMM kernels for
the NVIDIA Kepler architecture - GeForce GTX 680,”
2012, LAPACK Working Note 267. [Online]. Available:
http://www.netlib.org/lapack/lawnspdf/lawn267.pdf[12] Y. Li, J. Dongarra, and S. Tomov, “A note on auto-tuning
GEMM for GPUs,” in Proceedings of the 9th International
Conference on Computational Science (ICCS ’09), ser.
LNCS, vol. 5544. Baton Rouge, LA: Springer-Verlag,
May 2009, pp. 884–892.
[13] C. Jang. (Accessed Jun. 27, 2012) GATLAS GPU
Automatically Tuned Linear Algebra Software. [Online].
Available: http://golem5.org/gatlas[14] Khronos Group. (Accessed Jun. 27, 2012) OpenCL- The open standard for parallel programming ofheterogeneous systems. [Online]. Available: http://www.
khronos.org/opencl[15] AMD Inc., “AMD Accelerated Parallel Processing OpenCL
Programming Guide, rev2.2,” Jun. 2012.

[16] K. Matsumoto, N. Nakasato, T. Sakai, H. Yahagi, and S. G.
Sedukhin, “Multi-level optimization of matrix multiplication for GPU-equipped systems,” in Procedia Computer
Science, vol. 4. Elsevier B.V., Jun. 2011, pp. 342–351.
[17] V. K. Prasanna, N. Park, and B. Hong, “Tiling, block
data layout, and memory hierarchy performance,” IEEE
Transactions on Parallel and Distributed Systems, vol. 14,
no. 7, pp. 640–654, Jul. 2003.

[18] F. G. Gustavson, “New generalized data structures formatrices lead to a variety of high performance algorithms,”
in Proceedings of the 4th International Conference on
Parallel Processing and Applied Mathematics (PPAM ’01),
ser. LNCS, vol. 2328, 2002, pp. 418–436.
[19] AMD Inc. (Accessed Jun. 27, 2012) AMD Accelerated
Parallel Processing Math Libraries (APPML). [Online].
Available: http://developer.amd.com/libraries/appmathlibs