Artículo presentado al Congreso Español de Informática 2013Performance Analysis of the Multi-pass Transformation for Complex 3D-Stencils on GPU

López-Zapata, Emilio

Romero, Luis F.

Tabik, Siham

English

Repositorio Institucional Universidad de Málaga

1Performance Analysis of the Multi-passTransformation for Complex 3d-Stencils onGPUsS. Tabik 1, L. F. Romero 1, E. L. Zapata 1Abstract|Complex iterative 3d stencils based on aseries of multiple simpler stencils with dierent com-putation intensities cannot be handled properly usingstandard techniques on the GPU. This work demon-strates that decomposing these kind of stencils into asequence of up to a specic number of simpler stencilsand further optimizing each individual kernel providesthe best overall performance. We focus on the familyof PDE-based denoising methods, which can be refor-mulated as sequence of multiple stencils-based tasks.The performance results and analysis show that thereexists an optimal level of splitting-coalescence of thesestencils-based tasks that reaches the best compromisebetween better use of fast-memories and higher con-currency.Keywords| complex stencils, decomposing-coalescing tasks, GPU accelerators, multi-passoptimization.I. Motivation and Related WorksTHERE exists a large number of works that fo-cuses on optimizing naive 2d and 3d stencil ker-nels on multi-core and many-core systems. Theseworks can be divided into two classes. The rst classfocuses on applying transformations such as ghost-zone, time-tiling, or a mix of both [9], [7], [5]. Whilethe second category focuses on building auto-tuningenvironments which determine the subset of opti-mizations that improves the performance of the sten-cil on each specic architecture [3].Using ghost-zone optimization, which implies en-larging the size of the halo and performing redun-dant computation, reduces communications betweentiles on the GPU but only for iterative 2d stencilswith low computation intensity [7], [9]. Time tiling,which consists of tiling the matrix along the time di-mension, enhances data locality but again only forsimple memory-bound iterative 2-d stencil codes [5],[9].All these standard optimizations are essential toimprove the performance of very simple stencils butthey do not require or simply are not applicable tocomplex stencils, such as iterative compute-intensive3d-stencils.On the other hand, the multi-pass optimization,which consists of decomposing the CUDA-kernel intomultiple CUDA-kernels, has not been analyzed in awider context. There exist few works in the litera-ture that employed this optimization to improve aspecic application. Indeed, none of them has an-alyzed at what extent it is benecial to the perfor-mance for a whole family of algorithms. For instance,1Dpto. de Arquitectura de Computadores, Univ. Malaga,e-mail: stabik@uma.esLee et al. claimed that applying two-passes ap-proach to a memory-bound neuroimaging algorithmincreases concurrency. While, Abdelfattah et al. [1]used two-passes optimization to separate two stageswith dierent computation patterns in the symmet-ric matrix-vector product code.The main contributions of this work are: 1) Un-derstanding when the multi-pass optimization can bebenecial to complex stencils representative of PDE-based denoising methods and 2) presenting an inter-esting pattern of PDE-denoising methods that sub-stantially helps optimizing them for both many andmulti-cores.II. PDE-denoising methods: ComplexStencilsDenoising algorithms are essential tools in im-age processing in many bioimaging modalities. Themost sophisticated and powerful methods are thosethat solve Partial Dierential Equations (PDE). Thisfamily of algorithms, like for instance the ones de-scribed in [2], [4], [6], apply an iterative sequenceof mathematical kernels, formulated as stencils withdierent computation intensity, where the input im-age is partially or entirely read and updated severaltimes during each denoising iteration till obtainingthe nal ltered 3d-image after a total number ofabout 60 to 100 iterations.A. A Case Study: Anisotropic Nonlinear DiusionAlgorithmFor simplicity, let's focus on a particular casestudy, the Anisotropic Nonlinear Diusion (AND)algorithm to understand the features of the consid-ered family of methods. AND-algorithm can be for-mulated as a stream of four stages, where the noisyimage travels through these four methods of dierentdata dependencies, data access patterns and arith-metic intensity. That is, each single iteration consistsof: First a Gaussian smoothing is applied to the ini-tial 3d-image. Each pixel is updated using a con-volution along the X-, Y- and Z-axes, i.e., usingthe six nearest neighbors. The update of eachpixel needs 12 ops (6 reads+1 write). Which iscomparable to a 7-points stencil read from andwrote to the 3d-image as shown in Figure 1(a). Then, the structure tensor, ST, of the obtained3d smoothed image is calculated. Each pixelof ST, which consists of a structure of 6 oats,2YZk-1kk+1XYZk-1kk+1X(a) (b)Fig. 1. (a) 7-point stencil, i.e., each pixel needs 7 nearest neighbors including itself to be updated. (b) 19-point stencil, i.e.,each pixel needs 19 nearest neighbors to be updated.is calculated using its corresponding six nearestneighbors from the smoothed image. The dataaccess pattern is a 7-points stencil read from the3d-image of size Nx Ny Nz and wrote to ST ofsize Nx Ny Nz  6. Afterwards, the Diusion Tensor, DT, is calcu-lated. Each pixel of DT is calculated by cal-culating the eigenvalues and eigenvectors of thecorresponding symmetric 3x3-matrix using theiterative Jacobi method. The elements of thesymmetric 3x3-matrix are initialized using thevalues of the corresponding ST pixel. This iscomparable to 1-point stencil read and wrote toST. This stage is clearly compute intensive. Finally, each pixel of the 3d-image is updatedusing 7 neighbor points from ST and 19 neighborpoints from the 3d-image. Which is a multiple7-points and 19-points stencil that reads datafrom the 3d-image and ST and writes the resultin the 3d-image. Updating each pixel of the 3d-image needs 26 reads, one write and 68 ops.III. Splitting-Coalescence of Stencils onGPU: Multi-pass ApproachCommonly, parallel computing PDE-based meth-ods on clusters of muli-processors employs the SPMDmodel, where each processor denoises its own 3d-block of the image during the whole iterative process.At the end of each iteration processors that work onneighbor blocks intercommunicate their halos [10].An intuitive CUDA implementation of this fam-ily of methods can consider either encapsulating thewhole complex stencil into one pass or splitting itinto multiple passes. The performance of the onekernel version can be delimited by the smaller size offaster memory resources during the kernel executionwhile the multiple kernels version may be penalizedby higher scheduling overheads and trac betweeno-chip and on-ship memories, which is necessary forinter-kernels communication. However, more com-binations can be also considered by merging dier-ent successive stencils, 2'-passes and 3'-passes. Oneof the objectives of this work is to analyze to whatextent one should split or coalesce these stencils toobtain the best overall performance on GPU accel-erators.We developed six implementations of AND-method by spiting it into one-, two-, three- and fourpasses. Applying two- and three-passes optimizationhas two possible versions. Stencil computation isperformed along the z-dimension using space-tiling.First, the 2D tiles and necessary halos are loaded intothe shared memory, where each threadblock com-putes its corresponding output 2D tile, similar to thestrategy employed in [8]. Here is a brief descriptionof each implementation. The 1-pass approach, encapsulates the wholemethod, Gauss smoothing, structure tensor cal-culation, diusion tensor calculation and, thePDE solution, i.e., Gauss+ST+DT+PDE, intoone pass. Updating one 2D tile of the output 3d-image needs to upload to shared memory 4 tilesfrom the original image (i.e., the correspondinginput tile + 3 halos) and 3 tiles of ST (i.e., thecorresponding ST tile to be calculated + 2 ST-tile-halos). Recall that three tiles of the original3d-image are required to compute one ST tile.The 2d-tiles with Z = 0 and Z = Nz of theoutput 3D-image and ST are calculated using amirroring operation of the 2D-plan with Z = 2and Z = Nz   2 respectively. The 2-passes implements Gauss+ST+DT in onepass and PDE in a second pass. The rst passupdates each ST-tile using 1 tile of the input3d-image; the halo pixels are uploaded into reg-isters. The PDE pass needs one ST-tiles and 3tiles of the input 3d-image; the halo pixels areuploaded into registers. The 2'-passes implements Gauss+ST in one passand DT+PDE in a second pass. The rst passupdates each ST tile using one tile of the input3d-image; the halo pixels are uploaded into reg-isters. The DT+PDE pass needs 3 ST-tiles and4 tiles of the input 3d-image.3ImgiNxNyNzImgi+1NxNyNzST            PDESTNxNyNz*6ImgiNxNyNz Gauss            DT            STNxNyNz*6Fig. 2. Pattern of one AND-method iteration: A sequence of four stencil-based stages (in pink color). Gauss() smooths theimage, ST() computes the structure tensor, DT() computes the diusion tensor and PDE() update the image The 3-passes version implements Gauss+ST ina rst pass, DT in a second pass and, PDE in athird pass. The rst pass needs to upload intoshared memory one 2D-tile of the input imageand one tile of its corresponding ST tile; the halopixels are uploaded into registers. The 3'-passes approach implements Gauss in arst pass, ST+DT in a second pass and, PDEin a third pass. The ST+DT pass needs one tileof 3d-image and one ST-tile. Finally, the 4-passes version implements eachsingle stage into one single kernel. Gauss passneeds one 2D-tile of the input image. ST passneeds one tile of ST and one tile of the 3D-image.DT pass needs only to update one ST tile.IV. Experimental Results and DiscussionThis section provides and discusses the perfor-mance results of the AND-complex stencil whendecomposed into one, two, three and four sim-pler stencils, each pass is encapsulated into oneCUDA kernel, and its comparison with the Pthreads-implementation described in [10]. The evaluationsof the CUDA-implementations, written in CUDAC v4.0, were carried out on a single Fermi C2070,with 14 multiprocessors, 448 cores, 48 Kb of sharedmemory, 32768 registers per block and 1.25 GB ofDRAM. While the Pthreads-implementation of thecode was executed on 8-core Intel X7550 Beckton us-ing 8threads. The Pthreads-implementation is usedonly for comparison purposes .Figure 3 (a) shows the speedup of the dierentCUDA-passes implementations (i.e., 1-, 2-, 2', 3-, 3'-and 4-passes codes) with respect to the Pthreads-implementation for three 3d-images of sizes 128 128  128, 256  256  256 and, 512  512  512.The found optimal thread-block sizes for the CUDA-codes are 4  32, 8  32, 16  64, 16  16, 16  64and, 16 64 respectively.Analyzing the speedup from the 1-pass to the 4-pass versions shows that the 3-passes version, whichimplements Gauss+ST, DT and PDE into one passeach one, is the fastest version over all the im-plementations providing a speedup with respect toPthreads-implementation of up to 290 for the smallersize and 55 for the larger image. In addition, Fig-ure 4(a),(b),(c) shows that the 3-passes version over-comes all the passes versions for the three image sizesby up to 160 and 40 for the smaller and larger im-age respectively.A further analysis of registers and shared memoryutilization and performance proling of each pass ofthe six codes are provided in Table 1 and 2 respec-tively. Examining the results of the 3-passes codeand comparing it with the other versions shows thatimplementing DT into an independent kernel allowsthis pass reaching an optimal utilization of on-shipand o-ship memories. In particular, it has a betterregisters and shared memory utilization, i.e., it ap-plies less stress on this fast memories, from 63 (whichis the maximum number of registers that can be usedper threads) to 42, in comparison with the pass thatmerges ST+DT from the 3'-passes version and theone that merges DT+PDE in the 2'-passes version.In addition, implementing DT into one pass reachesa better use of o-ship memories by showing bet-ter L1 and L2 hit rates. Moreover, in this 3-passesversion DT reaches a higher concurrency, i.e., highernumber of active warps per cycle from 15 to 22 whencomparing 3'-passes version versus the 3-passes onefrom Table 2. Notice that DT represents more than40% of the total runtime of all the implementationsand therefore increasing its performance aects theoverall performance.Notice from Table 1 and 2 that merging ST withGAUSS does not aect the concurrency of ST pass.This means that the overhead produced by separat-ing ST and GAUSS into two dierent passes is higherthan the fact of sharing the same fast memories.On the other hand, the performance of the PDEpass is limited by registers as can be seen from Ta-ble 1 and it provides better performance when it doesnot have to share fast memories.In summary and from a programming point ofview, an eective use of the multi-pass approach aswe learnt from the considered case study consists of1) implementing stencils limited by registers into onepass, 2) merging stencils that have the same data ac-cess pattern and use the same data structures intoone pass and 3) implementing the most time con-suming stencils with high concurrency grad into onepass because merging it with passes with dierentdata utilization penalizes its concurrency.4TABLE IRegisters and shmem utilization per thread (estimated using nvcc 4.0) for 256256256 3D-image1-pass 2-passes 2'-passes 3-passes 3'-passes 4-passesReg count 63 63,59 39,63 39,42,63 26,61,63 26,38,42,63shmem 29736 41616,32376 9080,29728 9080,7776,11664 1304,9072,11664 4632,32368,27744,41616(bytes)0501001502002501-pass 2-pass 2'-pass 3-pass 3'-pass 4-passspeedup (TPthreads/TCUDA)128x128x128256x256x256512x512x5120,E+005,E+081,E+092,E+092,E+093,E+093,E+09#pixels/s128x128x128256x256x256512x512x512(a) (b)Fig. 3. (a) Speedup of the 1-, 2-, 2'- 3-, 3'- and 4- passes CUDA-versions (on Tesla C5070) with respect to the Pthreads-version(on a 8-core Intel X7550 Beckton) for three 3d-image sizes, 128 128 128, 256 256 256 and, 512 512 512. (b) Thenumber of processed pixels per second by the Pthreads- and CUDA-, 1-, 2-, 2'- 3-, 3'- and 4- passes versions respectively.4x 3x160x 158x66x1511011512-pass 2'-pass 3-pass 3'-pass 4-passT1-pass/Ti-pass128x128x1284x2x50x37x 38x1112131412-pass 2'-pass 3-pass 3'-pass 4-passT1-pass/Ti-pass 256x256x256(a) (b)4x 2x40x32x 32x1112131412-pass 2'-pass 3-pass 3'-pass 4-passT1-pass/Ti-pass512x512x512(c)Fig. 4. (a) Speedup of 2-, 3-, 3'- and 4- passes versions with respect to the 1-pass implementation, on Tesla C5070, for three3d-image sizes, (a) 128 128 128, (b) 256 256 256 and, (c)512 512 512.V. ConclusionsThis paper characterizes PDE-denoising images tomake their CUDA-parallel implementation as e-cient as possible with less programming eorts. Weshowed that 1) implementing register limited sten-cils into one single pass, 2) merging stencils that usethe same data structure and data access into onepass and 3) implementing time consuming stencilswith high concurrency into one pass improves sub-stantially the performance of complex stencils thatare initially bounded by registers and shared mem-5TABLE IIOne(Gauss+ST+DT+PDE), two(Gauss+ST+DT PDE), three'(Gauss+ST, DT, PDE) and (Gauss, ST+DT, PDE)and, four-passes(Gauss, ST, DT, PDE) for a 256256256 3D-image and a fix thread-block 1616time active warps/ L1 gld IPC inst/ L2 glb mem read L2 glb mem write(%) active cycles hit rate byte throughput throughput(%) (Gb/s) (Gb/s)1-pass:Gauss+ST+DT+PDE 83,71 7,98 85,28 0,47 3,52 16,4 16,472-passess:Gauss+ST+DT 44,64 16,02 46,24 0,84 1,34 67,69 84,3PDE 7,95 16,46 86,39 1,06 4,51 47,16 9,792'-passess:Gauss+ST 2,56 16,79 66,18 0,36 1,53 9,20 52,78DT+PDE 80,11 4,21 84,36 0,55 13,81 5,24 5,253-passes:Gauss+ST 9,8 23,74 69,82 0,37 1,52 9,18 53,42DT 47,3 22,04 72,46 0,67 1,23 70,36 67,75PDE 6.08 14,89 86,42 1 4,65 45,61 9,813'-passes:Gauss 1,44 31,81 69,47 1,15 2,06 84,11 52,87ST+DT 43,76 15,92 40,92 0,83 1,33 68,17 84,66PDE 7,76 15,58 86,26 1,05 4,64 45,74 9,854-passes:Gauss 1,06 30,34 68,29 1,13 2,03 83,60 52,66ST 9,49 23,70 67,26 0,33 1,43 8,91 51,87DT 44,5 22,59 72,59 0,67 1,22 70,54 67,81PDE 5,7 16,51 85,69 1,13 4,89 45,64 9,83ory. We also proposed a way of systematizing theapplication of this optimization to diminish the pro-gramming time and eort.References[1] Dongarra J. Keyes D. Abdelfattah, A. and H. Ltaief. Op-timizing memory-bound numerical kernels on gpu hard-ware accelerators, 10th international meeting on high-performance computing for computational science (vec-par 2012), riken advanced institute for computational sci-ence (aics), kobe, japan, july 17th-20th. 2012.[2] D. Barash. Fundamental relationship between bilateralltering, adaptive smoothing, and the nonlinear diusionequation. Pattern Analysis and Machine Intelligence,IEEE Transactions on, 2002.[3] K. Datta, M. Murphy, V. Volkov, S. Williams, J. Carter,L. Oliker, D. Patterson, J. Shalf, and K. Yelick. Sten-cil computation optimization and auto-tuning on state-of-the-art multicore architectures. In Proceedings of the2008 ACM/IEEE conference on Supercomputing, SC '08,2008.[4] J.J. Fernandez. Tomobow: feature-preserving noise l-tering for electron tomography. doi:10.1186/1471-2105-10-178. 2009.[5] J. Holewinski, L. Pouchet, and P. P. Sadayappan. High-performance code generation for stencil computations ongpu architectures. In Proceedings of the 26th ACM inter-national conference on Supercomputing, ICS '12, pages311{320. ACM, 2012.[6] S. Li J.J. Fernandez. Anisotropic nonlinear ltering ofcellular structures in cryo-electron tomography. Com-puting in Science & Engineering, 7(5):5461, 2005.[7] J. Meng and K. Skadron. Performance modeling andautomatic ghost zone optimization for iterative stencilloops on gpus. In Proceedings of the 23rd internationalconference on Supercomputing, ICS '09, pages 256{265,New York, NY, USA, 2009. ACM.[8] P. Micikevicius. 3d nite dierence computation on gpususing cuda. In Proceedings of 2nd Workshop on Gen-eral Purpose Processing on Graphics Processing Units,GPGPU-2, pages 79{84, 2009.[9] G. Rivera and Chau-Wen Tseng. Tiling optimiza-tions for 3d scientic computations. In Supercomputing,ACM/IEEE 2000 Conference, page 32, 2000.[10] S. Tabik, E. M. Garzon, I. Garca, and J. J. Fernandez.High performance noise reduction for biomedical multi-dimensional data. Digit. Signal Process., 17(4):724{736,July 2007.

Performance Analysis of the Multi-pass Transformation for Complex 3D-Stencils on GPUs

https://riuma.uma.es/xmlui/bitstream/10630/5894/1/paper-2.pdf

Performance Analysis of the Multi-pass Transformation for Complex 3D-Stencils on GPUs

Abstract

Similar works

Full text

Available Versions

Repositorio Institucional Universidad de Málaga