Repositorio Institucional de Documentos

CUDA implementation of integration rules within an hp-Finite Element code. Author: Adrián García 11 , Directors: David Pardo 2* , and Ricardo Celorrio 3+ * Departamento de Matemática Aplicada, UPV/EHU University + Departamento de Matemática Aplicada, Zaragoza's University. September 4, 2012 1 [email protected] 2 [email protected] 3 [email protected] 2 Contents 1 Introduction. 7 1.1 High Performance Computing. And Parallelization. . . . . . . 8 1.1.1 Amhdal'slaw........................ 8 1.1.2 Parallel computing. . . . . . . . . . . . . . . . . . . . . 8 1.2 WhyCUDA. ........................... 10 1.2.1 GPUEciency....................... 11 1.2.2 Architecture & Programming Model. . . . . . . . . . . 14 1.3 Finite Element Problem. . . . . . . . . . . . . . . . . . . . . . 16 1.3.1 Poisson's Equation Variational Formulation. . . . . . . 17 1.3.2 hp-FEM. ......................... 19 2 Problem and Proposed Solutions. 21 2.1 Transformed Space. . . . . . . . . . . . . . . . . . . . . . . . . 21 2.2 Shape Functions, Lagrange's Polynomials. . . . . . . . . . . . 22 2.3 Gaussian Quadrature, integration method. . . . . . . . . . . . 23 3 Code Implementation. 25 3.1 Lagrange's polynomials. . . . . . . . . . . . . . . . . . . . . . 25 3.2 Integration. ............................ 26 3.2.1 b Integration........................ 26 3.2.2 A Integration. ...................... 28 3.3 System integration. . . . . . . . . . . . . . . . . . . . . . . . . 32 3.4 MatrixAssembly.......................... 32 4 Results. 35 5 Conclusions And Future Work. 37 3 4 CONTENTS Abstract With the introduction in 2006 of CUDA architecture for Nvidia GPUs, a general purpose parallel computing architecture. A new programming model where scientist and engineers had found and answer to their claim of high performance computation with moderate cost. This thesis work begin with the motivations that underline parallel computation and high performance computation. Later it is discussed why CUDA architecture and General Porpoise Graphics Processor Units (GPGPU) had modied the high performance computation world and why it is turning into GPGPU computing. The author has chosen this new language because of its multiple and countless possibilities focusing in engineering applications. CUDA's novelty it is both a motivation and an obstacle because to the lack of support and documentation specically to engineering applications. This work has two main goals. First one is to implement a FEM integration code in a GPU architecture using CUDA language and a Nvidia graphic device. This code performance will be compared to a parallelized CPU code developed by the same author. Second, once the integration has been done the obtained solution is transform into a CRS compress matrix which is the beginning of another work, solve the system described by that particular matrix. The performance of matrix assembly in GPU it is also measured and compared to its respective code in parallelized CPU. 5 6 CONTENTS Chapter 1 Introduction. Since the early days of computers, computer graphics had been essential to make possible a better human-machine interaction. While some applications only need 2D graphics many other start using 3D graphics. This applications whose most representative are computer games, photography& video edition, have evolved in applications who demand a huge computation capability. As this applications require more and more resources standard CPUs were not enough to support their function. That was the main reason to build up a device specically dedicated to graphics, General Porpoise Graphics Processor Units (GPGPU / GPU) were born. GPU microprocessor were born to draw points inside a dened screen what is basically linear algebra operations. Within the years GPU advanced and 3D computer graphics became an integral part of many computer applications. The 3D to every application resulted in a huge market of aordable 3D graphics cards. Devices prices get lower and lower and its compute capability higher and higher. Few years ago GPUs outperformed the number of oat points operations per second that CPUs were capable, this was the trigger to a new era. GPUs could now be used to high performance applications within a moderate cost. That was particularly benecial to science and engineering world who is always looking for more compute capabilities. Science and engineering world turned into this kind of computation obtaining incredible results to many problems and realizing in some other areas that making use of large number of cores with less compute capability it is not always a panacea. The main problem of GPU computing it is closely related to Amhdal's law (section 1.1.1). This computational law says. The speedup of a program using multiple processors in parallel computing is limited by the sequential fraction of the program. In this articles two master students will implement a Finite Element code, which is a solution largely utilized in engineering to simulate structures. And compare the results of running this program in both CPU and GPU. 7 8 CHAPTER 1. INTRODUCTION. 1.1 High Performance Computing. And Parallelization. It was 2006 when Intel launched its Core Duo into the market. This was the beginning of a revolution, many years of investigation ended with a single CPU with two cores. Single processor units or mono-cores rapidly died because of the signicant computational advantages of utilizing several cores. CPUs had continued evolving and nowadays are available even to commodity computers to home utilization counting on several cores. But is this growth in the number of cores inside a single processor the solution to every problem in high computation problems? The answer is no. If we take a look to Moore 's Law, the number of transistor inside a microprocessor doubles every 18 months. Does it mean that microprocessors duplicates its velocity every 18 months? One more time the answer is no. CPU's velocity relies on internal clock frequency inside each microprocessor. Anyway it is a fact that with more microprocessor the program performance increases. This section will make a brief introduction to Amhdal 's law. And after introducing Amhdal 's law we will have another introduction to the techniques used in CPU parallel computing. 1.1.1 Amhdal's law. It was 1967 when Gene Ahmdal proposed the following observation. The speedup of a program using multiple processors in parallel computing is limited by the time needed for the sequential fraction of the program. Amdahl's law is a model for the relationship between the expected speedup of parallelized implementations of an algorithm relative to the serial algorithm, under the assumption that the problem size remains equal when parallelized. In case a part of the algorithm can not be parallelized, the program speedup will improve only because of the part that is possible to run in parallel. And the improvement will depend in how much of that algorithm can be parallelized. More technically, the law is concerned with the speedup achievable from an improvement to a computation that aects a proportion P of that computation where the improvement has a speedup of S. 1 (1 −P) + P S (1.1.1) 1.1.2 Parallel computing. Parallel computing is a form of computing where many cores execute the same program at a once. It has been widely used long time ago but it was not until 2006 that processors with multiple cores were accessible to non 1.1. HIGH PERFORMANCE COMPUTING. AND PARALLELIZATION. 9 Figure 1.1: The speedup of a program using multiple processors in parallel computing is limited by the sequential fraction of the program. For example, if 95% of the program can be parallelized, the theoretical maximum speedup using a parallel algorimth would be 20x as shown in the diagram. professional users. Before 2006 parallel programs had to run in physically separated computers, this still is a technique widely used nowadays. Carry out the computation simultaneously turns into an improvement of required time to execute the program. First approximation to parallel computers can be divided in two groups. Processes that run parallel inside a single machine and processes that require several computers to perform the task, to this group belong computer cluster, grids and MMP's. Closer this article main goal are multi-core, multi-processor machines and symmetric multiprocessing techniques. Who are inside the category of parallel execution using only one computer. Multi-core processor. This are nowadays processors, they count with several cores, inside a single silicon component with two or more cores. The improvement in performance gained by the use of a multi-core processor depends on the software algorithms used and their implementation. Possible gains are limited by the fraction of the algorithm that can be run in parallel. Best case, so-called embarrassingly parallel problems may realize speedup factors near the number of cores, or even more if the problem is split up enough to t within each core's cache, avoiding the use of main system memory which is much slower. But this is not usual neither real, this problems exist mostly as laboratories examples. 16 CHAPTER 1. INTRODUCTION. Coalesced Memory Access. When the kernels are been executed Blocks and Threads have access to device global memory (gure 1.4). Memory access it has always being important both in CPU and in GPU to ensure ecient code. In GPU programming case this is a particular problem, to be handled with care. A coalesced memory transaction is one in which all of the threads in a half-warp access global memory at the same time. This is oversimple, but the correct way to do it is just have consecutive threads access consecutive memory addresses. So, if threads 0, 1, 2, and 3 read global memory 0x0, 0x4, 0x8, and 0xc, it should be a coalesced read. In a matrix example, keep in mind that you want your matrix to reside linearly in memory. You can do this however you want, and your memory access should reect how your matrix is laid out. So, the 3x4 matrix below. 0123 4567 89ab could be done row after row, like this, so that (r, c) maps to memory (r·4 + c) . 0123456789ab Suppose you need to access element once, and say you have four threads. Which threads will be used for which element? Probably either. thread 0: 0, 1, 2 thread 1: 3, 4, 5 thread 2: 6, 7, 8 thread 3: 9, a, b or thread 0: 0, 4, 8 thread 1: 1, 5, 9 thread 2: 2, 6, a thread 3: 3, 7, b Which is better? Which will result in coalesced reads, and which will not? Either way, each thread makes three accesses. Let's look at the rst access and see if the threads access memory consecutively. In the rst option, the rst access is 0, 3, 6, 9. Not consecutive, not coalesced. The second option, it's 0, 1, 2, 3. Thus is consecutive and therefore coalesced . 1.3 Finite Element Problem. The physical concept on which the nite element method is based has its origins in the theory of structures. The idea of building up a structure by tting together a number of structural elements was used in the early truss 1.3. FINITE ELEMENT PROBLEM. 17 and framework analysis approaches employed in the design of bridges and buildings in the early 1900s. 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 simultaneous algebraic equations. The limitation on the number of equations that could be solved posed a severe restriction on the analysis. The introduction of the digital computer has made possible the solution of the large-order systems of equations. Finite Element methods are based on dierential equations to solve problems. Dierential equations arise in many areas of science and technology, areas so disperse as classical mechanics, electromagnetism, uid mechanics and basically any science and technology area. Dierential equations strength rely on their capability to link together a physical variable and its variation. This simple concept makes this mathematical tool one of the most essential mathematical knowledge to every scientist and engineer. 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 explicit solution formulas. Moreover most of the systems that involve dierential equations study do not have an exact solution form. When it is not possible to nd an explicit solution it may be numerically approximated using computing techniques. The theory of dynamical systems puts emphasis on qualitative analysis of systems described by dierential equations, while many numerical methods have been developed to determine solutions with a given degree of accuracy. 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 but diculties introduced by irregular geometry or other discontinuities render the problems intractable analytically. To obtain a solution simplifying assumptions must be performed, reducing the problem to one that can be solved, or a numerically approximated. Numerical methods provide approximate values of the unknown quantity only at discrete points in the region. In the nite element method, the region of interest is divided up into numerous connected subregions or elements within which approximate functions which are usually polynomials and are used to represent the unknown quantity. 1.3.1 Poisson's Equation Variational Formulation. (WP)(−∆u=f, x ∈ΩR2 u|Γ= 0 (1.3.1) To solve the problem written in weak formulation to Poisson's equation 18 CHAPTER 1. INTRODUCTION. it is necessary to nd a function which Laplacian equals, ∆u=∂2u ∂x2 1 +∂2u ∂x2 1 (1.3.2) Using Green's formula. ZΩ ∆u·vdx =−ZΩ ∇u· ∇vdx +ZΓ v∂u ∂ndγ (1.3.3) Where the gradient is equal ∆u=∂2u ∂x2 1 +∂2u ∂x2 2 =ux1x1+ux2x2 ∇u=∂2u ∂x2 1 ,∂2u ∂x2 2= (ux1, ux2) (1.3.4) The product of the gradients ∇u· ∇n= (ux1, ux2)(vx1, vx2) = ux1, vx1+ux2, vx2 (1.3.5) And the derivative respect to −→ n . ∂u ∂n =∇u·n=ux1·n1+ux2·n2 (1.3.6) to nd a variational formulation, rst a trial space is needed. (V P)(v∈ΩR l0(Ω) (1.3.7) If we take a function and apply scalar product −ZΩ (∆u)·vdx =ZΩ f·vdx;∀v∈V Green's function ZΩ ∇u∇vdx −ZΓ v∂u ∂ndγ = ZΩ ∇u∇vdx −0 = (1.3.8) This trial space must satisfy continuity and boundary conditions vΓ= 0 . It is known that if a solution exist is unique. Making the problem discontinuous, is equal as it is done in one dimension. A trial space Vh is needed, and this space has to satisfy. Vh⊂V s.t. dim(Vh)<∞ (1.3.9) This space will have a base ϕ1, ϕ2, . . . , ϕM in Vh . As proved before the system could be proposed in two formulations, and is easier to solve in variational formulation. The steps needed to solve the problem, are. 1.3. FINITE ELEMENT PROBLEM. 19 1. Find a suitable base in the trial space. un= M X j=1 ujϕj(x) (1.3.10) 2. Find the stiness matrix and load tensor . (V Ph)       un∈R a  M X j=1 uj∇ϕj(x),∇ϕi(x) =L(ϕi(x)) i= 1, . . . , M (1.3.11) And operations needed to perform this task are. A = (ϕj(x), ϕi(x)) ZΩ ∇ϕj· ∇ϕi b =L(ϕi(x)) ZΩ f·ϕidx (1.3.12) 1.3.2 hp-FEM. hp-FEM is a general version of the nite element method FEM . This numerical method is based on polynomials approximations that makes use of elements of variable size h and polynomial degree p . The origins of hp-FEM date back to the pioneering work of Ivo Babuska who discovered that the nite element method converges exponentially fast when the mesh is rened using a suitable combination of h-renements which is make by dividing elements into smaller ones and p-renements increasing polynomial order in shape functions. This exponential convergence makes the method one of the best possible choice when implementing a numerical simulation. 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 what is more important dierent elements may have dierent size h or dierent polynomial order approximation p and that is known as hp-adaptivity. hp-adaptivity as a combination of h-adaptivity splitting elements in space while keeping their polynomial degree xed and p-adaptivity increasing their polynomial degree. Related to work's problem, the code implemented must have this capability and be able to properly simulate dierent situations where h and p can vary as the user likes. Up to the rst version polynomial approximation is fully implemented up to 9th order polynomials. And h implementation forces h to be equal to the three directions in space, making this way an homogeneous grid. 20 CHAPTER 1. INTRODUCTION. Chapter 2 Problem and Proposed Solutions. The main goal of this article is to measure the ratio of execution times of two calculus algorithms for the same problem. One of them is a C code with CPU parallelelized implementation. The second code makes the same calculations inside a General Purpose Graphics Processor Unit (GPGPU) . Both codes runs under the same execution, as will be explained later in this work. The code simulates a time static and homogeneous three dimensional grid. Distance between elements h and polynomial order approximation p is xed at compilation time as well as the number of elements per side which will give the total number of elements. This allows to simulate dierent grid sizes and polynomial approximations to compare executions times. Solutions to the main diculties arisen when creating the algorithms will be now briey discussed. 2.1 Transformed Space. As seen in section 1.3.1. It is easier to solve the problem inside a transformed space. The space chosen is an homogeneous three dimensional cube in coordinates [−1 : 1],[−1 : 1],[−1 : 1] . This space has the advantage to perfectly t the selected integration method which is Gauss Quadrature dened between [−1 : 1] . Gauss quadrature can be utilised in [a:b] spaces, transforming that space. So the advantage to directly transform, [xa:xb][ya:yb][za:zb]→[−1 : 1][−1 : 1][−1 : 1] (2.1.1) is to make only one space transformation. To transform space is easy in this particular case where the grid is cubic and homogeneous. Let xc, yc, zc be the coordinates in the center of the element inside the weak formulation space. 21 22 CHAPTER 2. PROBLEM AND PROPOSED SOLUTIONS. Then the coordinates in variational formulation space to a point x, y, z , will be. ξ=x−xc hη=y−yc hζ=z−zc h (2.1.2) Where ξ, η, ζ dene the coordinates inside transformed space. 2.2 Shape Functions, Lagrange's Polynomials. Shape functions choice it is not trivial. Shape functions must have value ϕij = 1 when i=j and have value ϕij = 0 in any other case. Furthermore it is mandatory to take into account the calculus succession necessary to obtain a correct result. The program implements Lagrange's polynomials as shape functions. Given a set of points xi, yi Lagrange polynomial is the polynomial of the least degree that at each point xi assumes value yi . Given (x0, y0),...,(xj, yj),...,(xk, yk) (2.2.1) where no two xi are the same, the interpolation polynomial in the Lagrange form is a linear combination. `j(x) := Y 0≤m≤k m6=j x−xm xj−xm =(x−x0) (xj−x0)· · · (x−xj−1) (xj−xj−1) (x−xj+1) (xj−xj+1)· · · (x−xk) (xj−xk) (2.2.2) Given the initial assumption that no xi are equals, every xi−xj6= 0 , so the expression it always properly dened. As requested in equation ?? , Lagrange's polynomials satisfy that particular condition. For all i6=j , `j(x) includes the term (x−xi) so the numerator will be zero at x=xi `j6=i(xi) = Y m6=j xi−xm xj−xm =(xi−x0) (xj−x0)· · · (xi−xi) (xj−xi)· · · (xi−xk) (xj−xk)= 0 (2.2.3) On the other and, if i=j . `i(xi) := Y m6=i xi−xm xi−xm = 1 (2.2.4) There is one more main reason to use Lagrange's polynomials. It is related to high performance computing. Shape functions must be hierarchical. When to obtain the n−esime value of a polynomial is fn(x) = fn−1(x)·x . 2.3. GAUSSIAN QUADRATURE, INTEGRATION METHOD. 23 As an example let assume that we need to evaluate function f(x) . Then f(x) is a hierarchical function if can be evaluated like. f1(x)=1 f2(x) = x f3(x) = x2 f4(x) = x3 (2.2.5) And continue evaluating until the polynomial ends. As can be seen any polynomial like, f(x) = a+bx +cx2+. . . +nxn+1 can be hierarchically computed. And there is a huge mistake on trying that approach in high performance computing, and at any computing problem in general. The fact is when approaching a polynomial result with that method x it is very likely to reach extreme values, zero or innite. Inside our particular case. The program uses a transformed space between [−1 : 1] , and up to 9th order polynomials so succession xn will reach zero value returning a wrong result. By using Lagrange's polynomials this issue does not appear because there is no more potency to be evaluated. Instead subtractions are evaluated and multiplied and x value is no longer modied. On the other and Lagrange's polynomials introduce a diculty when its gradient is calculated. Because it is necessary to introduce an if condition which is not advisable at all, but it is a minor issue when compared to the evaluation point reaching to zero. This is an easy and systematic method of generating shape functions of any order now can be achieved by simple products of Lagrange polynomials in the two or more coordinates. Thus, in three dimensions, if we label the node by its column and row number, I, J and K we have. Na≡NIJK =ln I(ξ)lm J(η)lp K(ζ) (2.2.6) where n, m and p stand for the number of subdivisions in each direction. 2.3 Gaussian Quadrature, integration method. A quadrature rule is an approximation to a dened integral. To perform the integration task the program uses Gaussian Quadrature method. This method evaluates a weighted sum of the function evaluated in certain points. More specically if the evaluating function is a polynomial, Gaussian quadrature reach an exact solution within N= 5 evaluation points. Gaussian quadrature to one dimension. Z1 −1 f(x)dx ≈ N X i=1 wif(xi) (2.3.1) 24 CHAPTER 2. PROBLEM AND PROPOSED SOLUTIONS. As the program runs inside a three dimensional space, quadrature has to approximate a space in three dimensions. Z1 −1 f(x)dx Z1 −1 f(y)dy Z1 −1 f(z)dz ≈ N X i=1 N X j=1 N X k=1 wiwjwkf(xi)f(yj)f(zk) (2.3.2) This integration method oers an exact solution to a polynomial when evaluated up to ve points. So the program implements a ve-point Gauss quadrature function to obtain the best possible solution. It should be noted the dierences when implementing the calculus of elements aij and bi . If we look equation 1.3.12 the dierences between both implementations will be noticed. To implement bi algorithm shape function is evaluated itself. But when performing aij calculus shape function is not evaluated. It is the two shape's functions gradient product to be evaluated. And as said earlier that particularity introduces an if condition which is not desirable but inevitable. Despite the gradient calculus inconvenient Gaussian quadrature method oers a good eciency to high performance computing and it is remarkable that Gaussian quadrature is an exact solution when evaluating polynomials. Chapter 3 Code Implementation. In chapter 2 several problems were discussed. This chapter the solutions used while implementing the code. Sections are dedicated to each of the main issues that have been found while working in this program. 3.1 Lagrange's polynomials. Lagrange's polynomial theory was exposed in 2.2. To implement Lagrange's polynomials it is not a dicult task. As the grid point is known, it is direct to transform a given point and its nearest neighbours to transformed space. If we denote as k point to have value 1 , to implement Lagrange's polynomials inside a three dimensional space three more points are needed. This points 25 32 CHAPTER 3. CODE IMPLEMENTATION. 3.3 System integration. Once we have build the dierent algorithms needed to solve the system it is time to develop a piece of code that uses those algorithms to solve the problem proposed in equation 1.3.12. It is important to remark that A matrix is made by boxes. Each row is formed by several boxes which represents the interaction between two element. Inside a row we will nd the box that represents the interaction with the element itself, and the boxes who represents the element interaction with its nearest neighbours. for idx ←0 to Number of Elements do tSpace ←Get Transformed space associated to idx ; for li ←0 to 4 do lpol ←Get Lagrange0s polynomials ; end b←Integrate tSpace0with lpol0 ; A0←Integrate tSpace0with lpol0 ; for li ←1 to 4 do Ali ←Integrate tSpaceli lpolli ; end end Algorithm 5 : Calculate A matrix and b vector. A matrix is made of boxes 3.4 Matrix Assembly. Before facing system solving, it is necessary to assembly A matrix. Already we have the values of the integrations but they are not properly distributed. As mentioned before A is formed by boxes, this step turns this box-formed matrix into a Compressed Storage Row (CRS) matrix. The main problem to perform this task is to handle memory directions properly. There are several ways to do it and in this work I have chosen a solution that implies look for every value twice and write it once. This solutions has been taken because 3.4. MATRIX ASSEMBLY. 33 reading memory is much faster than writing. for idx ←0 to Number of Elements do Matrix is copied row by row for p←0 to p−order do for i←0 to Boxes Behind the diagonal do for j←0 to p−order do CRS[cont]←A[(idx −(i+ 1)) ∗p−order2∗4) + (i∗ p−order2)+(p∗p−order) + j ; cont + + ; end end for j←0 to p−order do CRS[cont]←A[idx ∗p−order2∗4) + (p∗p−order) + j ; cont + + ; end for li ←1 to Boxes after the diagonal do for j←0 to p−order do CRS[cont]← A[idx∗p−order2∗4)+(i∗p−order2)+(p∗p−order)+j ; cont + + ; end end end end Algorithm 6 : A matrix assembly into Compressed Row Storage format. 34 CHAPTER 3. CODE IMPLEMENTATION. Chapter 4 Results. Before expose and discuss obtained results is important to describe the hardware equipment that is going to be used. The program will run in a home computer with a CPU AMD Phenom(tm) II X6 1090T Processor which has a top frequency of 3.3Mhz . Mother board has 8GB as RAM memory. And nally the GPU which is an GeForce GTX 550 Ti Nvidia graphics card, all of this running under Linux Ubuntu 12.04 . This equipment presents two direct problem, given the fact that CPU is much better than GPU. And the fact that this GPU can only operate in oat which is single precision with no have value in engineering world. Anyway this work is about compare the rate between CPU and GPU executions. Previous chapters have widely talked about GPU its structure and its programming model. So the results are presented directly. Figure 4 represents the time required to execute integration process. Figure 4 represents the time taken to transform A from boxess to CRS format. Lastly Figure 4 represents the ration between the executions time. Every execution has been made with h= 5 and it is the number of elements and the shape function polynomial approximation the variables that are handled. This may not appear a good solution, but after many executions parameter h was almost irrelevant when measuring time execution. Polynomial order approximation has chosen to be 3 and 6 . Higher polynomials order make the program to fully occupy device and cause problems to the OS graphics environment which is running inside the same GPU so no data could be acquired to higher polynomial approximation order in this precise equipment. 35 36 CHAPTER 4. RESULTS. Figure 4.1: This chart represents executions times to dierent grid sizes. As we can see the GPU implementation it is not as god as it should. And it runs much slower that parallel CPU. Figure 4.2: Once more CPU parallel performance outcomes when task involves reading and writing. In this particular case the incredible diculty of making A vectors coalesced severely penalizes the GPU. Figure 4.3: This particular chart express the time ratios between the dierent executions that have been made. Chapter 5 Conclusions And Future Work. This work represents rst step towards the implementation of a fully operational hp-FEM solver. It is true that in this work results does not support the idea of a functional hp-FEM solver but despite the poor results exist are a number of factors that worth to analyse. CUDA is a new programming language and a new programming structure this translates into a huge lack of support. There are few specialized books and even few of those books are oriented to high performance computation. Most of those CUDA specialized books address issues related to graphic user interface and image processing improvement not engineering problems. In conjunction with this rst problem author's CUDA inexperience surely has derived into a poorly optimised implementation. As it is easy to see in gure 4 vectors are clearly not coalescence which is explained in section 1.2.2. When it is referenced to integration problem it is Amhdal's law who has tricked the results. Implemented algorithm is good enough to have parallelism in 6 cores as CPU has. Problem appear when the number of cores grown but the velocity of themselves descends. And with that number of cores a completely new algorithm is needed. The nal conclusion about CUDA is the long road of development that still has this new programming language and the bright future is GPU computing. Although in this particular job performance was not achieved expected there countless studies have shown that the feasibility of code implentación graphics cards in order to make high performance computing. 37 38 CHAPTER 5. CONCLUSIONS AND FUTURE WORK. Acknowledgment. I want to thank and dedicate this work to all the people who have been close to me during the making of it. To my parents, to Rocio, Raúl and Javi. BIFI's people specially sysadmins that gave full support and access to their computers. I'm specially grateful to Asier Lacasta without your help it would have been impossible to do this job. Finally I want to thank David Pardo and Ricardo Celorrio by guide me with this dicult work. 39 40 CHAPTER 5. CONCLUSIONS AND FUTURE WORK. Bibliography [1] Jens Krüger and Rüdiger Westermann, Linear Algebra Operators for GPU Implementation of Numerical Algorithms. Computer Graphics and Visualization Group, Technical University Munich, 10.1.1.1.3310. [2] K. Fatahalian, J. Sugerman, and P. Hanrahan, Understanding the Eciency of GPU Algorithms for Matrix-Matrix Multiplication. Graphics Hardware (2004), Stanford University 10.1.1.1.6823. [3] Stephanie Winner *, Mike Kelley *** , Brent Pease **, Bill Rivard*, and Alex Yen *** Hardware Accelerated Rendering Of Antialiasing Using A Modied A- buer Algorithm. * 3Dfx Interactive, San Jose, CA USA, ** Bungie West, San Jose, CA USA, *** Silicon Graphics Computer Systems, Mountain View, CA USA 10.1.1.46.6965 [4] Andreas Schilling, A New Simple and Ecient Antialiasing with Subpixel Masks. Computer Graphics, vol. 25, no. 4, July 1991 (SIGGRAPH '91 Proceedings), pp. 133141 10.1.1.59.3971 [5] Je Bolz, Ian Farmer, Eitan Grinspun, Peter Schröder Sparse Matrix Solvers on the GPU: Conjugate Gradients and Multigrid. Caltech 10.1.1.112.5723 [6] Tor Dokken Trond, R. Hagen Jon, M. Hjelmervik The GPU as a high performance computational resource. SINTEF ICT, Applied Mathematics P.O. Box 124 Blindern 0314 Oslo, Norway 10.1.1.133.5648 41