Analytical Modeling is Enough for High Performance BLIS

A Analytical Modeling is Enough for High Performance BLIS TZE MENG LOW, The University of Texas at Austin FRANCISCO D. IGUAL, Universidad Complutense de Madrid TYLER M. SMITH, The University of Texas at Austin ENRIQUE S. QUINTANA-ORTI, Universidad Jaume I We show how the BLAS-like Library Instantiation Software (BLIS) framework, which provides a more detailed layering of the GotoBLAS (now maintained as OpenBLAS) implementation, allows one to analytically determine optimal tuning parameters for high-end instantiations of the matrix-matrix multiplication. This is of both practical and scientific importance, as it greatly reduces the development effort required for the implementation of the level-3 BLAS while also advancing our understanding of how hierarchically layered memories interact with high performance software. This allows the community to move on from valuable engineering solutions (empirically autotuning) to scientific understanding (analytical insight). Categories and Subject Descriptors: G.4 [Mathematical Software]: Efficiency General Terms: Algorithms, Performance Additional Key Words and Phrases: linear algebra, libraries, high-performance, matrix multiplication, analytical modeling ACM Reference Format: ACM V, N, Article A (January YYYY), 19 pages. DOI:http://dx.doi.org/10.1145/0000000.0000000 1. INTRODUCTION The field of dense linear algebra (DLA) was among the first to realize, understand, and promote that standardizing (de facto) an interface to fundamental primitives supports portable high performance. For almost four decades, those primitives have been the Basic Linear Algebra Subprograms (BLAS) [Lawson et al. 1979; Dongarra et al. 1988; Dongarra et al. 1990]. The expectation was and still remains that some expert will provide to the community a high performance implementation of the BLAS every time a new architecture reaches the market. For this purpose, many hardware vendors currently enroll a numerical libraries group and distribute well-maintained libraries (e.g., Intel’s MKL [Intel 2015], AMD’s ACML [AMD 2015], and IBM’s ESSL [IBM 2015] libraries), while over the years we have enjoyed a number of alternative Open Source solutions (e.g., GotoBLAS [Goto and van de Geijn 2008b; 2008a], OpenBLAS [OpenBLAS 2015], ATLAS [Whaley and Dongarra 1998]). Nevertheless, these solutions all require considerable effort in time and labor for each supported target architecture. In order to reduce the exertion of developing high performance implementations, ATLAS and its predecessor PHiPAC [Bilmes et al. 1997b] introduced autotuning to the field of high performance computing. The fundamental rationale behind these two libraries is that optimizing is too difficult and time-consuming for a human expert and that autogeneration in combination with autotuning is the answer. Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies show this notice on the first page or initial screen of a display along with the full citation. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, to republish, to post on servers, to redistribute to lists, or to use any component of this work in other works requires prior specific permission and/or a fee. Permissions may be requested from Publications Dept., ACM, Inc., 2 Penn Plaza, Suite 701, New York, NY 10121-0701 USA, fax +1 (212) 869-0481, or p[email protected]. c YYYY ACM 0000-0000/YYYY/01-ARTA $15.00 DOI:http://dx.doi.org/10.1145/0000000.0000000 ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY. A:2 Tze M. Low et al. The work in [Yotov et al. 2005] shows that the ATLAS approach to optimizing the general matrix-matrix product (gemm) can be, by and large, analytically modelled. Thus, it reveals that autotuning is unnecessary for the operation that has been tauted by the autotuning community as the example of the success of autotuning. The problem with that work ([Yotov et al. 2005]) is that the ATLAS approach to optimizing gemm had been previously shown to be suboptimal, first in theory [Gunnels et al. 2001] and then in practice [Goto and van de Geijn 2008b]. Furthermore, ATLAS leverages an inner kernel optimized by a human expert, which still involves a substantial manual encoding. Precisely, the presence of this blackbox kernel makes analytical modeling problematic for ATLAS, since some of the tuning parameters are hidden inside. GotoBLAS (now supported as OpenBLAS) also builds on a substantial inner kernel that is implemented by a human, turning the analysis equally difficult without a complete understanding of the complex details embedded into the inner kernel. (Until recently, no paper that detailed the intricacies of the GotoBLAS inner kernel existed.) BLIS (BLAS-like Library Instantiation Software) [Van Zee and van de Geijn 2014] is a new framework for instantiating the BLAS. A key benefit of BLIS is that it is a productivity multiplier, as the framework allows developers to rapidly unfold new high-performance implementations of BLAS and BLAS-like operations on current and future architectures [Van Zee et al. 2014]. From the implementation point of view, BLIS modularizes the approach underlying GotoBLAS2 (now adopted also by many other implementations) to isolate a microkernel that performs gemm upon which all the level-3 BLAS functionality can be built. Thus, for a given architecture, the programmer only needs to develop an efficient microkernel for BLIS, and adjust a few key parameter values1that optimize the loops around the micro-kernel, to automatically instantiate all the level-3 BLAS. Importantly, BLIS features a layering that naturally exposes the parameters that need to be tuned. While BLIS simplifies the implementation of BLAS(-like) operations, a critical step to optimize BLIS for a target architecture is for the developer to identify the specific parameter values for both the micro-kernel and the loops around it. Conventional wisdom would dictate that empirical search must be applied. In this paper, we adopt an alternative model-driven approach to analytically identify the optimal parameter values for BLIS. In particular, we apply an analysis of the algorithm and architecture, similar to that performed in [Yotov et al. 2005], to determine the data usage pattern of the algorithm ingrained within the BLIS framework, and to build an analytical model based on hardware features in modern processor architectures. The specific contributions of this work include: —An analysis of the best known algorithm. We address the algorithm underlying BLIS (and therefore GotoBLAS and OpenBLAS) instead of the algorithm in ATLAS. This is relevant because it has been shown that, on most architectures, a hand-coded BLIS/GotoBLAS/OpenBLAS implementation (almost) always yields higher performance than an implementation automatically generated via ATLAS [Van Zee et al. 2014], even if ATLAS starts with a hand-coded inner kernel. —A more modern model. We accommodate a processor model that is more representative of modern architectures than that adopted in [Yotov et al. 2005]. Concretely, our model includes hardware features such as a vector register file and a memory hierarchy with multiple levels of set-associative data caches. This considerably improves upon the analytical model in [Yotov et al. 2005] which considered one level only of fully associative cache and ignored vector instructions. 1The micro-kernel itself is characterized by three additional parameters, which the programmer has to consider when implementing it. ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY. Analytical Modeling is Enough for High Performance BLIS A:3 —A more comprehensive model. Our analytical model is more comprehensive in that it includes the parameter values that also characterize the micro-kernel. These are values that a developer would otherwise have to determine empirically. This is important because the best implementations provided by ATLAS often involve loops around a hand-coded (black-box) kernel. Since the internals of the hand-coded kernel are not known, the model in [Yotov et al. 2005] was not able to identify the globally optimal parameter values. —Competitive with expert-tuned implementations. Unlike previous work comparing the performance attained against auto-generated implementations in C obtained by ATLAS [Yotov et al. 2005; Kelefouras et al. 2014], which are typically not competitive with those that use hand-coded inner kernels (so-called “ATLAS unleashed” in [Yotov et al. 2005]), this paper shows that analytically-obtained parameter values can yield performance that is competitive with manually-tuned implementations that achieve among best-in-class performance. Hence, this paper provides an answer to the question posed in [Yotov et al. 2005] —i.e. whether empirical search is really necessary in this context—, by demonstrating that analytical modeling suffices for high performance BLIS implementations which then shows that careful layering combined with analytical modeling is competitive with GotoBLAS, OpenBLAS, and vendor BLAS (since other papers show BLIS to be competitive with all these implementations). In this paper, we restrict ourselves to the case of a single-threaded implementation of gemm in double precision. How to determine optimal parameter values for a multithreaded implementation is orthogonal to this paper, and is at least partially answered in [Smith et al. 2014]. 2. APPROACHES TO IDENTIFY THE OPTIMAL PARAMETER VALUES FOR GEMM It has been argued that empirical search is the only way to obtain highly optimized implementations for DLA operations [Demmel et al. 2005; Bilmes et al. 1997a; Whaley and Dongarra 1998], and an increasing number of recent projects (Build-To-Order BLAS [Belter et al. 2010] and AuGEM [Wang et al. 2013]) now adopt empirical search to identify optimal parameter values for DLA algorithms. The problem with empirical-based approaches is that they unleash a walloping search space, due to the combination of a large number of possible values for a substantial set of parameters. Therefore, even with the help of heuristics, generating, executing, and timing all variants of the program, with their unique combinations of parameter values, requires a considerable amount of time (sometimes measured in days). In addition, empirical search has to be performed on the target machine architecture. This implies that a high performance implementation of BLAS for a new machine architecture cannot be developed until the programmer is granted access to the target machine, and then it still requires a significant amount of time. This is especially critical on state-of-the-art embedded architectures, where cross-compiling processes are necessary, greatly complicating or even invalidating in practice the application of empirical-based approaches. A refinement of the empirical approach is iterative optimization [Knijnenburg et al. 2002; Kisuki et al. 2000; 2000]. As the name suggests, instead of a single run that explores the entire search space, multiple optimization passes are executed. There exist different techniques to perform iterative optimization, but they all share the following similarities: An initial sampling run is performed while a monitoring tool captures specific performance information such as the number of cache misses. Using the captured information and the application of some heuristics, the search space is trimmed/refined, and new program variants (with new sets of parameter values) are generated for further exploration in subsequent passes. Although iterative optimization represents an improvement over empirical search, the information obtained, e.g. from the performance counters, may not be sufficient to limit the ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY. A:4 Tze M. Low et al. 500000 1e+06 1.5e+06 2e+06 2.5e+06 3e+06 3.5e+06 4e+06 4.5e+06 kc mc TLB misses 32 64 96 128 160 192 224 256 288 320 352 384 416 448 480 512 544 576 608 640 64 128 192 256 320 384 448 512 576 640 2e+06 4e+06 6e+06 8e+06 1e+07 kc mc L2 misses 32 64 96 128 160 192 224 256 288 320 352 384 416 448 480 512 544 576 608 640 64 128 192 256 320 384 448 512 576 640 1e+07 1.5e+07 2e+07 2.5e+07 3e+07 3.5e+07 4e+07 kc mc L1 misses 32 64 96 128 160 192 224 256 288 320 352 384 416 448 480 512 544 576 608 640 64 128 192 256 320 384 448 512 576 640 8 10 12 14 16 18 20 kc mc GFLOPS 32 64 96 128 160 192 224 256 288 320 352 384 416 448 480 512 544 576 608 640 64 128 192 256 320 384 448 512 576 640 Fig. 1: Performance metrics when empirically testing different parameter values. Clockwise from top-left: GFLOPS, and misses in the L1, TLB and L2 caches. A heuristic for reducing the search space is to find the parameter values that minimize the cache misses and TLB misses. From the above data, this still leaves a relatively large space that needs to be searched. search space. Consider the graphs in Figure 1, which report the GFLOPS (i.e., billions of floating-point arithmetic operations, or flops, per second), and L1, L2 and TLB cache misses observed for the execution of a tuned implementation of dgemm from BLIS on a single core of an Intel Xeon E3-1220 processor (3.1 GHz), while varying two of the optimization parameters only (mcand kc, to be introduced later). The top-left plot illustrates the challenge for empirical search: generate and evaluate all different combinations of an exponentially growing search space. (BLIS, e.g., needs to optimize 5–6 different parameters, which is a reduced quantity when compared with other solutions including ATLAS.) In contrast to this, iterative optimization runs a coarse-grained grid of experiments that can help to identify contour lines which roughly delimit the area which simultaneously minimizes the amount of cache misses for all three caches (L1, L2, and TLB). However, even armed with that information, we can observe from the three cache miss plots that such region still comprises an enormous number of cases that nonetheless have to be individually evaluated to finally reveal the optimal combination for the architecture (mc=96 and kc=256 for this particular case). A third approach to identify the optimal parameter values is to use analytical models as advocated by some in the compiler community. Concretely, simple analytical models have been previously developed in the compiler domain in order to determine optimizing parameter values for tile sizes and unrolling factors [Bacon et al. 1994; Yotov et al. 2005; Kelefouras et al. 2014]. These models are based on the data access pattern and exploit a fair knowledge ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY. Analytical Modeling is Enough for High Performance BLIS A:5 of the hardware features commonly found in modern processor architectures. Previous work has shown that parameter values derived analytically can deliver performance rates that are similar to those attained using values determined through empirical search [Yotov et al. 2005]. More importantly, the time it takes to analytically identify the optimal parameter combination is often a fraction of that required by empirical search. Let us distinguish this from autotuning by calling it auto-fine-tuning. We adopt this in our quest for the parameter values that optimize the algorithm in BLIS, with the aim to make even auto-fine-tuning unnecessary for many architectures. 3. EXPERT IMPLEMENTATION OF GEMM The gemm algorithm embedded within the BLIS framework is described in [Van Zee and van de Geijn 2014], and the same approach is instantiated in BLAS libraries optimized by DLA experts such as GotoBLAS [Goto and van de Geijn 2008b; 2008a] and its successor, OpenBLAS [OpenBLAS 2015]. Importantly, BLIS exposes three loops within the inner kernel used by the GotoBLAS and OpenBLAS, which then facilitates the analytical determination of the parameters. For completeness, we provide a brief overview of the algorithm next. In addition, we identify the data usage pattern and highlight the parameters that characterize the algorithm. 3.1. An expert’s implementation of gemm Without loss of generality, we consider the special case of the matrix-matrix multiplication, C:= AB +C, where the sizes of A,B,Care m×k,k×n, and m×n, respectively. In the following elaboration, we will consider different partitionings of the m,nand k-dimensions of the problem. For simplicity, when considering a (vertical or horizontal) partitioning of one of the problem dimensions, say q, into panels/blocks of size (length or width) r, we will assume that qis an integer multiple of r. BLIS implements gemm as three nested loops around a macro-kernel2, plus two packing routines. The macro-kernel is implemented in terms of two additional loops around a microkernel. The micro-kernel is a loop around a rank-1 (i.e., outer product) update, and it is typically implemented in assembly or with vector intrinsics. In BLIS, the remaining five loops are implemented in C. Pseudo-code for the gemm algorithm is given in Figure 2. The outermost loop (Loop 1, indexed in the figure by jc) traverses the n-dimension of the problem, partitioning both Cand Binto column panels of width nc. The second loop (indexed by pc) processes the k-dimension, partitioning Ainto column panels of width kcand the current panel of Binto blocks of length kc. The third loop (indexed by ic) iterates over the m-dimension, yielding partitionings of the current column panels of Cand Ainto blocks of length mc. The following loops comprise the macro-kernel. Let us define Cc=C(ic:ic+mc−1, jc:jc+nc−1). Loop 4 (indexed by jr) traverses the n-dimension, partitioning Ccand the packed block Bcinto micro-panels of width nr. Loop 5 (indexed by ir) then progresses along the m-dimension, partitioning the current micro-panel of Ccinto micro-tiles of length mrand the packed block Acinto micro-panels of the same length. The innermost loop (Loop 6, indexed by pr), inside the micro-kernel, computes the product of a micro-panel of Acand micro-panel of Bcas a sequence of kcrank-1 updates, accumulating the result to a micro-tile of Cc, which can be described mathematically as Cc(ir:ir+mr−1, jr:jr+nr−1) += Ac(ir:ir+mr−1,0 : kc−1)·Bc(0 : kc−1, jr:jr+nr−1). At this point it is worth emphasizing the difference between Ccand Ac, Bc. While the former is just a notation artifact, introduced to ease the presentation of the algorithm, the latter two correspond to actual buffers that are involved in data copies. 2The macro-kernel is also known as the inner kernel in GotoBLAS. ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY. A:6 Tze M. Low et al. ic jc nc mc ic kc mc c p jc nc kc c pBC A Loop 1 for jc= 0,...,n−1in steps of nc Loop 2 for pc= 0,...,k−1in steps of kc B(pc:pc+kc−1, jc:jc+nc−1) →Bc// Pack into Bc Loop 3 for ic= 0,...,m−1in steps of mc A(ic:ic+mc−1, pc:pc+kc−1) →Ac// Pack into Ac Loop 4 for jr= 0,...,nc−1in steps of nr// Macro-kernel Loop 5 for ir= 0,...,mc−1in steps of mr Loop 6 for pr= 0,...,kc−1in steps of 1 // Micro-kernel Cc(ir:ir+mr−1, jr:jr+nr−1) += Ac(ir:ir+mr−1, pr) ·Bc(pr, jr:jr+nr−1) endfor endfor endfor endfor endfor endfor Fig. 2: High performance implementation of gemm by an expert. A B c c r m nr Fig. 3: Packing in the BLIS and GotoBLAS implementations of gemm. 3.2. The importance of packing It is generally known that accessing consecutive memory locations (also known as in-stride access) is usually faster than non-consecutive memory accesses. By packing elements of Ain Loop 3 and Bin Loop 2 in a special manner into Acand Bc, respectively, we can ensure that the elements of these two matrices will be read with in-stride access inside the micro-kernel. ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY. Analytical Modeling is Enough for High Performance BLIS A:7 Loop Reused data Size of reused data Reuse factor 6 Micro-tile, Crmr×nrkc 5 Micro-panel, Brkc×nrmc/mr 4 Packed block Acmc×kcnc/nr 3 Packed block Bckc×ncm/mc 2 Column panel of C m ×nck/kc 1 Matrix A m ×k n/nc Table I: Analysis of data reuse performed in the different loops of the BLIS implementation of gemm. Concretely, each mc×kcblock of Ais packed into Ac. Elements are organized as row micro-panels of size mr×kc, and within each micro-panel of Ac, the elements are stored in column-major order. Each kc×ncblock of Bis packed into Bc. In this case, elements are packed into column micro-panels of size kc×nr, and each column micro-panel is stored into row-major order. The reorganization of the entries of Aand Binto blocks of Acand Bcwith the packing layout illustrated in Figure 3 ensures that these elements are accessed with unit stride when used to update a micro-tile of C. Packing Acand Bcalso has the additional benefit of aligning data from Acand Bcto cache lines boundaries and page boundaries. This enables to use instructions for accessing aligned data, which are typically faster than their non-aligned counterparts. A third benefit of packing is that data from Acand Bcare preloaded into certain cache levels of the memory hierarchy. This reduces the time required to access the elements of Ac and Bcwhen using them to update a micro-tile of C. Since Acis packed in Loop 3, after Bchas been packed in Loop 2, the elements of Acwill likely be higher up in the memory hierarchy (i.e. closer to the registers) than those of Bc. By carefully picking the sizes for Ac and Bc, the exact location of Acand Bcwithin the memory hierarchy can be determined. 3.3. Data usage pattern Consider Loop 6. This innermost loop updates a micro-tile of C, say Cr, via a series of kc rank-1 updates involving a micro-panel Ar, from the packed matrix Ac, and a micro-panel Br, from the packed matrix Bc. At each iteration, mrelements from Arand nrelements from Brare multiplied to compute mrnrintermediate results that will be accumulated into the micro-tile of C. Between two iterations of Loop 6, different elements from Arand Br are used, but the results from the rank-1 updates are accumulated into the same micro-tile. Hence, the micro-tile Cris the data reused between iterations (Reused data) in Loop 6. In addition, because kcrank-1 updates are performed each time Loop 6 is executed, Cris said to have a reuse factor of kctimes. Since the columns of Arand the rows of Brare used exactly once each time through Loop 6, there is no reuse of these elements. Identifying pieces of data that are reused is important because it tells us which data portions should be kept in cache in order to allow fast access and reduce cache trashing. An analysis similar to that with Loop 6, performed on the remaining loops to identify the data that is reused in each one of them, the size of the reused data, and the number of times the data is being reused, yields the results summarized in Table I. Notice that the size of the reused data becomes smaller as the loop depth increases. This nicely matches the memory hierarchy, where faster memory layers (caches) feature lower capacity than slower caches. As such, gemm maps smaller reused data to faster layers in the memory hierarchy, as depicted in Figure 4. ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY. A:8 Tze M. Low et al. 0 0 0 1 1 1 00 00 11 11 00 00 00 00 11 11 11 11 00 00 11 11 Cr BA br BC Memory L3 cache L2 cache Registers L1 cache ar C Pack B Pack A micro−kernelmicro−kernel CC C A i = 0 r Load from Load from micro−kernel Load from micro−kernel Load from Fig. 4: Packing in the BLIS and GotoBLAS implementations of gemm. 3.4. Characteristic parameters that drive performance As described in the previous subsection and captured in Table I, the gemm algorithm that underlies BLIS is characterized by the following five parameters: mc, kc, nc, mrand nr, which correspond to the block/panel/tile size for each of the loops around the micro-kernel. In addition, this set of parameters also determines the dimension of the reused data and the reuse factor of each of the loops. Loop 6 is defined by three of the five parameters (mr, nrand kc), while Loops 3, 4 and 5 are characterized by four of the five parameters. In addition, three of the four parameters that characterize Loops 3, 4 and 5 carry beyond to the next inner loop as well. Now, when the parameters that characterize Loop 5 have been identified, the only unknown parameter for Loop 4 is nc. This observation suggests that the parameters should be identified from the innermost loop outwards. In the next section, we leverage this observation to optimize this collection of parameters analytically. ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY. Analytical Modeling is Enough for High Performance BLIS A:9 4. MAPPING ALGORITHM TO ARCHITECTURE An effective mapping of the micro-kernel and the loops around it to the architecture is fundamental to attain high performance with the BLIS approach. This means that we need to identify values for the characteristic parameters that are tuned for the target machine architecture. As discussed in the previous section, we start by mapping the micro-kernel to the architecture and work our way outwards in order to identify the optimal values for the different characteristic parameters. 4.1. Architecture model To develop an optimization model for mapping the gemm algorithm onto an architecture, we first need to define a model for our target architecture. We make the following assumptions regarding our hypothetical processor: —Load/store architecture and vector registers. Data have to be loaded into the processor registers before computation can be performed with them. A total of Nreg vector registers exist, where each can hold Nvec elements (floating-point numbers) of size Sdata. (In a machine with no vector registers, Nvec = 1.) In addition, we assume that memory instructions can be issued in parallel with floating-point arithmetic instructions. —Vector instructions. The processor floating-point arithmetic units have a throughput of Nvfma vector (or SIMD) fused multiply-add instructions (vfma) per clock cycle (i.e., 2NvecNvfma flops per cycle). A single vfma instruction combines Nvec regular (nonvector) instructions and produces Nvec scalar results. Furthermore, each vfma has a latency of Lvfma cycles, which is the minimum number of cycles between the issuance of two dependent consecutive vfma instructions. On an architecture without vfma instructions, Lvfma is computed by adding the latencies for a multiply and an addition instruction. —Caches. All data caches are set-associative, and each cache level Liis characterized by four parameters as follows: –CLi: size of cache line, –WLi: associativity degree, –SLi: size, and –NLi: number of sets, where SLi=NLiCLiWLi. A fully-associative cache can be modelled by setting NLi= 1 and WLi=SLi/CLi. The replacement policy for all caches is least-recently-used (LRU) [Hennessy and Patterson 2003]. For simplicity, we assume that the size of a cache line is the same for all cache levels. 4.2. Parameters for the body of the inner-most loop: mrand nr Recall that the micro-kernel is characterized by three parameters: mr,nrand kc, where the former two determine the size of the micro-tile of Cthat is reused across every iteration. In addition, these two parameters also define the number of elements of Arand Brthat are involved in the rank-1 update at each iteration of Loop 6. In this subsection, we discuss how mrand nrcan be identified analytically. 4.2.1. Strategy for finding mrand nr.The main strategy behind our approach to identify mr, nris to choose them “large enough” so that no stalls due to the combination of dependencies and latencies are introduced in the floating-point pipelines during the repeated updates of the micro-tile Cr. In addition, the smallest values of mr, nrwhich satisfy this condition should be chosen because this implies that Croccupies the minimum number of registers, releasing more registers for the entries of Ar, Br. This in turn enables larger loop unrolling factors and more aggressive data preloading [Hennessy and Patterson 2003] to reduce loop overhead. ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY. A:16 Tze M. Low et al. 24 120 216 312 408 504 600 696 792 888 984 1080 1176 1272 1368 1464 1560 1656 1752 1848 1944 16 18 20 22 24 mc GFLOPS BLIS on AMD Kaveri (m=n=k=4000). GFLOPS 24 120 216 312 408 504 600 696 792 888 984 1080 1176 1272 1368 1464 1560 1656 1752 1848 1944 2e+07 4e+07 6e+07 8e+07 1e+08 mc L2 misses BLIS on AMD Kaveri (m=n=k=4000). L2 misses Fig. 6: GFLOPS and L2 cache misses (left and right plots, respectively) on the AMD Kaveri , for a matrix multiplication with all three operands square of dimension m=n=k= 4,000, the optimal kc= 128, and varying mcin steps of 24. The blue and green markers respectively identify the empirical optimal value mE c= 1,368 and the analytical suboptimal ˜mA c= 1,788. The same comments of the previous case apply to the suboptimal case mc= 512 (red line), with the only exception of the high number of L1 cache misses when kcis very small. The reason is that, for this large value of mc, the packing of Acproduces this effect. 5.3.2. Quantification of the analytical model deviations. The AMD Kaveri case. Let us analyze the discrepancy between the analytical and empirical optimal values for mcon the AMD Kaveri. In particular, the model fixes the analytical optimal at mA c= 1,792, while the experimentation found the empirical optimal at mE c= 1,368. Unfortunately, BLIS enforces that mc is an integer multiple of both mrand nr(on this architecture, 4 and 6 respectively) so that the analytical optimal cannot be experimentally tested. To avoid exceeding the capacity of the cache, in the following analysis we will thus consider the analytical suboptimal value at ˜mA c= 1,788, which corresponds to the closest integer multiple of both 4 and 6 smaller than mA c. Figure 6 relates the performance (GFLOPS) and L2 cache misses on the AMD Kaveri, using the optimal value kc= 128, for a wide range of values for mcthat include both mE cand ˜mA c. The first point to note in the right-hand side plot is that the analytical suboptimal correctly determines the dimension from which the number of L2 cache misses initiates an exponential growth due to capacity constraints (specifically, this occurs from mc≥1,800 > mA c>˜mA c). Compared with the empirical optimal, the case ˜mA cincreases the L2 cache misses by 8.85% which, in principle, seems a large value but cannot be appreciated in the plot due to the scale. However, relative to the maximum number of cache misses in the plot (i.e, normalized against the misses for mc= 2,016), the misses for ˜mA conly incur an increase of 0.78% with respect to the figure for mE c. With these parameters, the slightly larger number of L2 cache misses results in a small decrease of performance: from 23.95 GFLOPS for mE cto 23.48 GFLOPS for ˜mA c, i.e. −1.97%. Let us consider now the empirical suboptimal ˜mE c= 1,656, which corresponds to the largest possible value of mcthat still is within the “noise” of the best performance. From a practical point of view, the difference in performance between mE cand ˜mEis negligible (23.946 GFLOPS for the latter, which represents a decrease of −0.017% with respect to mE cand can be considered to be in the noise level). Taking into account the dimensions of Ac, for ˜mE cthis block occupies 1,656 KBytes, i.e. 85.55% of the L2 cache or, in this 16-way associative cache, close to 13 lines of each 16-line set (the exact value is 12.93). We can observe that this is not far from the analytical optimal mA c, which dictates that ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY. Analytical Modeling is Enough for High Performance BLIS A:17 Acshould occupy 87.50% of the L2 cache, or 14 lines of each 16-line set. The conclusion from this analysis is that the discrepancies between the analytical values ˜mA c, mA cand the empirical values ˜mE c, mEseem to be due to problems with the alignment of Acin the 16way associative cache, probably related to the packing of Acproducing unexpected conflict misses, and not to a deviation in the analytical model. A second source of aligment problems may be that, for this architecture, nr(= 6) is not a power of 2. Therefore, it may occur that each micro-panel of Bcdoes not overwrite the previous micro-panel of this block, in turn causing conflict misses since the streaming of this matrix does not overwrite the micro-panel of Bin the L2 cache. 6. CONCLUSION AND FUTURE DIRECTIONS We have developed an analytical model that yields optimal parameter values for the gemm algorithm that underlies BLIS as well as other manually-tuned high performance dense linear algebra libraries. The key to our approach is the recognition that many of the parameters that characterize BLIS are bounded by hardware features. By understanding how the algorithm interacts with the processor architecture, we can choose parameter values that leverage, more aggressively, hardware features such as the multi-layered cache memory of the system. We compare the values obtained via our analytical approach with the values from manually-optimized implementations and report that they are similar if not identical. Unlike similar work that compares with ATLAS, this is the first paper, to the best of our knowledge, that demonstrates that an analytical approach to identifying parameter values can be competitive with best-in-class, expert-optimized implementations. This demonstrates that a high performance implementation for the gemm algorithm can be obtained without relying on a costly empirical search. We believe that more can be done in terms of enhancing the analytical model. One possible direction we are exploring is to extend the analytical model to more complicated linear algebra operations such as those in LAPACK [Anderson et al. 1999]. In general these operations, e.g. matrix factorizations, are implemented as blocked algorithms, and the gemm operation is often a sub-operation in the loop body. It would be interesting to determine how to identify optimal block sizes, given that the parameters that favor performance for the gemm operation are known. A second direction to enhance the analytical model is to incorporate data movement considerations. Currently, the parameters for the micro-kernel (mrand nr) are determined by only considering the latency of the floating-point arithmetic instructions. However, this makes the assumption that bandwidth is large enough. On low-power systems, this assumption may no longer hold. In such scenario, mrand nrmay have to be larger, so that the time for computation is sufficiently long to hide the time it takes to load the next elements of Arand Brinto the processor registers. Finally, our current analytical model assumes double precision arithmetic. With complex arithmetic, a micro-tile of the same size incurs four times as many flops, and involves twice as much data. These changes necessarily mean that the formulas in this paper have to be updated. Nonetheless, we believe that a similar analysis of the algorithm and the hardware features yields an analytical model for a micro-kernel that operates with complex arithmetic. A question becomes whether our analysis can be used to better design future architectures. This question is at least partially answered in [Pedram et al. 2012b; Pedram et al. 2012a], which examines how to design specialized hardware (both compute core and entire processor) for linear algebra computation. The models used for such purpose have much in common with our model for determining the parameter values for the micro-kernel. This suggests that our model for determining the parameter values for loops around the microkernel can potentially be leveraged to either determine the ideal cache size and/or cache replacement policy. ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY. A:18 Tze M. Low et al. Acknowledgments This research was sponsored in part by NSF grants ACI-1148125/1340293 and CCF0917167. Enrique S. Quintana-Ort´ı was supported by project TIN2011-23283 of the Ministerio de Ciencia e Innovaci´on and FEDER. Francisco D. Igual was supported by project TIN201232180 of the Ministerio de Ciencia e Innovaci´on. This work was partially performed during their visit to The University of Texas at Austin (UT), funded by the JTO visitor applications programme from the Institute for Computational Engineering and Sciences (ICES) at UT. REFERENCES AMD. 2015. AMD Core Math Library. http://developer.amd.com/tools-and-sdks/cpu-development/ amd-core-math-library-acml/. (2015). E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, Jack J. Dongarra, J. Du Croz, S. Hammarling, A. Greenbaum, A. McKenney, and D. Sorensen. 1999. LAPACK Users’ guide (3rd ed.). SIAM. David F. Bacon, Susan L. Graham, and Oliver J. Sharp. 1994. Compiler Transformations for High-performance Computing. ACM Comput. Surv. 26, 4 (Dec. 1994), 345–420. DOI:http://dx.doi.org/10.1145/197405.197406 G. Belter, J. G. Siek, I. Karlin, and E. R. Jessup. 2010. Automatic Generation of Tiled and Parallel Linear Algebra Routines: A partitioning framework for the BTO Compiler. In Proceedings of the Fifth International Workshop on Automatic Performance Tuning (iWAPT10). 1–15. Jeff Bilmes, Krste Asanovic, Chee-Whye Chin, and Jim Demmel. 1997a. Optimizing Matrix Multiply Using PHiPAC: A Portable, High-performance, ANSI C Coding Methodology. In Proceedings of the 11th International Conference on Supercomputing (ICS ’97). ACM, New York, NY, USA, 340–347. DOI:http://dx.doi.org/10.1145/263580.263662 Jeff Bilmes, Krste Asanovi´c, Chee whye Chin, and Jim Demmel. 1997b. Optimizing Matrix Multiply using PHiPAC: a Portable, High-Performance, ANSI C Coding Methodology. In Proceedings of International Conference on Supercomputing. Vienna, Austria. Jim Demmel, Jack Dongarra, Victor Eijkhout, Erika Fuentes, Antoine Petitet, Rich Vuduc, R. Clint Whaley, and Katherine Yelick. 2005. Self adapting linear algebra algorithms and software. In Proceedings of the IEEE. 2005. Jack J. Dongarra, Jeremy Du Croz, Sven Hammarling, and Iain Duff. 1990. A Set of Level 3 Basic Linear Algebra Subprograms. ACM Trans. Math. Soft. 16, 1 (March 1990), 1–17. Jack J. Dongarra, Jeremy Du Croz, Sven Hammarling, and Richard J. Hanson. 1988. An Extended Set of FORTRAN Basic Linear Algebra Subprograms. ACM Trans. Math. Soft. 14, 1 (March 1988), 1–17. Kazushige Goto and Robert van de Geijn. 2008a. High Performance Implementation of the Level-3 BLAS. ACM Trans. Math. Software 35, 1 (July 2008), 4:1–4:14. http://doi.acm.org/10.1145/1377603.1377607 Kazushige Goto and Robert A. van de Geijn. 2008b. Anatomy of a High-Performance Matrix Multiplication. ACM Trans. Math. Software 34, 3 (May 2008), 12:1–12:25. http://doi.acm.org/10.1145/1356052. 1356053 John A. Gunnels, Greg M. Henry, and Robert A. van de Geijn. 2001. A Family of High-Performance Matrix Multiplication Algorithms. In Computational Science - ICCS 2001, Part I (Lecture Notes in Computer Science 2073), Vassil N. Alexandrov, Jack J. Dongarra, Benjoe A. Juliano, Ren´e S. Renner, and C.J. Kenneth Tan (Eds.). Springer-Verlag, 51–60. John L. Hennessy and David A. Patterson. 2003. Computer Architecture: A Quantitative Approach. Morgan Kaufmann Pub., San Francisco. IBM. 2015. Engineering and Scientific Subroutine Library. http://www-03.ibm.com/systems/power/ software/essl/. (2015). Intel. 2015. Math Kernel Library. https://software.intel.com/en-us/intel-mkl. (2015). Vasilios Kelefouras, Angeliki Kritikakou, and Costas Goutis. 2014. A MatrixMatrix Multiplication methodology for single/multi-core architectures using SIMD. The Journal of Supercomputing (2014), 1–23. DOI:http://dx.doi.org/10.1007/s11227-014-1098-9 T. Kisuki, P.M.W. Knijnenburg, M.F.P. O’Boyle, and H. A. G. Wijshoff. 2000. Iterative Compilation in Program Optimization. (2000). Peter M. W. Knijnenburg, Toru Kisuki, and Michael F. P. O’Boyle. 2002. Iterative Compilation.. In Embedded Processor Design Challenges (Lecture Notes in Computer Science), Ed F. Deprettere, Jrgen Teich, and Stamatis Vassiliadis (Eds.), Vol. 2268. Springer, 171–187. http://dblp.uni-trier.de/db/conf/samos/ samos2002.html#KnijnenburgKO02 ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY. Analytical Modeling is Enough for High Performance BLIS A:19 C. L. Lawson, R. J. Hanson, D. R. Kincaid, and F. T. Krogh. 1979. Basic Linear Algebra Subprograms for Fortran Usage. ACM Trans. Math. Soft. 5, 3 (Sept. 1979), 308–323. OpenBLAS 2015. http://www.openblas.net. (2015). Ardavan Pedram, Andreas Gerstlauer, and Robert A. van de Geijn. 2012a. On the Efficiency of Register File versus Broadcast Interconnect for Collective Communications in Data-Parallel Hardware Accelerators. Computer Architecture and High Performance Computing (SBAC-PAD), 2012 24th International Symposium on (October 2012). Ardavan Pedram, Robert A. van de Geijn, and Andreas Gerstlauer. 2012b. Codesign Tradeoffs for HighPerformance, Low-Power Linear Algebra Architectures. IEEE Trans. Comput. 61 (December 2012), 1724–1736. DOI:http://dx.doi.org/10.1109/TC.2012.132 Tyler M. Smith, Robert van de Geijn, Mikhail Smelyanskiy, Jeff R. Hammond, and Field G. Van Zee. 2014. Anatomy of High-Performance Many-Threaded Matrix Multiplication. In IPDPS ’14: Proceedings of the International Parallel and Distributed Processing Symposium. To appear. Field G. Van Zee, Tyler Smith, Bryan Marker, Tze Meng Low, Robert A. van de Geijn, Francisco D. Igual, Mikhail Smelyanskiy, Xianyi Zhang, Michael Kistler, Vernon Austel, John Gunnels, and Lee Killough. 2014. The BLIS Framework: Experiments in Portability. ACM Trans. Math. Soft. (2014). In review. Field G. Van Zee and Robert A. van de Geijn. 2014. BLIS: A Framework for Generating BLAS-like Libraries. ACM Trans. Math. Soft. (2014). To appear. Qian Wang, Xianyi Zhang, Yunquan Zhang, and Qing Yi. 2013. AUGEM: Automatically Generate High Performance Dense Linear Algebra Kernels on x86 CPUs. In Proceedings of SC13: International Conference for High Performance Computing, Networking, Storage and Analysis (SC ’13). ACM, New York, NY, USA, Article 25, 12 pages. DOI:http://dx.doi.org/10.1145/2503210.2503219 R. Clint Whaley and Jack J. Dongarra. 1998. Automatically Tuned Linear Algebra Software. In Proceedings of SC’98. Kamen Yotov, Xiaoming Li, Mar´ıa Jes´us Garzar´an, David Padua, Keshav Pingali, and Paul Stodghill. 2005. Is Search Really Necessary to Generate High-Performance BLAS? Proceedings of the IEEE, special issue on “Program Generation, Optimization, and Adaptation” 93, 2 (2005). ACM Journal Name, Vol. V, No. N, Article A, Publication date: January YYYY.