Repositorio Institucional de Documentos

CUDA implementation of the solution of a system of linear equations arising in an hp-Finite Element code. Author: Javier Osés Villanueva 1+ , Directors: David Pardo 2* , and Ricardo Celorrio 3+ * Departamento de Matemática Aplicada y Estadística e I.O., University of the Basque Country UPV/EHU, and Ikerbasque. + Departamento de Matemática Aplicada, University of Zaragoza. February 1, 2013 1 ja[email protected] 2 [email protected] 3 [email protected] Abstract The FEM has proven to be one of the most ecient methods for solving dierential equations. Designed to run on dierent computer architectures, technological improvements have led over the years to the fast solution of larger and larger problems. Among these technological improvements, we emphasize the development of GPU (Graphic Processor Unit). Scientic programming in graphics cards was extremely dicult until 2006 the company NVIDIA developed CUDA (Compute Unied Device Architecture). It is a programming language designed for generic computing which does not require knowledge of traditional graphics programming. GPUs are capable of performing a large number of operations simultaneously. This capability makes them very attractive for use in FEM. One of the parts of the FEM which requires large computational capacity is the solution of systems of linear equations. In this work, an algorithm for solving systems of linear equations in CUDA has been implemented. It will be applied as a part of a hp-FEM code that tries to solve Laplace equation. The aim of this study is to compare the performance of an an implementation of a solver in CUDA vs. a C implementation and check if CUDA has advantages over traditional programming. For that purpose, we select an algorithm suitable for GPU programming. The iterative algorithms have properties that ts to CUDA programming architecture. However, the use of these algorithms require from double precision arithmetic to minimize round-o eects. Nowadays, only high performance GPUs are able to work in double precision. FEM matrices are sparse and the use of compression format for the system matrix is needed. Exist multiple compression formats and we select one which better ts to the matrix structure that FEM generates in our problem. The implementation in CUDA introduces improvements in execution times compared to traditional programming in C. Recent works has proved that it can be obtained programs that works until 80 times faster. But, this result can not be generalized because the improvements depends on dierential equation, boundary conditions, mesh generation, FEM, model of GPU, version of CUDA(now 5.0), and of course implementation. Contents 1 Introduction 3 2 Parallel programming in CUDA 6 2.1 Parallel computing . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.1 AmdahlLaw ....................... 7 2.2 CUDA............................... 8 2.2.1 Basic CUDA concepts . . . . . . . . . . . . . . . . . . 8 2.2.2 GPU architecture . . . . . . . . . . . . . . . . . . . . . 9 2.2.3 Information ux in GPU . . . . . . . . . . . . . . . . . 12 2.2.4 Programming features in CUDA . . . . . . . . . . . . 14 3 Model problem and Finite Element Formulation 17 3.1 FEM................................ 17 3.1.1 Variational Formulation of Laplace Equation. . . . . . 18 3.1.2 hp-FEM.......................... 20 4 Linear equation solvers 21 4.1 Gaussian elimination . . . . . . . . . . . . . . . . . . . . . . . 21 4.2 LU decomposition . . . . . . . . . . . . . . . . . . . . . . . . 22 4.3 Cholesky decomposition . . . . . . . . . . . . . . . . . . . . . 23 4.4 Conjugate Gradient . . . . . . . . . . . . . . . . . . . . . . . . 24 4.5 Preconditioned Conjugate Gradient . . . . . . . . . . . . . . . 25 4.5.1 Jacobi Preconditioner . . . . . . . . . . . . . . . . . . 25 4.5.2 SSOR Preconditioner . . . . . . . . . . . . . . . . . . 26 4.6 The most suitable algorithm . . . . . . . . . . . . . . . . . . . 26 5 Implementation 28 5.1 Algorithms ............................ 28 5.1.1 CG Algorithm . . . . . . . . . . . . . . . . . . . . . . 28 5.1.2 PCGAlgorithm...................... 29 5.2 Sparse matrix formats . . . . . . . . . . . . . . . . . . . . . . 30 5.3 Basic operations for CG and PCG . . . . . . . . . . . . . . . 31 5.3.1 Dotproduct........................ 31 5.3.2 Vector addition . . . . . . . . . . . . . . . . . . . . . . 33 1 CONTENTS 2 5.3.3 Matrix vector multiply in CRS . . . . . . . . . . . . . 34 5.3.4 SAXPY (Single-precision real Alpha X Plus Y) . . . . 35 5.4 Structure of the code . . . . . . . . . . . . . . . . . . . . . . . 36 6 Results 37 7 Conclusions and future work 38 Chapter 1 Introduction The nite element method FEM is a numerical technique for nding approximate solutions of partial dierential equations associated with physical problems on dierent types of geometry. Since its inception in 1950 until today, the use of FEM has spread of continuum mechanics to many elds such as heat transfer, uid mechanics, and even the study of biological systems. The FEM converts a problem dened in terms of linear second order partial dierential equations in a linear system of equations. For a basic introduction on FEM, see [53] [28]. Currently exist dierent types of FEM, including: AEM(Applied Element Method) GFEM(Generalized Finite Element Method) hp-FEM hpk-FEM XFEM(Extended nite element method) S-FEM(Smoothed nite element method) Spectral methods Meshfree methods Discontinuous Galerkin methods Finite element limit analysis Stretched grid method. 3 CHAPTER 1. INTRODUCTION 4 In this master thesis, we chose hp-FEM [12]. Although its computational implementation is somewhat challenging, it is currently one of the most ecient methods due to its adaptability and the fast convergence of their solutions. From a computational point of view, the problem with FEM solution consists of the following tasks: 1.- Pre-processing: 1.-Denition of geometry. 2.-Mesh generation. 3.-Boundary conditions. 2.- Calculation: 1.-Generation of the basis functions. 2.-Numerical integration. 3.-Solving the system of linear equations. 3.- Post-processing: 1.-Determination of approximation errors. The steps that have greater computational cost are the numerical integration and the solution of the system of linear equations. In this master thesis, we will focus on the ecient solution of systems of linear equations. There are two types of solvers: direct and iterative. Direct solvers have explicit expressions for the solution. These methods provide an exact solution (up to round-o errors). In contrast to this, iterative methods calculate an approximate solution in every step and the accuracy of the solution increases with the number of the iterations, achieving a superior computational speed than direct solvers. An advantage of the iterative solvers is that can provide suciently accurate solutions in a reduced number of iterations. A disadvantage of iterative methods that the former only over direct methods is to calculate approximations to the solution. The hp-FEM method applied to Laplace equations generates symmetric positive denite matrices with a specic pattern of sparsity. In practical applications the dimension of the stiness system can be quite large. Systems of dimension 1000000 are commonly solved by commercial software. Larger ones, up to tens or hundreds of millions, are being solved on supercomputers. To solve large systems it is necessary to have a high computational capacity and a code to run in parallel. Initially, this capability was only available CHAPTER 1. INTRODUCTION 5 to large computer centers and parallel programming was extremely complex. Today, this has changed due to technological advances in GPUs. Initially designed as graphics accelerators were able to perform very specic tasks [49]. Now, however are programmable devices with capacity to perform many operations simultaneously. In 2006, NVIDIA, one of the leading companies in the development of graphics cards created CUDA. This new language makes it easier to parallel program in parallel using GPUs because it was created as a extension of C language. For more information see [33,3537,39]. The main objective of this master thesis is to implement a solver of systems of linear equations in CUDA in a hp-FEM method applied to Laplace equation. A second objective is to determine whether the implementation of this program in CUDA provides signicant advantages versus its implementation in sequential CPU. In the last years have appeared many works related with GPUs [5], CUDA, FEM, and linear equation solvers. But commonly each publication focuses in very specic themes [13, 15, 22]. However, recently works have been published in which CUDA it is commonly used in FEM methods [22, 34, 42, 43]. In contrast, this master thesis focuses in the CUDA implementation of an iterative solver in a hp-FEM method, which is undeveloped. In the second chapter we are going to expose the basic ideas about CUDA. In the third chapter we will explain de model problem and the FEM formulation applied to Laplace equation. In the forth chapter we will see the most common solvers. In the fth chapter it is explain how have been implemented the algorithm. The sixth chapter shows the results obtained. In the last chapter we will describe the conclusions. Chapter 2 Parallel programming in CUDA The rst section of this chapter discusses basic ideas on parallel computing. The second section explains what is CUDA and its most important characteristics. 2.1 Parallel computing The most basic unit of operation is a transistor, which is able to perform operations in binary. Its main characteristics are size, operation frequency, power consumption, and heat dissipation. From the beginning of computers, the objective has been to increase the number of operations that an integrated circuit can perform. Normally, the two ways to increase the number of operations per unit time are to increase the number of transistors or to increase the operating frequency. As technology has evolved, the transistors have been reduced in size and the operating frequencies are getting higher. According to Moore's law [47] [29], we can expect that the number of transistors on an integrated circuit is doubles every 18 months. Current technology allows to work with frequencies from 1.5 to 4 GHz. A common size of one transistor is measured in nanometers and an integrated circuit can easily have hundreds of millions of transistors. This tendency cannot be sustained indenitely because the current silicon technology has physical limits on the scale and frequency of operation. It also happens that grouping large amounts of transistors increases heat generation. Therefore, it is necessary to have cooling systems in order to ensure correct functioning of the circuits, which also increases the power consumption of the device. The hardware limitations have forced developers to nd another way to increase performance. The solution is easy, use many computers working together and simulta- 6 CHAPTER 2. PARALLEL PROGRAMMING IN CUDA 7 neously. But to perform this task requires the development of dierent kind of software. Traditionally, software has been created for serial computing. To solve a problem, we construct an algorithm that it is implemented in a serial instruction stream. These instructions are executed in the central processing unit of a computer (CPU). When an instruction is completed, the following one is executed. Parallel computing (see [17] [44]) is a programming technique in which many instructions are executed simultaneously. A great way to deal with problems is to divide it into smaller problems that can be solved simultaneously. There are dierent types of parallel computing, according to the instructions to be parallelized: bit-level operations, instructions, data or tasks. Parallel computers can be classied according to their hardware into two groups. The rst group is composed of multicore and multiprocessing computers with multiple processing elements in a single machine, see( [32] [9] [21] [8]). The second group are the Clusters, the MPP(Massively Parallel Processing) and Grids that use multiple computers to work on the same task. The GPUs are in the rst group (see [41], [19], [50], [54]). 2.1.1 Amdahl Law As explained above, we can assume that the ability to parallelize can be used to improve processing speed indenitely. But it is not. To increase speed by parallelization of a code, one needs to know which parts of this can be parallelized. Amdahl's Law gives us a way to calculate the maximum improvement of a code when a part of this is improved. [24] A=1 (1 −Fm) + Fm Am (2.1.1) A is the acceleration or velocity gain achieved in the entire programming code due to the improvement of one of its sub-codes. Am is the acceleration or velocity gain achieved into the sub-codes improved. Fm is the fraction between the time of execution of the sub-code improved and the time of execution of the complete code. In simple terms, Amdahl's Law says that it is the algorithm the one that decides the speed improvement not the number of processors. At the end, you reach a situation in which the algorithm cannot be parallelized anymore. [1] CHAPTER 2. PARALLEL PROGRAMMING IN CUDA 14 Fifth Step: Executing Kernel in the device Before running the kernel, you need the number of threads and blocks that will need to perform the operation. We create two variables of type dim3(three dimensional arrays). If we are to operate with vectors of size DIM and blocks the size MAX_BLOCK_SIZE, we need at least DIM/MAX _ BLOCK _ SIZE+ 1 blocks and DIM threads. 1 dim3 blocks ((DIM/MAX_BLOCK_SIZE)+1 ,1 ,1) ; dim3 threads (MAX_BLOCK_SIZE,1 ,1) ; After that, we can invoke the kernel: Kernel_name<<<blocks , threads>>>(kernel_arguments ) ; Sixth Step: Copying values from the device to the host With the same function of the forth step we can copy the result of the kernel from the device to the host. 1 cudaMemcpy(h_Vector , d_Vector ,10 ∗ sizeof(float) , cudaMemcpyDeviceToHost) ; h_Vector :destination d_Vector :origin 10*sizeof(oat ): size of memory to copy cudaMemcpyDeviceToHost: type of instruction Seventh Step: Free memory from the device After nishing the kernel execution, it is necessary to free the GPU memory. This is done with the instruction cudaFree(pointer). It is also a good practice to free the CPU memory. 1 cudaFree (d_Vector) ; free (h_Vector) ; 2.2.4 Programming features in CUDA In summary, in order to make a program in CUDA, we want algorithms to exhibit the following features: Parallelizable: All algorithms should follow a sequence of steps and in general, not all the steps can be parallelized. Typically, algorithms with many conditionals and perform operations when the dierence between the number of inputs and outputs are very large. One example is the dot product of two vectors. This operation has 2N inputs (for dimension N) and only one output value. CHAPTER 2. PARALLEL PROGRAMMING IN CUDA 15 Simplicity in the operations: We have seen in the previous sections that one SM contains 8 SP an 2 SFU. In one clock cycle the SM can perform eight operations simultaneously (additions and multiplications). On the other hand, SFU can only perform two operations simultaneously, which are powers, logarithms, and, trigonometric functions. It means that the basic operations performed by an algorithm implemented in CUDA should be additions and multiplications. Minimizing the number of accesses between CPU and GPU: To perform operations in the GPU, it is necessary copy the information from the CPU host to the GPU device. But done this requires time. For large amount of data, the copy process generates a bottleneck. This eect can reduce considerably any gain of speed. For this reason, it is necessary to minimize the transfer of information between the CPU and the GPU. Minimizing dispersion in the execution of threads: To maximize the capabilities of the GPU, the number of operations assigned to the threads should be as similar as possible. A CUDA Kernel is executed as many times as threads are, and it cannot nish until each of all threads are completed. It means, that the execution time of a kernel will increase with the number of calculations of the largest thread. See the following example: We can perform matrix vector multiplication with a kernel. We will compute it using one thread for each row.     1 1 0 0 1 1 0 0 1 1 0 0 1 1 1 1        1 1 1 1     (2.2.1) With this choice, the rst three threads have to wait until the last one ends. The execution time of the kernel doubles for only one thread. Reduce the use of structures: The current version of CUDA supports the use of structures, but it is not fully optimized. For this reason, the use of structures should be limited as much as possible. If structures are used, they should be programmed in a simple way. In contrast with other programming languages, the use of structures of structures is highly inecient in CUDA. Avoid recursion: In the traditional programming is very common see expression like a[i] = a[i] + b[i]i= 0...N . It means that the value of the sum a[i] + b[i] will be CHAPTER 2. PARALLEL PROGRAMMING IN CUDA 16 stored in a[i] once the addition has completed. By dierence, in CUDA you have to wait until all threads nish a instruction to perform another one. Recursion expressions cannot be used in CUDA. Chapter 3 Model problem and Finite Element Formulation The objective is to implement in CUDA an algorithm to solve a system of linear equations arising from a hp-FEM method, and test it in the Laplace equation. Once done, we will compare the execution times of the program implemented in C in parallel with the implementation in CUDA. 3.1 FEM The nite element method has its origins in the theory of structures. By knowing the characteristics of individual structural elements and combining them, the governing equations for the entire structure could be obtained. This process produces a set of algebraic equations. The limitation on the number of equations that could be solved is one of the major constraints of the method. The introduction of the digital computer has made possible the solution of large systems of equations. Nowadays, Finite Element methods are used to solve dierential equations problem in many areas of science and technology, including mechanics, electromagnetism, uid mechanics and many others. Dierential equations are mathematically studied from several dierent perspectives, mostly concerned to the set of functions that satisfy the equation. Only the simplest dierential equations have analytical solutions. Moreover, most of the systems that involve dierential equations do not have a known exact solution form. When it is not possible to nd an explicit solution, it may be numerically approximated using computational techniques. This makes the coupling of dierential equations and high performance computing an incredible tool to approximate solutions in engineering problems. While the governing equations and boundary conditions can usually be written to these problems diculties, introduced by irregular geometry or other discontinuities, render makes the problems intractable analytically. 17 CHAPTER 3. MODEL PROBLEM AND FINITE ELEMENT FORMULATION 18 To obtain a solution, simplifying assumptions must be considered, reducing the problem to one that can be numerically approximated. Numerical methods provide approximate values of the unknown quantity in the region of interest (computational domain). In the FEM, this region of interest is divided into numerous connected subregions or elements within which approximate functions which are usually polynomials which are used to represent the unknown quantity, see [53] [28]. 3.1.1 Variational Formulation of Laplace Equation. To solve the problem, we will rst consider the Laplace equation with Dirichlet boundary conditions. (SP)(−∆u=f, in Ω⊂R3, u|Γ= 0, on Γ = ∂Ω, (3.1.1) Where Ω is a bounded open domain in the space R3={x= (x1, x2, x3) : xi∈R} with boundary Γ . Notice that: ∆u=∂2u ∂x2 1 +∂2u ∂x2 2 +∂2u ∂x2 3 (3.1.2) Dening bn= (n1, n2, n3) as the normal(outward)to Γ . d−→ x denotes the element of volume in R3 , and ds the surface element along Γ . In general, we nd that, ZΩ ∂v ∂xi wd−→ x+ZΩ v∂w ∂xi d−→ x=ZΓ vwbnids, i = 1,2,3. Integrating by parts in the three coordinates: ZΩ∇v·∇wd−→ x=ZΓ v∂w ∂n ds −ZΩ v∆wdx, where the normal derivative in the outward normal direction to the boundary Γ , is: ∂w ∂n =∂w ∂x1 n1+∂w ∂x2 n2+∂w ∂x3 n3. u satises 3.1.1 and the solution to the variational problem is u∈V . Since we only consider Dirichlet boundary conditions, the term associated with Neumann boundary conditions vanishes. Thus, we have: a(u, v) = (f, v)∀v∈V, where a(u, v) = ZΩ∇u·∇vd−→ x , (f, v) = ZΩ fvd−→ x . CHAPTER 3. MODEL PROBLEM AND FINITE ELEMENT FORMULATION 19 Ω can be subdivided into a set Th=C1, . . . , Cm of non-overlapping cubes Ci , Ω = m [ i=1 Ci=C1[C2[. . . [Cm, We can now dene Vh as follows, Vh={ v: v is continous on Ω, v|k is linear in K ∈Th, v|Γ= 0}, The space Vh consists of all continuous functions that are linear on each cube and vanish on Γ . Dening the basis of Vh as follows. ϕj(Ni) = δij ≡(1 if i=j 0 if i6=ji, j = 1, . . . , M. (3.1.3) Thus, the support of ϕj consist of the cubes with the common node nj . The function vh∈Vh has now the expression, vh(x) = M X j=1 ηjϕj(x), ηj=v(nj), for x∈Ω (3.1.4) It is possible to formulate the nite element method for 3.1.1 starting from the variational formulation 3.1.3. a(uh, vh) = (f, vh)∀vh∈Vh. Where the stiness matrix is a MxM matrix, whose elements are dened as: aij =a(ϕi, ϕj), and b= (bi) is a size M vector which elements are dened as bi= (f, ϕi). Stiness matrix A is usually computed by summing the contributions of the dierent cubes: a(ϕi, ϕj) = X C∈Th aC(ϕi, ϕj), where, aC(ϕi, ϕj) = ZC∇ϕi·∇ϕjd−→ x , Then, stiness matrix and load vector are dened in their discrete version as:          A=X C∈ThZC∇ϕi·∇ϕjd−→ x , b=ZΩ fϕid−→ x , (3.1.5) CHAPTER 3. MODEL PROBLEM AND FINITE ELEMENT FORMULATION 20 Where A matrix, is a sparse matrix. Each diagonal entry of the matrix represents the interaction of each basis function with itself, while o-diagonal entries correspond to interactions between dierent basis functions. Basis functions that do not share support result in a zero contribution to the stiness matrix. Thus, the resulting matrix is highly populated with zeros. 3.1.2 hp-FEM The hp-FEM (see [25], [12], [7])is a general version of the FEM. This numerical method is based on polynomials approximations that makes use of elements of variable size h and polynomial degree p . This method converges exponentially fast, when the mesh is rened using a suitable combination of h-renements and p-renements. The h-renements are performed by dividing elements into smaller ones. The p-renements are obtained by increasing the polynomial order in shape functions (see [53]). This exponential convergence (see [10]) makes the method one of the best possible choice when implementing a numerical simulation. The hp-FEM eciency relies on the capability of approximate functions with larger polynomial order, or smaller piecewise-linear elements. This capability is also extended to all the elements inside the grid. And more importantly, dierent elements may have dierent size h and and/or dierent polynomial order approximation p , which is known as hp-adaptivity. The implemented code have this capability and be able to properly simulate dierent situations where h and p can vary as requested by the user. Chapter 4 Linear equation solvers In this chapter we are going to study the basic characteristics of the most common linear equation solvers. After that we will justify the selection one of them to be implemented in CUDA. To learn more about GPU implementation of linear algebra, see [31]. The problem is expressed mathematically by this way. We want to solve the following linear system of n equations with n unknowns x1, x2, . . . , xn :          a11x1+a12x2+. . . +a1nxn=b1 a21x1+a22x2+. . . +a2nxn=b2 . . . an1x1+an2x2+. . . +annxn=bn (4.0.1) In matrix form, we have: Ax =b (4.0.2)      a11 a12 ... a1n a21 a22 ... a2n . . .. . ..... . . an1an2... ann          x1 x2 . . . xn     =     b1 b2 . . . bn      (4.0.3) For the existence of solution det(A)6= 0 . 4.1 Gaussian elimination This linear algebra procedure rst perform a forward elimination. Gaussian elimination reduces a given system to ones triangular. Second, it performs a backward elimination to solve the linear system [4]. First step : eliminate x1 in the rows:2,...,n by linear combination of rows. Second step : eliminate x2 in the rows:3,...,n by linear combination of rows. 21 CHAPTER 4. LINEAR EQUATION SOLVERS 22 Iteratively after n-1 steps      u11 u12 ... u1n 0u22 ... u2n . . .. . ..... . . 0 0 ... unn          x1 x2 . . . xn     =     k1 k2 . . . kn      (4.1.1) Gaussian Elimination algorithm: forward elimination and triangular form k= 1, . . . , n−1                            a(k+1) ij =a(k) ij i= 1, . . . , k j = 1, . . . , n a(k+1) ij = 0 i=k+ 1, . . . , n j = 1, . . . , k a(k+1) ij =a(k) ij −a(k) ik a(k) kj a(k) kk i=k+ 1, . . . , n j =k+ 1, . . . , n b(k+1) i=b(k) ii= 1, . . . , k b(k+1) i=b(k) i−a(k) ik b(k) k a(k) kk i=k+ 1, . . . , n (4.1.2) By backward elimination, it means that we star by eliminating xn in the rows:1,...,n-1 by linear combination of rows. By iteration of this process the matrix only have ones in the main diagonal, so we get the solution x. 4.2 LU decomposition We will study a direct method for solving linear systems: the LU decomposition. Given a matrix A, the aim is to build a lower triangular matrix L and an upper triangular matrix which has the following property: diagonal elements of L are unity and A=LU. [6] For the resolution of linear system : Ax=b, the system becomes LUx =b⇔Ly =b(1), Ux =y(2). (4.2.1) L=     1 l21 1 . . .. . .... ln1ln2··· lnn     U=     u11 u12 ··· u1n u22 ··· u2n .... . . unn      (4.2.2) We solve the system (1) to nd the vector y, then the system (2) to nd the vector x. The resolution is facilitated by the triangular shape of the matrices. CHAPTER 4. LINEAR EQUATION SOLVERS 23 Ly =b⇔     y1=b1/l11 yi=1 lii (bi− i−1 X j=1 lijyj)∀i= 2,3, . . . , n. (4.2.3) Ux =y⇔     xn=yn/unn xi=1 uii (yi− n X j=i+1 uijxj)∀i=n−1, n −2,...,1. (4.2.4) Before starting the system solution L and U must be found, which typically are solved by Gaussian elimination. 4.3 Cholesky decomposition Given a symmetric positive denite matrix A, the aim is to build a lower triangular matrix L which has the following property: the product of L and its transpose is equal to A. [2] A=LLT (4.3.1)   a11 a21 a31 a21 a22 a32 a31 a32 a33 =  l11 0 0 l21 l22 0 l31 l32 l33   l11 l21 l31 0l22 l32 0 0 l33  (4.3.2)   a11 a21 a31 a21 a22 a32 a31 a32 a33 =  l2 11 l21l11 l31l11 l21l11 l2 21 +l2 22 l31l21 +l32l22 l31l11 l31l21 +l32l22 l2 31 +l2 32 +l2 33   (4.3.3) For the diagonal elements (lkk) of L there is a calculation pattern: l11 =√a11 l22 =qa22 −l2 21 . . . lkk =v u u takk − k−1 X j=1 l2 kj (4.3.4) For the elements below the diagonal (lik, where i > k) there is also a CHAPTER 5. IMPLEMENTATION 30 i←0 r←b−Ax d←M−1r δnew ←rTd δo←δnew while i<imax and δnew > δoε2 do end q←Ad α←δnew dTq x←x+αd if i is divisible by iiter then r←b−Ax else r←r−αq s←M−1r δo←δnew δnew ←rTs β←δnew δo d←s+βd i←i+ 1 Algorithm 2 : PCG Algorithm 5.2 Sparse matrix formats This section explains the compression formats for sparse matrices. There is not a best method to compress sparse matrices because the eciency of the method depends on the sparsity pattern. COO The coordinate format is a simple storage scheme. The matrix information is stored in three arrays. The element jth not null is stored in jth position of the array data[j]. The other arrays contain the coordinates of the element row[j],col[j].     1−200 0 3 1 0 2 0 −5 1 0 2 0 4     (5.2.1) data[] = 1−2 3 1 2 −5 1 2 4 (5.2.2) row[] = 001122233 (5.2.3) col[] = 011202313 (5.2.4) CHAPTER 5. IMPLEMENTATION 31 It can be seen that in the row array indexes are repeated as many times as elements are not null, so we need more memory to store repeated information. However, the size of the three arrays are the same. It is an advantage in GPU computing because it produces full coalescence. CRS The compressive row storage format is a improvement of COO storage scheme. The matrix information is stored in three arrays of dierent size. Non null elements are stored sequentially in the array data[](from left to the right in rows and up to down in columns), in consequence the size of data is the number of non null elements in the matrix. The vector col[] stores the position of the element jth in the row. The third is an integer vector that stores where the row starts in vector data, where the last element is the number of total non null elements.     1−200 0 3 1 0 2 0 −5 1 0 2 0 4     (5.2.5) data[] = 1−2 3 1 2 −5 1 2 4 (5.2.6) col[] = 011202313 (5.2.7) cout[] = 02479 (5.2.8) In this format cout contains the start of the row k in kth position and the end in k+ 1th . This is an advantage when you have to iterate an instruction in rows. In contrast, rows have not got the same number of non null elements. It means that there will be dierent times of computation when we multiply a matrix by a vector. More information about sparse matrix formats can be found in [26] [40] [30] [51] 5.3 Basic operations for CG and PCG In this section we are going to explain all mathematical operations required to perform the CG and PCG method and its implementation in CUDA. All of they are well known. 5.3.1 Dot product If you have two real vectors a , b of size N. We can compute its dot product as: ~a ·~ b= N X i=1 ai·bi (5.3.1) CHAPTER 5. IMPLEMENTATION 32 A c code to perform dot product could be: float C_DotProduct( float ∗ a , float ∗ b , int N) 2 { 4 int i ; float dot ; 6 dot =0.0; for( i =0;i<N; i++) 8 dot+=(a [ i ] ∗ b[ i ] ) ; 10 return dot ; } If we think for a moment we can realise that the dot product is an operation with 2N inputs and one output. It means that at one moment the program execution must be serialized(a waste of time in CUDA). We can parallelize the pairwise sum, but the addition of the result is in serial. 1 __global__ float CUDA_DotProduct( float ∗ a , float ∗ b , float ∗ a_dot_b , int N) { 3 int i ; float dot ; 5 int row = blockIdx . x ∗ blockDim . x + threadIdx . x ; i f ( row < N && row >= 0) 7 a_dot_b[ row ] = a [ row ] ∗ b[ row ] ; 9 __syncthreads() ; dot =0.0; 11 for( i =0;i<N; i++) dot+=(a [ i ] ∗ b[ i ] ) ; 13 return dot ; 15 } To exploit parallelism of the GPU the idea is that the addition of terms is done by pairs,and this process is repeated until the result is nd. In theory it is a good idea but has some details. The dimension of the vector does not have to be a multiple of two, so there will be to keep in mind that we will have additional terms, we have not added. Also we must remember that information can be easily shared among the threads in a Block. We can use this property to perform the addition by blocks and continue the sum adding the results of each Block. There is another algorithm that calculates the dot product which is fully explained in this article. [27] In the program we have used the function cublasSdot()of the CUBLAS library, which performs the dot product highly ecient. In the ( [38]) its exposed how it works. CHAPTER 5. IMPLEMENTATION 33 Figure 5.1: Addition of terms in dot product done by pairs. Picture taken from [27] 5.3.2 Vector addition In some steps of the algorithm like r=b−Ax we need to compute a sum of vectors of size N. This operation involves three arrays allocated in the device memory. In contrast with the dot product we have 2N input values and N output values. In this case the idea of the parallelism is very simple: Computing each element as a sum of the corresponding terms in pairs. −→ a+−→ b=−−−→ a+b (5.3.2)      a1 a2 . . . an     +     b1 b2 . . . bn     =     a+b1 a+b2 . . . a+bn      (5.3.3) We must note that despite its simplicity not all sums are performed simultaneously. So before using the result of the sum in another operation, we must be sure that all the terms have been calculated. It can be done by using the command syncthread() . 1 void C_Vector_Addition( float ∗ a , float ∗ b , float ∗ a_add_b, int N) { 3 int i ; 5 for( i =0;i<N; i++) a_add_b[ i ]=a [ i ]+b[ i ] ; 7 } CHAPTER 5. IMPLEMENTATION 34 1 __global__ void CUDA_Vector_Addition( float ∗ a , float ∗ b , float ∗ a_add_b, int N) { 3 int row = blockIdx . x ∗ blockDim . x + threadIdx . x ; i f ( row < N && row >= 0) 5 a_add_b[ row ] = a [ row] ∗ b[ row ] ; 7 __syncthreads() ; } This code can't be used to solve equations like: a=a+b (typically done in c)because you can't have accessed at positions to write until all operations are done. Change this way of thing is one of the most typical mistakes when programming in CUDA. [14] 5.3.3 Matrix vector multiply in CRS Although the matrix-vector operation is well known, it should be noted that the matrix is stored in CRS format and we have to choose properly the instructions to perform for each thread of code. For a matrix of dimension N this operation implies N2+N input data and N output data.      a11 a12 ... a1N a21 a22 ... a2N . . .. . ..... . . aN1aN2... aNN          x1 x2 . . . xN     =    PN i=1 a1ixi PN i=1 a2ixi . . . PN i=1 aNixi      (5.3.4) Each term of the multiply vector is a sum of products of the same row.     =⇒=⇒=⇒=⇒ =⇒=⇒=⇒=⇒ =⇒=⇒=⇒=⇒ =⇒=⇒=⇒=⇒        =⇒ =⇒ =⇒ =⇒    =    =⇒ =⇒ =⇒ =⇒     (5.3.5) We can compute each element with one thread. The kernel to compute matrix vector operation with a matrix storage in a two dimensional array. __global__ void Matrix_Vector_GPU( float ∗ ∗ A, // Matrix 2 float ∗ vector_in , // Input vector float ∗ vector_out , // Output vector 4 int N) //Dimension { 6 int i ; // iteration index int row = blockIdx . x ∗ blockDim . x + threadIdx . x ; 8 i f ( row < N && row >= 0) for( j =0; j < N; j++) 10 vector_out [ row]+= A[ row ] [ j ] ∗ vector_in [ j ] CHAPTER 5. IMPLEMENTATION 35 12 __syncthreads() ; } The kernel to compute matrix vector operation with a matrix storage in CRS format. 1 __global__ void Matrix_Vector_CRS_GPU( float ∗ data , // Matrix in CRS format int ∗ col_position , 3 int ∗ cout , float ∗ vector_in , // Input vector 5 float ∗ vector_out // Output vector int N) //Dimension 7 {int i ; //index 9 int start , end ; // iteration limits int row = blockIdx . x ∗ blockDim . x + threadIdx . x ; 11 float sum; i f ( row < DIM && row >= 0) 13 { sum = 0.0; 15 start = cout [ row ] ; // cout [ k ] t e l l us where the k row begins end = cout [ row+1];//cout [ k+1] t e l l us where the k row ends 17 for( i = start ; i < end ; i++) 19 sum+= data [ i ] ∗ vector_in [ col_position [ i ] ] ; 21 vector_out [ row]=sum; } 23 __syncthreads() ; } 5.3.4 SAXPY (Single-precision real Alpha X Plus Y) At each step of the algorithm we have to recalculate x, r and p. All these calculations are of the form: x(m)←− x(m−1) +α(m)p(m−1) (5.3.6) r(m)←− r(m−1) −α(m)Ap(m−1) (5.3.7) p(m)←− r(m)+β(m)p(m−1) (5.3.8) All of them depend on the same in the previous step. These operations are basically additions, and the kernel required to perform this operation is simple. __global__ void Add_Vectors( float ∗ a , // f i r s t input vector 2 float ∗ b , // second input vector float ∗ saxpy , // addition vector 4 float sign1 , CHAPTER 5. IMPLEMENTATION 36 float sign2 , 6 int N) { 8 int row = blockIdx . x ∗ blockDim . x + threadIdx . x ; i f ( row < N && row >= 0) 10 {saxpy [ row ] = sign1 ∗ a [ row ] + sign2 ∗ b[ row ] ; 12 } __syncthreads() ; 14 a [ row]=saxpy [ row ] ; 16 } Despite its simplicity is kind of operations require an intermediate vector to store the result and copy it to the original vector.The syncthreads command is essential because the elements of the a vector will not be overwritten until all them are calculated. One of the premises in GPU computing is to reduce the stored memory. To take full advantage of the GPU capabilities is necessary to consider these details. Nvidia provides libraries in which these functions are optimized [38]. SAXPY (Single-precision real Alpha X Plus Y) multiplies the vector x by the scalar α and adds it to the vector y overwriting the latest vector with the result. Hence, the performed operation is y[j] = αx[j] + y[j] . [13] 5.4 Structure of the code The code has been structured as follows: 1.- To generate a sparse symmetric positive denite matrix A. It's the kind of matrices generated by hp-FEM. A is the stiness matrix. 2.- To generate a solution vector −→ xsol . 3.- To generate a vector −→ b with the matrix-vector product A−→ xsol =−→ b 4.- Using an algorithm implemented in CUDA the system of equations A−→ x=−→ b is solved. The execution time tcuda is measured. 5.- Using an algorithm implemented in C parallel the system of equations A−→ x=−→ b is solved. The execution time tparallel is measured. 6.- We compare execution times tcuda tparallel The vector −→ xsol allows us to check the obtained solutions. Chapter 6 Results The results exposed in this chapter was obtained with a version of the program that works with static memory. This fact, only allows to arrive to dimension 40000 (approximately the size of the RAM memory of the CPU in double precision). This limitation is a problem to show the advantages of GPU computing. In this order of size the CUDA implementation of the solver works between 2-3 times faster in GPU than CPU. For lower size (the order of thousands) CUDA does not improves CPU implementation, the GPU programs spent the most of time copying within CPU and GPU. For higher dimension, the implementation of the program with dynamic memory provokes segmentation faults in the execution. With dynamic memory sizes of millions can be easily reach. Another fact relevant is that the hp-FEM in Laplace problem generates an structure of matrix that is usually solved in a reduce number of iterations (for moderate tolerances). For dimension 40000, the program nish in less than 4000 iterations. 37 Chapter 7 Conclusions and future work The CUDA implementation of a iterative solver versus CPU implementation only exhibits improvements at higher dimension. At lower dimensions it is most practical CPU implementation. At this dimension, all the speed improvements in the basic operations are neglected by time of CPU to GPU copying process. In our future work we want to solve the problem with dynamic memory implementation. The solve of this problem will provide us a very powerful to solve systems of linear equations. Also, we want to be able to implement a preconditioner. It will allow us to face higher dimension (tens of millions). 38 Bibliography [1] Amhdal's law, wikipedia. http://en.wikipedia.org/wiki/Amdahl's_law. [2] Cholesky decomposition, math-linux. http://www.mathlinux.com/spip.php?article43. [3] Conjugate gradient method, math-linux. http://www.mathlinux.com/spip.php?article54. [4] Gaussian elimination, math-linux. http://www.mathlinux.com/spip.php?article53. [5] General-purpose computation on graphics hardware. http://gpgpu.org/. [6] Lu decomposition, math-linux. http://www.mathlinux.com/spip.php?article51. [7] M. Ainsworth and J. Oden , A procedure for a posteriori error estimation for hp nite element methods , Computer Methods in Applied Mechanics and Engineering, 101 (1992), pp. 7396. Elsevier. [8] S. Akhter and J. Roberts , Multi-Core Programming , vol. 33, 2006. Intel Press. [9] G. Almasi and A. Gottlieb , Highly parallel computing , (1988). Menlo Park, CA (USA); Benjamin-Cummings Pub. Co. [10] I. Babu²ka and B. Guo , The h, p and h-p version of the nite element method: basis theory and applications, advances in engineering software . [11] M. Baskaran and R. Bordawekar , Optimizing sparse matrix-vector multiplication on gpus , IBM Research Report, (2008). [12] C. Baumann and J. Oden , A discontinuous hp nite element method for the euler and navierstokes equations , International Journal for Numerical Methods in Fluids, 31 (1999), pp. 7995. Wiley Online Library. [13] N. Bell and M. Garland , Ecient sparse matrix-vector multiplication on cuda , NVIDIA Corporation, NVIDIA Technical Report NVR- 2008-004, (2008). 39