Masaryk University Faculty of Informatics Bachelor Thesis Utilization of OpenCL Language for Acceleration of Scientific Calculations Richard Sirný 2010 Bibliographic Identification Bachelor thesis title: Utilization of OpenCL Language for Acceleration of Scientific Calculations Author full name: Richard Sirný Study program: Applied Informatics Study field: Applied Informatics Supervisor: Mgr. Jan Kmuníček, Ph.D. Supervisor specialist: RNDr. Petr Kulhánek, PhD. Year of defend: 2010 Keywords: electrostatic energy calculation, electrostatic interactions, GPGPU, OpenCL, parallel computing Bibliografická identifikace Téma bakalářské práce: Využití jazyka OpenCL k akceleraci vědecko-technických výpočtů Jméno autora: Richard Sirný Studijní program: Aplikovaná informatika Studijní obor: Aplikovaná informatika Školitel: Mgr. Jan Kmuníček, Ph.D. Konzultant: RNDr. Petr Kulhánek, Ph.D. Rok obhajoby: 2010 Klíčová slova: elektrostatické interakce, GPGPU, OpenCL, paralelní výpočty, výpočet elektrostatického potenciálu ii Declaration Hereby I declare, that this thesis is my original authorial work, which I have worked out by my own. All sources, references and literature used or excerpted during elaboration of this work are properly cited and listed in Bibliography section. iii Acknowledgment I would like to thank to my family, friends and teachers. Very special thanks goes to Jan Kmuníček and Petr Kulhánek whose support, leadership and helpful attitude helped me greatly when writing this work. iv Abstract Computational chemistry problems can be computationally exhaustive. Introduction of parallelism into these problems can lead to performance boost. The aim of this thesis is to assess the potential of OpenCL programing language for parallelization of scientific calculations. The applicability of OpenCL is demonstrated on calculation of electrostatic interactions in molecular systems. The theoretical part of the thesis provides introduction to parallel computing and description of the OpenCL, which was recently released as a standard for parallel programming across different parallel architectures. It also discusses other parallel technologies such as CUDA, ATI Stream and OpenMP. Subsequently, existing methods of electrostatic interaction calculation in molecular systems are reviewed. To demonstrate the performance and various features of OpenCL programming language, a basic algorithm for the calculation of electrostatic energy was implemented in OpenCL and tested on NVIDIA GeForce 9600 GT graphic card. The obtained results suggest that OpenCL is likely to be well applicable for electrostatic energy calculation. In the end, the prospects of OpenCL and observations made during the problem implementation are discussed. Abstrakt Úlohy řešené ve výpočetní chemii mohou mít značnou výpočetní složitost. Jedním ze způsobů, jak zrychlit výpočet těchto úloh, je řešit některé jejich části paralelně. Cílem této práce je zhodnotit možnosti jazyka OpenCL pro paralelizaci vědecko-technických výpočtů. Použitelnost jazyka je demonstrována na úloze výpočtu elektrostatických interakcí v molekulárních systémech. Teoretická část práce obsahuje základní úvod do problematiky paralelních výpočtů a popis jazyka OpenCL, nového standardu pro paralelní programování na různých počítačových architekturách. Je zahrnut i popis alternativních technologií jako je CUDA, ATI Stream a OpenMP. Dále jsou diskutovány existující metody pro výpočet elektrostatických interakcí v molekulárních systémech. Výkon a některé důležité vlastnosti jazyka OpenCL jsou demonstrovány na implementaci základního algoritmu pro výpočet elektrostatické potenciální energie. Testování výkonosti algoritmu je provedeno na grafické kartě NVIDIA GeForce 9600 GT. Výsledky testů naznačují velmi dobrou použitelnost jazyka OpenCL pro výpočet elektrostatické potenciální energie. Na závěr jsou zhodnoceny možnosti jazyka OpenCL a diskutovány postřehy získané během implementace úlohy. v Table of Contents Introduction............................................................................................................................1 Theory....................................................................................................................................2 Parallel Computing............................................................................................................2 OpenCL...........................................................................................................................10 Computation of Electrostatic Interactions in Molecular Systems...................................15 Results..................................................................................................................................18 Implementation................................................................................................................18 Performance.....................................................................................................................24 Conclusion............................................................................................................................30 Bibliography.........................................................................................................................32 Appendix A...........................................................................................................................35 OpenCL Kernels..............................................................................................................35 Appendix B...........................................................................................................................40 C/C++ Algorithm for Potential Energy...........................................................................40 Appendix C...........................................................................................................................42 Program Arguments.........................................................................................................42 Appendix D..........................................................................................................................43 Measurements..................................................................................................................43 vi Introduction Computational chemistry is a branch of chemistry that utilizes computers for modeling and simulation of molecular systems. Introduction of parallelism to some problems that are dealt with in this discipline can lead to significant speedup, especially when solutions that exploit massive parallelism available in graphic cards are used. An example of such solution is CUDA by NVIDIA or ATI Stream by AMD/ATI. Unfortunately, these solutions restrict the users to one platform, i.e. CUDA cannot be used on AMD graphic cards and vice versa. OpenCL is a standard for parallel programming across various parallel architectures including graphic cards. In other words, it enables to create program for the graphic cards developed by different companies. The standard was released on December 9, 2008. First implementations appeared in 2009. The standard is young and the implementations are slowly reaching their maturity. Hence, the thesis aims to assess the potential of employing OpenCL in acceleration of computational chemistry calculations, namely the computation of the potential energy due to electrostatic interactions in a molecular system. An introduction to parallel computing and the description of OpenCL standard is provided in chapter Theory. It also provides discussion about electrostatic interactions in molecular system and the way they are calculated. Chapter Results describes the implementation of the electrostatic energy calculation. It discusses the performance of the implementation as well. Finally, applicability of OpenCL is summarized and assessed in the last chapter Conclusion. 1 Theory Parallel Computing Introduction A computational problem can be usually solved in a serial or parallel way. Serial solution consists of discrete steps; each of them being done after the preceding one is finished. Most programs targeted at personal computers have serial nature: they basically consist of a set of instructions which are executed one after another on CPU. The other way is to solve a problem in a parallel manner. In that case, the problem is divided into parts and some of them are solved at the same time, i.e. in parallel. Of course, there are other ways of describing parallel computing. For example, as suggested in [1] "parallel computing is the simultaneous use of multiple compute resources to solve a computational problem." Application of parallel computing can lead to improving performance of many problems. For example, matrix multiplication or addition show significant performance increase. The reason is simple – subtasks of these operations demonstrate large degree of independence and therefore they can be done at the same time in parallel. However, problems with tasks that heavily depend on each other are not likely to perform well on parallel architectures. The scope for parallel processing varies greatly. It is possible to look at parallelism in the context of a processor, a computer or a set of computers connected through a network. Different techniques can be combined to work in concert. For example, a problem can be solved in parallel by multi-processor computers connected via network. The part of problem, which is to be solved in the context of single computer, can also take advantage of parallelism. To conclude, the topic of parallel computing is very broad. Due to the extent, the parallel computing will be further discussed in the context of single computer only. This attitude does not impose any restrictions on the thesis since the use of OpenCL standard fits well into this category. It should also be noted that there are many models and categories of parallelism, but only those that the author considers most relevant to the topic will be discussed. The concept of parallelism can be employed for several reasons. Clearly, a performance boost can be expected, if the problem consists of parts which demonstrate considerable degree of independence. In addition, according to [1]: "both physical and practical reasons pose significant constraints to simply building ever faster serial computers." Thus, paral- 2 lelism represents a way of avoiding limits of serial computing. Additionally, economical benefits can also be expected – a high performance processor can be more expensive than several relatively cheap standard processors which if combined achieve the same perform- ance. Important Issues Computation speedup may be considered as the main reason for introducing parallelism into a problem. Therefore, it is necessary to investigate speedup related topics, i.e. speedup limits or problem classification with respect to parallelism. Speedup limits In 1967 [2] Amdahl suggested a way for determining the speedup upper bound. The argument (also known as Amdahl's law) is based on the size of serial portion in a problem. The potential speedup can be defined [3] as speedup= sp s p N where s is the time needed for execution of serial part by a serial processor, p is is the time needed for execution of parallel part by a serial processor and N is the number of processors on which the problem is solved. With infinitely large N and simplification that s+p = 1, the reduced formula speedup= 1 s shows that the maximum speedup is limited by the serial portion of the program. As can be seen in Figure 1, one of the law's consequences is that increasing the number of processors does not necessarily yield the better performance. 3 Amdahl's law has been criticized by Gustafson in his paper [3] which serves as source of the following text. He states that the law is based on the assumption that the size of a problem is fixed – in other words: "p is independent of N." However, he adds that it "is virtually never the case" as the size of a practical problem "scales with the number of processors". Furthermore, Gustafson's findings suggest that the serial portion of a program remains relatively constant as "times for vector start-up, program loading, serial bottlenecks, and I/O that make up the … (serial) component of the run do not grow with problem size". In contrast with Amdahl, Gustafson proposes that both serial and parallel execution time should be measured in context of parallel system. Thus, the potential speedup can be defined as speedup=N s⋅1−N  where N is the number of processors and s is the serial time spent on the parallel system. To conclude, the new approach, which is called Gustafson-Barsis' law, does not limit the speedup as severely as does Amdahl's law. Problem Classification As mentioned above, to solve a problem in parallel it is necessary to divide it into parts which can be calculated simultaneously. However, these parts can be partially dependent on each other; therefore, they may need to synchronize their execution. The synchroniza- 4 Figure 1: Illustration of the impact of Amdahl's law on maximum speedup (source: Wikipedia) tion usually complicates the implementation and also increases overheads. Therefore, it is convenient to classify problems with respect to the degree of interconnection manifested by their parts. According to [4], a problem demonstrate • fine-grained parallelism if the subtasks communicate frequently; • coarse-grained parallelism if the communication among subtasks occurs in moderate rate; • embarrassing parallelism if their parts "rarely or never have to communicate." Types of Parallelism There are many classifications that address the topic of parallel computing. In following text, only topics related to this work will be discussed. Flynn's Taxonomy Flynn [5] suggested classification of computer architectures based on the concept of streams. A stream "means a sequence of items (instructions or data) as executed or operated on by a processor". The four categories result from the use of single or multiple streams of data or instructions: • SISD (single instruction stream, single data stream) Single instruction stream operates on single stream of data. Barney states [1] that "this is the oldest and even today, the most common type of computer". • SIMD (single instruction stream, multiple data streams) The same instruction stream operates on multiple streams of data. For example, modern graphics processors fall into this category. • MISD (multiple instruction streams, single data stream) Single stream of data is processed by multiple instruction streams. Few machines that would fit in this category existed [1], and this model is considered impractical [6]. • MIMD (multiple instruction streams, multiple data streams) According to Duncan [6, p. 6], this category "involves multiple processors autonomously executing diverse instructions on diverse data". Examples are modern supercomputers, networked parallel computer clusters or multi-core PCs [1]. Data Parallelism Data parallelism is characteristic [7] with "similar (not necessarily identical) operation sequences or functions being performed on elements of a large data structure". SIMD 5 concept is closely related to data parallelism. The data parallelism usually scales well [7] with the size of data parallel problem. Task Parallelism As the name suggests, task parallelism is focused on execution of different tasks which can have the same or different functionality. Another definition is suggested in [8]: "The term task parallelism … refers to the explicit creation of multiple threads of control, or tasks, which synchronize and communicate under programmer control." Unlike with data parallelism, the exploitable degree of task parallelism does not grow much with the size of a problem [7]. CPU versus GPU Traditionally, both CPUs and GPUs have been developed for different purposes. A CPU serves as a computational unit and provides a way to operate other computer components; its aim is to execute a task as fast as possible [9]. A GPU is specialized hardware for rendering computer graphics. Therefore, the GPU needs [9] to "be capable of processing a maximum of tasks, or to be more accurate, one task for a maximum of data in a minimum period of time". Bearing this in mind, it is no surprise that the architectures of both devices have developed differently in order to suit their needs. As a result, the GPUs show massive performance in certain types of computations, usually include significant number of SIMD processors and have native support for parallelism built in hardware. The CPUs, on the other hand, perform well in all kinds of computations, but tend to have rather moderate support of parallelism built in hardware. Figure 2 demonstrates parallelism available in GPUs. 6 Figure 2: Comparison of CPU and GPU architecture shows large availability of parallelism in GPU (source: NVIDIA) Selection of pre-OpenCL technologies OpenMP OpenMP (Open Multi-Processing) is an API specification [10] that supports multi-platform parallel programming in C, C++ and Fortran. The specification defines behavior of collection of compiler directives, library routines, and environment variables, thus forming the API. The fork-join execution model is employed by OpenMP (see Figure 3). The base idea of the model is that a program consists of master thread and regions in which the master thread forks into multiple threads. These multiple threads join back to the master thread once the execution of parallel region is finished. The multiple threads have their own private memory and can also access the memory shared among all of them. The parallel regions of a program are marked in the source code by compiler directives. ATI Stream ATI Stream is a technology developed by AMD/ATI which comprises "a set of advanced hardware and software technologies that enable AMD graphics processors (GPU), working in concert with the system’s central processor (CPU), to accelerate many applications beyond just graphics" [11]. The technology provides both software and hardware background which allows development of applications that take advantage of data parallelism available in AMD GPUs. The hardware side of technology is represented [12] by ATI Stream processor. The processor includes number of programmable stream cores which can execute user-defined functions, called kernels, on streams of data. The stream cores are physically grouped into thread processors which, in turn, are grouped into units called SIMD engines. Thread pro- 7 Figure 3: The fork-join execution model (source: OpenMP Tutorial by Blaise Barney) cessors within the same SIMD engine execute the same instruction sequence. Different engines, on the other hand, may execute different instruction sequences. The software includes compilers, stream processor drivers, libraries and support programming tools [12]. Brook+ is [13] a a high-level C based language that allows use of graphics chips for computation. A compiler for this language is part of the included software. The ATI Compute Abstraction Layer (CAL) is a device driver library for the GPU. According to [13], thanks to CAL "much of the hardware peculiarities of GPUs can be abstracted away". There are various advantages and disadvantages of the solution. On one hand, the performance boost can be expected in many applications. On the other hand, the use of the technology is limited only to the graphic processors manufactured by AMD. It should also be noted that the programming style differs considerably from the traditional one. For example, the memory access pattern has to be taken into account, since if done incorrectly it may cause significant performance degradation [12]. CUDA CUDA (Compute Unified Device Architecture) is a solution developed by NVIDIA that allows [14] utilization of NVIDIA GPUs' parallel capabilities for general purpose computing. Data parallel problems are likely to be addressed well by the technology. According to [15], the technology "allows developers to use C as a high-level programming language". The architecture offers support or plans to support in the future different languages and APIs such as C++, OpenCL, DirectX Compute, FORTRAN, etc. CUDA parallel programming model is based on user-defined functions called kernels, user-defined functions, that are executed in parallel by different threads. Each NVIDIA GPU contains certain number of multiprocessors which can execute multiple threads in parallel. The thread execution is managed automatically by hardware thread manager. Hence, the programmer does not have to treat parallel threads execution explicitly. The threads are grouped into one-dimensional, two-dimensional or three-dimensional blocks which are in turn part of one-dimensional or two-dimensional grid. When a thread block is to be executed, it is partitioned in warps which consists of 32 parallel threads that are created, managed, scheduled, and executed by the multiprocessor. Execution of threads in the same block can be synchronized via intrinsic function that "acts as a barrier at which all threads in the block must wait before any is allowed to proceed" [15, p. 8]. Each thread in the block has its own memory as well as access to block-shared memory and global memory. The global memory access is slowest while the thread memory access is fastest. There are two ways [15, p. 15] of writing programs for CUDA which are, however, mutually exclusive. The first is use of C for CUDA which is a C-based high-level program- 8 ming language with extensions that allow to define kernels as C functions. CUDA driver API, which represents the second option, is lower-level C API. Driver API provides high level of control, though it is more difficult to program with. The programming has two levels: managing NVIDIA graphic devices and execution parameters which is done on the host (CPU) and defining kernels that are executed on the NVIDIA graphic device or devices. When programming for CUDA, it is necessary [14] to consider the structure of data. The size of thread blocks can have significant impact on the performance. Hence, understanding of the structure may lead to optimal block size parameters. The advantages of the technology are a) the massive parallel power it opens for programmers, b) good support by NVIDIA and c) large number of existing solutions in CUDA which allow a beginner to grasp faster the subject. The older generations of graphic cards lack the double-precision float support which, however, is not always necessary. Even if the double-precision arithmetic is available the use of single-precision arithmetic is likely to lead to better performance [16, p. 35]. Another important issue concerning arithmetic is that it may not be as precise as it is on CPUs [15, pp. 115–117]. When using CUDA, the programmer is restricted only to NVIDIA graphic cards which can be viewed as a drawback. Due to the dependency of performance upon the structure of data and block size parameters, rather unfamiliar approach has to be considered. For example, data restructuring or different memory access pattern can greatly improve performance [15]. 9 OpenCL General Description OpenCL (Open Computing Language) is an open standard which aims [17, p. 11] to enable parallel programming across heterogeneous parallel processing platforms. The OpenCL specification [17] defines the standard, and a device can be said to support the standard, if an implementation that complies with the specification exists. Both data and task parallel programming models are supported, though the design was primarily driven [17, p. 26] with respect to the data parallel model. Programming for OpenCL consists of two parts. In the first part an API provided by standard is used to set execution parameters and specify the devices that should be used for computation. The second part lies in programming of kernels which are functions executed on a device or devices. The standard provides a language for the purpose of kernel programming. OpenCL C programming language is based [17, p. 126] on the C99 dialect of the C programming language. A set of restrictions and extensions is employed in order to reflect features specific to parallel programming. Due to combination of API and C-based programming language, both low-level abstraction and portability among various devices can be achieved. Another interesting feature of the standard is that it can interoperate with OpenGL and other graphics APIs. The concept of OpenCL was created by Apple, and in June 2008, with support of AMD, Intel and NVIDIA, it was proposed [18], [19] to the Khronos Group consortium as the basis for a new standard. Working group charged with the creation of a new standard was formed [20] in the same month. The group included [20] representatives of computer industry – processor developers as AMD, NVIDIA, Intel and IBM, software developers and others. The specification was released on December 9, 2008, ca. six months after the formation of the working group which illustrates interest of the industry in the standard. Architecture In order to understand how OpenCL works, it is necessary to investigate the way the standard views platform, execution of OpenCL program and memory hierarchy. The OpenCL specification [17] serves as the main information source for the following text. 10 Platform Model Parallel processors can vary greatly in their architecture. Therefore, in order to allow use of different processors, the standard introduces certain degree of abstraction into the view of the devices it works with. As Figure 4 shows, the platform is modeled as a host which is connected to one or more devices. A program that uses API provided by standard to manage the devices and kernel execution is executed on the host. A device is modeled as a set of one or more compute units which consists of one or more processing elements. For example, a GPU or multicore CPU can correspond to the device. Execution Model Execution of an OpenCL program can be divided into two parts. Host program which manages the devices and kernel execution is executed on the host while kernels are executed on one or more devices. Host program may query platform and devices for information that are necessary for creating context which is an environment for execution of kernels. The context comprises a set of devices that are to be used for kernel execution, memory available to the devices and mechanisms of memory management and kernel execution scheduling. Once the context is properly set, it is possible to create a program that will be executed in the context. The program consists of a set of kernels, optional auxiliary functions called by the kernels and con- 11 Figure 4: OpenCL Platform Model stants. The next step is program linking and compilation. OpenCL also allows usage of already compiled programs. After kernel parameters are set, kernels are executed. Kernel execution takes place in one-dimensional, two-dimensional or three-dimensional space which is defined when a kernel is submitted for execution. A point in space is called a work-item. Each work-item executes the same code and operates upon different or same data. A set of work-items forms a work-group. The work-items in the same work-group are executed concurrently on the processing elements of a single compute unit. They also have access to the work-group shared memory. It is also possible to synchronize them via special function that has to be executed by each work-item before execution can continue. Every work-item can be identified by its unique global ID or by a combination of its local ID, which identifies a work-item within its work-group, and work-group unique ID. Memory Model OpenCL memory model recognizes four types of memory that is available to kernels. The memory differs in access time, accessibility for kernels and the host and in their relationship to the OpenCL device model which is shown in Figure 5. Due to the abstract view of the device, it is only possible to make an assumption1 about the access time which may not always hold though and is not part of the specification. However, it has to be noted as the user's work with memory is greatly affected by it. The four types of memory defined in the model are as follows: • Global Memory The memory can be accessed by a host and every work-item. Access time is likely to be rather high. Access may be cached if supported by a device. • Constant Memory Constant memory is virtually the same as global memory, however, the data stored in it cannot by modified from kernels. Access time may be slightly lower than the one of global memory. • Local Memory This memory is unique for every work-group and cannot by accessed by a host. All work-items of the same work-group share access to the memory. Global memory access time should be considerably higher than local memory access time. • Private Memory The memory is unique for every work-item and cannot be accessed by anyone else. Access time can be expected to be lower than in all other cases. 1 The assumption is based on author's observations. 12 Optimization Due to the abstraction inherent to OpenCL platform model, it is rather difficult to present a strategy suitable for architecture of every device. However, as certain groups of devices share common features, consideration of these features may show a good strategy for these groups. For example, when using GPU as OpenCL device, the PCI bus may represent [21] the bottleneck for traffic between the device and the host. Therefore, minimization of data transfers may improve the performance. Since the OpenCL implementation is likely to make use of existing architecture, an optimization common to several architectures may be employed. For example, in case of GPUs, appropriate memory access pattern can increase performance. 13 Figure 5: Relationship of the OpenCL device architecture and memory model. Advantages and Disadvantages The main foreseeable problem lies in differences between CPU and GPU architectures. While GPUs may demonstrate optimal performance when executing large number of threads, multi-core CPUs, due to the lack of hardware-based massive parallelism support, tend not to cope well with large number concurrently running threads. Therefore, considerable performance variations per architecture can be encountered when executing a kernel designed to be massively parallel. Logically, diversification of OpenCL programs to programs for CPU and GPU might be a result. Architecture-specific optimization may also contribute to the possible trend. The main point of previous text is that OpenCL should not be thought of as a miraculous standard that will allow easy exploitation of parallelism without regard to the device architecture. Despite aforementioned issues, OpenCL holds advantages too. Due to the similarities in GPU architectures, the standard provides a way to create high-performance parallel programs for GPUs that are independent of graphics card manufacturer. The same can be applied to multi-core CPUs. Despite the possibility of need for different approach to CPUs and GPUs, a programmer does not need to learn new programming language or master new API. 14 Computation of Electrostatic Interactions in Molecular Systems Overview Computational chemistry is a branch of chemistry that "deals with the modeling and the computer simulation of systems such as biomolecules, polymers, drugs, inorganic and organic molecules, and so on" [22, p. 1]. The discipline often employes approximation or omission of certain properties to reduce the complexity of simulation. There are several methods which differ in the precision, the theory they are based on and their applicability. For example due to their complexity, rather precise methods derived from quantum mechanics tend to be used in relatively small systems [22, p. 6]. Molecular mechanics, on the other hand, omits some properties of the system and models [22, p. 7], [23, p. 165] the system as a simple algebraic expression for its total energy. As a result, molecular mechanics is often used to model systems with large number of atoms. Moreover, in some cases the method provides answers as accurate as those of quantum mechanics methods [23, p. 165]. The formula for total energy of a molecular system treated by molecular mechanics has usually the following form: E=Eintramolecular forcesEintermolecular forces (1) Energy due to the intermolecular forces arises from the non-bonded interactions between atoms in the system. The non-bonded interactions are van der Waals interactions and electrostatic interactions. Using Coulomb's potential, the electrostatic potential energy for system with N particles is defined as E= 1 40 ∑ i=1 N ∑ j=i1 N qi qj ri , j (2) where εₒ is vacuum permittivity, q is electric charge and r is distance between the interacting atoms. Due to the term 1/r in the Equation (2), electrostatic interactions are relevant over long distances. The computation of the total energy of a molecular system is rather difficult due to the electrostatic interactions term. As the Equation (2) shows, its computational complexity is O(N²). Hence, the straightforward approach with usual computational resources can be troublesome if applied to large systems. Therefore, several techniques that approximate the Equation (2) were developed. Most common techniques [23, p. 334] are the Ewald summation method, the reaction field method and the fast multipole method. 15 The Ewald Summation Method The method was developed [23, p. 334] by Ewald to study the energetics of ionic crystals. The system is viewed as a composition of regular cells. The cell array is considered to have a spherical shape in the limit. The distribution of particles in one cell is the same as in others. This approach is also known as periodical boundary condition [23, p. 317]. For simplicity, the assumed shape of cell will be a cube of side L containing N charges. One cell is declared to be the central cell. As a result, position of each image cell can be expressed by vector n = (iL, jL, kL) where i, j, k = 0, ±1, ±2, etc. Thus, the electrostatic potential energy due to the central cell can be expressed [23, p. 335], [24] as E= 1 2 ∑ ∣n∣≠0 ' ∑ i=1 N ∑ j=1 N qi q j 40∣ri , jn∣ (3) where the prime on first summation denotes that the interaction i = j is not included for n = 0, r is the position vector of a particle and ri , j=ri−r j . In Ewald method, the slowly convergent series of the Equation (3) is recast [23, pp. 334– 339], [24] into two series; each of which converges more rapidly. Additionally, one series deals with rapid variation of electrostatic interactions at small distances and the other takes care of slow decay at long distances. Hence, the Ewald sum can be expressed as E=Eshort−rangeElong−rangeEcorrection (4) Moreover, the summation of long-range term is performed in reciprocal space, hence the Fourier transform may be employed [23, p. 338]. Depending on approach to the summations, different computational complexity can be achieved. Direct summations with various parameters lead to O N 2  or O N N  complexity [23, p. 338]. However, the problem can be modified so that the fast Fourier transform is used to calculate the long-range term which with few adjustments to short-range term calculation leads to O N lnN  complexity [23, p. 338]. Due to the similarity between OpenCL and CUDA, it can be assumed that if an algorithm for a problem exists for CUDA it can be transformed to OpenCL. Since the implementation [25] for CUDA exists, it is highly probable that it can be achieved by OpenCL as well. Moreover, graphic cards can perform well even in calculations that manifest O N 2  complexity [26]. The Reaction Field Method The reaction field method [23, pp. 339–341] constructs a sphere around the molecule with given radius (also known as cutoff distance). The interactions with molecules within the sphere are calculated explicitly whereas the interactions with molecules outside the sphere 16 are modeled as a single interaction with a homogenous medium of given dielectric con- stant. Effective Calculation Of Short-range Electrostatic Interactions Both Ewald summation method and the reaction field method need to evaluate short-range electrostatic interactions. If the short-range surroundings of each atom was determined during the calculation of electrostatic interactions there would be really no benefit as the problem would still show O N 2  complexity. In other words, if the distance between the interacting atom and each other atom was to be computed then there would be no significant performance gain, and the electrostatic interactions between atoms could be calculated as well with minimal overhead [23, p. 325]. Hence, there has to be some memory overhead in order to quickly determine which atoms are within the cutoff. For example, in neighbor lists technique [23, p. 325] all atoms within the cutoff distance are stored in a list, together with the atoms that are slightly further away than the cutoff distance. Due to rather irregular memory access pattern, implementation of the neighbor-list technique is likely to be challenging on GPUs. However, there is an article [27] suggesting that such implementation is feasible. The Fast Multipole Method In the fast multipole method [23, pp. 341–343], [28], the simulation space, which is assumed to be a large cubic cell, is recursively subdivided into uniform cubic cells. The cell hierarchy is stored in an octree (a tree whose nodes have eight children or none). Electrostatic interactions of an atom with atoms inside its cell and atoms in nearby cells are calculated directly. The faraway interactions are treated as interactions between the cell containing the atom and the faraway cells. Furthermore, the cells that are very faraway are grouped into larger cells, which are then used instead the smaller cells in the cell-cell interactions. The grouping is done simply via the cell hierarchy. Based on the parameters of computation, the O N complexity can be achieved [28]. According to [23, p. 342], the algorithm "requires an amount of bookkeeping to keep track of the hierarchy of the cells, which means that up to certain size of problem the exact N² algorithm is faster". However, the implementation for GPU architecture exists [29], but the authors note that hierarchical structure "is not easily implemented on data-parallel pro- cessors". 17 Results Implementation Introduction The capabilities and features of OpenCL will be demonstrated on the calculation of the electrostatic potential energy of a molecular system. The chosen molecular system can be seen in Figure 6; it comprises atoms of restriction enzyme BsoBI surrounded by water molecules and counter ions. The restriction enzymes can act as molecular scissors, cutting DNA at specific spots [30]. Thus, they can be used to manipulate DNA. 18 Figure 6: Restriction enzyme BsoBI. The blue lines show the shape of water box surrounding the enzyme. One of the most important aspects of OpenCL is its capability to be executed on GPUs from different manufacturers. Therefore, the implementation of the energy calculation is focused on GPU platform. It was shown that the problems with O N 2  computational complexity can be solved efficiently on GPUs [26], therefore, the direct approach to the calculation of the electrostatic energy is taken. A brute-force algorithm in C/C++ for CPU is also implemented to provide basic comparison of CPU and GPU performance. It should be noted that there might be better performing algorithms for CPU, but this is not an issue as the purpose of the algorithm is to provide rough comparison with GPU platform. The OpenCL algorithm uses single-precision arithmetic, which has broader support on graphic cards and is considerably faster than double-precision arithmetic. The C/C++ algorithm uses double-precision arithmetic. For further details of the C/C++ algorithm see Appendix B. The implementation consists of two parts – a host program and an OpenCL kernel. The host program loads the given molecular system and sets the parameters of kernel execution such as work-group size, level of optimization features and profiling. It can also calculate the electrostatic potential on CPU. To enable performance comparison, both C/C++ and OpenCL algorithms support profiling. Only the time of computation is measured, i.e. the data preparations are not included. The profiling in OpenCL is managed via the functions that it provides. The C/C++ algorithm employs the standard C Time Library. Algorithm Description First, a parallelism that can be mapped on the execution model has to be identified in order to implement it on a graphic card. Fortunately, the calculation of electrostatic potential shows coarse-grained data parallelism. The Equation (2) can be also expressed as E=∑ j=2 N q1 q j r1, j ∑ j=3 N q2 qj r2, j ...∑ j= N N qN −1 qj rN −1, j  ∑ j=N 1 N qN q j r N , j (5) Thus, the potential energy can be viewed as a sum of energies arising from the interactions of the j-th atom with (j+1)-th atom, (j+2)-th atom, etc. Hence, each thread2 will compute this atom energy based on its global ID, i.e. the contribution of the i-th atom will be computed by thread with global ID value (i-1, due to indexing from zero). In a sense, the thread can also be viewed as the atom as it stores its coordinates and charge whose electrostatic contribution it calculates. Another step is to define the function of a work-group. The global memory access is a costly operation, hence it is necessary to minimize it. In a naive 2 The OpenCL term "work-item" would be more accurate, however, "thread" is likely to be more familiar to the reader. 19 approach, each thread, which represents an atom, would directly load information about all interacting atoms from the global memory. This would, however, result in a large number of global memory accesses, which could significantly degrade the performance. Hence, a smarter way to deal with this issue has to be devised. Due to the fact that local and global memory accesses can be synchronized within a work-group, a thread in a work-group can be used to load only some of the necessary data while loading of the rest is handled by the other threads in the work-group. With these principles in mind, the algorithm can be described (from the perspective of a thread in a work-group) as follows: • load the coordinates and charge of atom due to this thread and store them both in private memory and shared memory • compute the electrostatic energy due to the interactions of the atom owned by this thread with atoms in the same work-group whose global ID is larger than that of this atom • in batches load and store into the shared memory the coordinates and charges of atoms that have higher global ID and are represented by threads in other workgroups except for the last one; when a batch is loaded, compute the electrostatic energy due to the interactions of the atom this thread represents with all the loaded atoms • if not in the last work-group, compute the energy due to the interactions with atoms in the last work-group • based on the work-item ID store the computed energy contribution into the global memory The algorithm described above is illustrated in Figure 7. Once the contributions are calculated, the parallel summation based on [31] takes place. The summation kernel is executed several times until the summation result is obtained. See Appendix A for the source code of discussed kernels. 20 Optimization Features The optimization features incorporated into the kernel include coalesced memory access to the global memory, use of local memory, loop unrolling, problem size adjustment, separation of boundary situations. Coalesced Memory Access In order to perform global memory access in optimized fashion, both AMD and NVIDIA suggest that the access to memory addresses should be done sequentially across threads in a work-group [15], [32]. The implementation stores charges and coordinates of the atoms in separate arrays which are accessed in the aforementioned way. 21 Figure 7: The principle of the algorithm demnostrated on the execution of the first work-group. Each square represents both atom with its properties and a thread. The active threads/atoms whose contributions are being calculated are marked with blue color. The red rectangle shows the current group of atoms with which the atoms from the first work-group interact. Local Memory As mentioned before, the access to local memory is usually very fast. Therefore, the local memory is used to store the information about the atoms interacting with those that the threads in given work-group represent. As Listing 1 shows, all threads in the work-group access the same address in the shared memory when calculating interactions with the currently loaded atoms. This approach leads to optimal access on NVIDIA platform as the stored value is broadcast to all the accessing threads. Unfortunately, the AMD solution does not contain [32] such mechanism hence the access may not be optimal. Loop Unrolling Evaluation of a condition in a for loop statement can cause unnecessary overhead and thus be a costly operation. Due to this fact, both AMD [32] and NVIDIA [15] suggest to unroll the for loop, example of which is shown in Listing 1 (value of EI_LOOP_UNROLL macro is set in the host program). In the implementation, the loop unrolling factor is chosen by user. The permitted values are 2, 4 and 8 and the result of work-group size modulo loop unrolling factor has to be zero. If the user sets value to 1 there will not be any loop unrolling. Problem Size Adjustment The problem size is another issue that has to be considered. Generally, the GPU hardware is likely to partition the work-group in some way. Therefore, it is a good idea [16, p. 48], [32] to keep the work-group size a power of two or at least multiple of 32 or 64. Consequently, the problem size has to be adjusted to be a multiple of the work-group size. 22 // compute interactions with the loaded atoms for(unsigned int j=0; j=2 j++; locEI+=locCharge*sharedCharge[j]/sqrt( d2(locX,sharedX[j])+d2(locY,sharedY[j])+d2(locZ,sharedZ[j])); #endif … } Listing 1: Calculation of interactions with loaded atoms Therefore, the implementation rounds the problem size up to the nearest multiple of workgroup size. As a result, it is necessary to ensure that work done by these excess threads does not affect the calculation. Separation of Boundary Situations There are two situation that require special treatment. First is the calculation of interactions with the atoms represented by the same work-group. An atom cannot interact with itself and the atoms that have lower global ID. The second situation arises from the problem size adjustment as discussed in previous section. Of course, both situation could be handled by conditional statements in the main cycle. However as the situations does not occur often, much of computational power would be wasted on evaluation of the conditional statements. Therefore, both situations are excluded from the main cycle and treated separately. 23 Performance Basic Considerations The configuration of computer that was used for performance tests of the implementation is shown in Table 1. Unfortunately, no matching AMD/ATI graphic card was available for performance comparison. However, since the main optimization features due to the GPU platform are very similar, it can be assumed that the implementation would perform similarly on AMD/ATI graphic card. Of course, minor performance increase might be obtained from platform-specific optimizations. Table 1: Hardware and software configuration Hardware CPU Intel Core 2 Duo CPU E8400 @ 3.00GHz RAM 4096 MB GPU NVIDIA GeForce 9600 GT 512MB Software Operating System Mandriva Linux 2010.0 64bit NVIDIA Driver Version 195.36.15 First, the performance will be tested on the aforementioned molecular system (enzyme BsoBI). In order to further test the implementation, tests on a system with different number of atoms will follow. However, the system will only consist of water molecules whose number will be gradually changed. The CPU performance will refer to the performance of C/C++ algorithm, whereas GPU performance will refer to the OpenCL implementation. As noted before, the electrostatic energy computation is handled by two OpenCL kernels. Originally, the execution time of the both kernels was measured. However, early test results showed that the execution time of the summation kernel is negligible in comparison with the execution time of the main kernel. Hence, only the total execution time will be considered in the tests. The early results also showed that the optimal work-group size is a power of two, and that from the available work-group sizes, 128 and 256 give very similar results. In the end, the main kernel work-group size was chosen to be 128 and the summation kernel work-group size was chosen to be 256. 24 Test Results Enzyme System The system containing the enzyme was tested with different loop unrolling factors to show its impact on performance and to find optimal loop factor for further testing. The results are shown in Table 2. For further details about results see Appendix D. Table 2: Test results of the enzyme system with 117,555 atoms GPU Results Loop Unrolling Factor Average Time (s) Standard Deviation 1 1.712 1.30E-005 2 1.400 1.56E-004 4 1.221 2.35E-004 8 1.688 3.53E-005 Calculated Electrostatic Energy: -7,849,986 kcal/mol CPU Results Average Time (s) Standard Deviation 124.657 1.42E-001 Calculated Electrostatic Energy: -7,849,986.153 kcal/mol Results Comparison Loop Unrolling Factor GPU Time (s) CPU Time (s) CPU/GPU Time Ratio 1 1.712 124.657 72.8 2 1.400 124.657 89.0 4 1.221 124.657 102.1 8 1.688 124.657 73.9 As shown in Figure 8, the optimal loop unrolling factor is four for the problem solved on GeForce 9600 GT hardware. The performance increase is likely due to faster execution of the work-groups. Due to the fact that graphic cards favor large number of concurrent computations, four is likely the optimal loop unrolling factor for systems with large number of atoms. 25 Depending on the loop unrolling factor, the GPU calculation is from 70 to 100 times faster than the CPU calculation as can be seen in Figure 9. Despite the fact that the used CPU algorithm is not probably the fastest one, the gained boost suggests that even if faster algorithm was used the performance gain would still be significant. Another issue that has to be addressed is the algorithm precision. As shown in Table 2, the single-precision OpenCL calculation gives -7,849,986 kcal/mol as a result, whereas double-precision CPU algorithm calculates result -7,849,986.153 kcal/mol. If the doubleprecision result is truncated to single-precision it is the same as the single-precision OpenCL result. Thus, it is possible to say that the single-precision result for given problem size is as precise as the one calculated in double precision. 26 Figure 8: Relationship of loop unrolling factor and performance 0 1 2 3 4 5 6 7 8 9 0.9 1.1 1.3 1.5 1.7 1.9 2.1 Loop Unrolling Factor Time Figure 9: CPU/GPU Ratio 1 2 4 8 0 20 40 60 80 100 120 Loop Unrolling Factor CPU/GPURatio Water System First, a water system larger than the enzyme system was tested in order to examine whether optimal loop factor determined in previous test is really suitable even for larger systems. Subsequently, the performance of the implementation was tested on systems with varying number of atoms. The detailed results are provided in Appendix D. Table 3: Effects of different loop unrolling factors on performance for system containing 541,956 atoms GPU Results Loop Unrolling Factor Average Time (s) Standard Deviation 1 35.992 2.77E-004 2 29.410 2.25E-003 4 25.651 9.87E-004 8 35.589 4.21E-004 As mentioned before, the loop unrolling factor which leads to the best performance should hold for systems with different size. As can be seen in Table 3, the optimal factor for a system with 541,956 atoms is the same as in previous case – four. This fact further supports the conjecture. As a result, this loop unrolling factor will be used for all subsequent tests. As Table 4 shows, the performance increase over CPU is stable with various number of atoms in the system. The results also suggests that for given hardware configuration, the reasonable problem size is ca. 140,000 atoms. Since the used graphic card is of the older generation, it is likely that the execution on newer graphic cards will perform well on systems with far more than 140,000 atoms. Another issue that has to be mentioned is precision. As can be seen in Table 4, the difference between energy calculated on CPU and GPU cannot be accounted to the fact that GPU result is the same as CPU result but truncated to single precision. The average error 0.0024% due to GPU energy likely stems from the use of single-precision arithmetic and the reduced precision of division and square root operations [31, pp. 53–57]. 27 Table 4: Test results of water systems with various number of atoms (electrostatic energy is given in kcal/mol) Time Comparison Number Of Atoms GPU Time (s) CPU Time (s) CPU/GPU Time Ratio 43,239 0.170 16.850 99.3 140,460 1.739 178.827 102.9 239,922 5.047 524.007 103.8 339,216 10.066 1,052.733 104.6 447,627 17.505 1,836.933 104.9 541,956 25.651 2,671.853 104.2 630,828 34.737 3,623.580 104.3 Calculated Electrostatic Energy Comparison Number Of Atoms GPU Energy CPU Energy Error 43,239 -3,073,207 -3,073,207.50 0.0016% 140,460 -9,991,924 -9,991,926.11 0.0021% 239,922 -17,072,500 -17,072,500.43 0.0003% 339,216 -24,143,448 -24,143,450.51 0.0010% 447,627 -31,864,542 -31,864,534.63 0.0023% 541,956 -38,582,628 -38,582,615.02 0.0034% 630,828 -44,909,884 -44,909,857.80 0.0058% Average Error: 0.0024% In order to devise a way to determine how well the implementation scales with the problem size, there has to be some assumptions made about the relationship between computation time and problem size. Due to the O N 2  complexity of the algorithm, the relationship can be viewed as t≈kN 2 (6) where k denotes the coefficient that depends on the used architecture and problem implementation. If the coefficient k is constant for different problem sizes N then the implementation is likely to scale well with the problem size. Due to the involved power of two, it is, however, difficult to see whether it scales quadratically or not, if the devised relationship is captured in graph. In order to fix this problem, the relationship can be transformed to linear dependency on the problem size: t≈k N (7) 28 Figure 10 shows the discussed relationship. Since the observed trend is linear, the implementation on the given GPU likely scales well with the problem size. Results Discussion It was shown that the implementation provides good precision for systems with up to 300,000 or 400,000 atoms. With respect to the age of GeForce 9600 GT, the implementation is likely suitable for systems containing up to 600,000 atoms; possibly less if there is need for high precision. The implementation also demonstrates stable performance boost over CPU and good scalability. In a broader sense, the implementation has shown important features of OpenCL and general purpose graphic cards programming. For example, loop unrolling or memory access pattern, which does not usually have to be considered when programming for CPUs, are critical on GPUs. The results also show that very good scalability and performance increase can be achieved in comparison with CPU if the algorithm is implemented properly Conclusion Problems of computational chemistry can be computationally exhaustive. One way to deal with this fact is to introduce parallelism to the problem. There are different technologies that allow parallel processing in context of single computer, namely ATI Stream, CUDA, 29 0 200000 400000 600000 800000 0 1 2 3 4 5 6 7 Number Of Atoms SquareRootOfTime OpenMP and newly OpenCL, a standard for parallel computing on heterogeneous platforms. This thesis discusses the features of this technologies. However, it is mainly focused on the usage of OpenCL for calculation of electrostatic interactions in molecular systems. To illustrate the way the OpenCL is used, a brute-force algorithm for calculation of electrostatic interactions on GPU with O N 2  complexity was implemented. Important optimization techniques such as loop unrolling, usage of shared memory and coalesced global memory access were demonstrated. The tests were conducted on molecular systems with the number of atoms ranging from ca. 40,000 to 630,000. Even though the implemented OpenCL algorithm uses single-precision arithmetic the calculated results lead to average error 0.0024% in comparison with double-precision brute-force CPU algorithm. Moreover, the testing also proved OpenCL implementation to be 100 times faster than basic CPU algorithm. To conclude, the results showed good scalability with problem size and promising performance for problems with O N 2  complexity. Future Works Since only the basic O N 2  algorithm was implemented, next move is to implement some of the described algorithm with O N complexity. In addition, there are many comparisons that could be done. For example, if AMD/ATI and NVIDIA graphics cards were obtained it would be very interesting to determine whether an implementation that performs optimally on both platforms can be achieved. Another interesting topic might be OpenCL performance and usability on IBM Cell processor. OpenCL Assessment Based on Personal Experience It is a fact that the same OpenCL code can be executed on both CPUs and GPUs. However, due to the architecture-specific features, the algorithms are likely to be aimed at the respective architectures or at the groups of architectures with similar features. The use of OpenCL on traditional multi-core processors (i.e. dual-cores, quad-cores, etc.) is quite debatable. Since there is matured OpenMP standard that allows parallel processing. In addition, in contrast with OpenMP standard, OpenCL requires fair amount of preparative coding in order to set and start the computation and the programming model is quite different from the traditional one. On the other hand, there might be cases in which OpenCL is probably more suitable, IBM's Cell processor being an example. The one of the main advantages of OpenCL is probably its ability to run on the GPUs by different manufacturers. Hence, if the goal is to make a program that can tap into the graphic card's performance and run on large number of computers then OpenCL is probably the best option. However, the situation may change for scientific calculations. If per- 30 formance is the main goal and the code is required to run really fast on a particular computer NVIDIA's CUDA might be a better option. Firstly, since it is a proprietary solution it is likely to react to new trends and evolve faster which can lead to easier program development and better performance. Secondly, NVIDIA offers considerable support and as of 2010 programming for CUDA is widespread among research community. Another issue to address is the user friendliness and ease of use. The provided API is quite robust and offers broad functionality. On the other hand, this approach also leads to great deal of coding before a kernel can be executed. There are experimental C++ bindings available, however, their support is not official. Hopefully, future versions of the standard will address these issues. From other features for example, the built-in profiling support is very convenient. When programing for GPU, it is rather unfortunate that there are virtually none debugging options. This issue is addressed by NVIDIA Parallel Nsight which is however limited only to the NVIDIA graphic cards. It will also be probably primarily focused on CUDA rather then on OpenCL. In conclusion, the standard will surely find its application and its use can be recommended in scientific calculations. Although, the user should be always aware of other applicable technologies and the differences in architectures. 31 Bibliography [1] B. Barney, "Introduction to Parallel Computing," Jan. 26, 2009. [Online]. Available: https://computing.llnl.gov/tutorials/parallel_comp/#WhyUse. [Accessed: Dec. 18, 2009]. [2] G.M. Amdahl, "Validity of the single processor approach to achieving large scale computing capabilities," Proceedings of the April 18-20, 1967, spring joint computer conference, Atlantic City, New Jersey: ACM, 1967, pp. 483-485. [Online]. Available: ACM, http://portal.acm.org/. [Accessed Dec. 10, 2009]. [3] J.L. Gustafson, "Reevaluating Amdahl's law," Commun. ACM, vol. 31, 1988, pp. 532- 533. [Online]. Available: ACM, http://portal.acm.org/. [Accessed Dec. 20, 2009]. [4] A. Saudi, "Parallel Computing: Introduction," April 2008. Available: http://www.azalisaudi.com/para/Para-Week1-Intro.pdf. [Accessed: Dec. 21, 2009] [5] M.J. Flynn, "Some Computer Organizations and Their Effectiveness," Computers, IEEE Transactions on, vol. C-21, 1972, pp. 948-960. [Online]. Available: IEEE Xplore, http://ieeexplore.ieee.org/. [Accessed Dec. 21, 2009]. [6] R. Duncan, "A survey of parallel computer architectures," Computer, vol. 23, 1990, pp. 5-16. [Online]. Available: IEEE Xplore, http://ieeexplore.ieee.org/. [Accessed Dec. 21, 2009]. [7] D.E. Culler, J.P. Singh, and A. Gupta, Parallel computer architecture, Gulf Professional Publishing, 1999, pp. 124-125. [E-book] Available: Google Books [8] G. Perrin and A. Darte, The data parallel programming model, Springer, 1996, pp. 10- 13. [E-book] Available: Google Books [9] D. Triolet, "Nvidia CUDA: preview," BeHardware, Mar. 21, 2007. [Online]. Available: http://www.behardware.com/articles/659-1/nvidia-cuda-preview.html. [Accessed: Dec. 22, 2009]. [10] OpenMP Architecture Review Board, "OpenMP Application Program Interface," openmp.org, May 2008. [Online]. Available: http://www.openmp.org/mpdocuments/spec30.pdf. [Accessed: Dec. 22, 2009]. [11] AMD, "ATI Stream SDK v2.0," Advanced Micro Devices. [Online]. Available: http://developer.amd.com/GPU/ATISTREAMSDK/Pages/default.aspx. [Accessed: Dec. 22, 2009]. [12] AMD, "Technical Overview: ATI Stream Computing," Advanced Micro Devices, Mar. 20, 2009. [Online]. Available: http://developer.amd.com/gpu/ATIStreamSDK/pages/Publications.aspx. [Accessed: Dec. 23, 2009]. 32 [13] A. Bailey, "Get Extraordinary Performance By Exploiting the GPU," Advanced Micro Devices, Aug . 13, 2008. [Online]. Available: http://developer.amd.com/documentation/articles/pages/GetExtraordinaryPerformanc eByExploitingtheGPU.aspx. [Accessed: Dec. 23, 2009]. [14] T. Halfhill, "Parallel Processing With CUDA," Microprocessor Report, Jan. 28, 2008. [Online]. Available: http://www.nvidia.com/docs/IO/55972/220401_Reprint.pdf. [Accessed: Jan. 22, 2010]. [15] NVIDIA, "NVIDIA CUDA Programming Guide," NVIDIA, Aug. 26, 2009. [Online]. Available: http://developer.download.nvidia.com/compute/cuda/2_3/toolkit/docs/NVIDIA_CUD A_Programming_Guide_2.3.pdf. [Accessed: Jan. 22, 2010]. [16] NVIDIA "Optimization: NVIDIA OpenCL Best Practices Guide Version 2.3," NVIDIA, Aug. 31, 2009. [Online]. Available: http://developer.download.nvidia.com/compute/cuda/3_0/toolkit/docs/NVIDIA_Open CL_BestPracticesGuide.pdf. [Accessed: Jan. 22, 2010]. [17] Aaftab Munshi, ed., "The OpenCL Specification," Khronos Group, Oct. 6, 2009. [Online]. Available: http://www.khronos.org/registry/cl/specs/opencl-1.0.48.pdf. [Accessed: Nov. 10, 2009]. [18] Apple, "OpenCL: Taking the graphics processor beyond graphics," Apple, August 2009. [Online]. Available: http://images.apple.com/macosx/technology/docs/OpenCL_TB_brief_20090903.pdf. [Accessed: Jan. 24, 2010]. [19] Achim Basermann, "OpenCL: The Open Standard for Heterogeneous Parallel Programming," July 2009. Available: http://www.slideshare.net/onyame/opencl-theopen-standard-for-heterogeneous-parallel-programming. [Accessed: Jan. 24, 2010]. [20] Khronos Group, "Khronos Press Releases: Khronos Launches Heterogeneous Computing Initiative," June 16, 2008. Available: http://www.khronos.org/news/press/releases/khronos_launches_heterogeneous_comp uting_initiative/. [Accessed: Nov. 10, 2009]. [21] NVIDIA, "NVIDIA OpenCL - Best Practices Guide," NVIDIA, Aug. 16, 2009. [Online]. Available: http://www.nvidia.com/content/cudazone/CUDABrowser/downloads/papers/NVIDIA _OpenCL_BestPracticesGuide.pdf. [Accessed: Jan. 26, 2010]. [22] K.I. Ramachandran, G. Deepa, and K. Namboori, Computational Chemistry and Molecular Modeling: Principles and Applications, Springer, 2008. [23] A.R. Leach, Molecular modelling, Pearson Education, 2001. [24] A.Y. Toukmaji and J.A. Board, "Ewald summation techniques in perspective: a 33 survey," Computer Physics Communications, vol. 95, Jun. 1996, pp. 73-92. [Online]. Available: ScienceDirect, http://sciencedirect.com/. [Accessed Apr. 27, 2010]. [25] M.J. Harvey and G. De Fabritiis, "An Implementation of the Smooth Particle Mesh Ewald Method on GPU Hardware," Journal of Chemical Theory and Computation, vol. 5, 2009, pp. 2371-2377. [Online]. Available: ACS Publications. [Accessed Apr. 27, 2010]. [26] L. Nyland, M. Harris, and J. Prins, "Fast N-Body Simulation with CUDA," GPU Gems 3, Addison Wesley Professional, 2007. [27] W. Liu, B. Schmidt, G. Voss, and W. Müller-Wittig, "Accelerating molecular dynamics simulations using Graphics Processing Units with CUDA," Computer Physics Communications, vol. 179, Nov. 2008, pp. 634-641. [Online]. Available: ScienceDirect, http://sciencedirect.com/. [Accessed May 1, 2010]. [28] D. Blackston and T. Suel, "Highly portable and efficient implementations of parallel adaptive N-body methods," Proceedings of the 1997 ACM/IEEE conference on Supercomputing (CDROM), San Jose, CA: ACM, 1997, pp. 1-20. [Online]. Available: ACM, http://portal.acm.org/. [Accessed May 2, 2010]. [29] N.A. Gumerov and R. Duraiswami, "Fast multipole methods on graphics processors," J. Comput. Phys., vol. 227, 2008, pp. 8290-8313. [Online]. Available: ACM, http://portal.acm.org/. [Accessed May 2, 2010]. [30] David S. Goodsell, "Molecule of the Month: Restriction Enzymes," pdb.org, August 2000. [Online]. Available: http://www.rcsb.org/pdb/static.do? p=education_discussion/molecule_of_the_month/pdb8_1.html. [Accessed May 11, 2010]. [31] NVIDIA, "OpenCL Programming for the CUDAArchitecture Version 2.3," NVIDIA, Aug. 31, 2009. [Online]. Available: http://www.nvidia.com/content/cudazone/download/OpenCL/NVIDIA_OpenCL_Pro grammingOverview.pdf. [Accessed May 4, 2010]. [32] AMD "Technical Notes: ATI Stream SDK v2.01 Performance and Optimization," Advanced Micro Devices, 2010. [Online]. Available: http://developer.amd.com/gpu/ATIStreamSDK/assets/ATI_Stream_SDK_Performanc e_Notes.pdf. [Accessed May 4, 2010]. 34 Appendix A OpenCL Kernels #define d2(a,b) (((a)-(b))*((a)-(b))) // (a-b)^2 /*Each thread represents an atom with unique ID and computes and sums its all electrostatic interactions with all atoms that have higher ID than the atom. Summation of results of all threads is electrostatic enegry of given system. */ __kernel void computeEI( __global float * X, __global float * Y, __global float * Z, __global float * charge, unsigned int atomCount, __local float * sharedX, __local float * sharedY, __local float * sharedZ, __local float * sharedCharge, __global float * EI){ size_t gID = get_global_id(0); size_t lID = get_local_id(0); unsigned int groupCount = get_num_groups(0); unsigned int groupSize = get_local_size(0); // first, the interactions with atoms, which other // threads in the work-group represent, are calculated sharedX[lID] = X[gID]; sharedY[lID] = Y[gID]; sharedZ[lID] = Z[gID]; sharedCharge[lID] = charge[gID]; barrier(CLK_LOCAL_MEM_FENCE); // load the characteristics of the atom which the thread represents float locX = sharedX[lID]; float locY = sharedY[lID]; float locZ = sharedZ[lID]; float locCharge = sharedCharge[lID]; float locEI = 0; if (gID < atomCount) { for(unsigned int j=0; jlID){ locEI+= locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); } // loop unrolling may lead to faster execution #if defined(EI_LOOP_UNROLL) && EI_LOOP_UNROLL>=2 j++; if ((get_group_id(0)*groupSize+j)lID){ locEI+= locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); } #endif #if defined(EI_LOOP_UNROLL) && EI_LOOP_UNROLL>=4 j++; if ((get_group_id(0)*groupSize+j)lID){ locEI+= locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); } j++; if ((get_group_id(0)*groupSize+j)lID){ locEI+= locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); } #endif #if defined(EI_LOOP_UNROLL) && EI_LOOP_UNROLL>=8 j++; if ((get_group_id(0)*groupSize+j)lID){ locEI+= locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); } j++; if ((get_group_id(0)*groupSize+j)lID){ locEI+= locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); } j++; if ((get_group_id(0)*groupSize+j)lID){ locEI+= locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); } j++; if ((get_group_id(0)*groupSize+j)lID){ locEI+= locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); } #endif } } barrier(CLK_LOCAL_MEM_FENCE); 36 // calculates interactions with other atoms that // have higher IDs than atoms represented by threads // in this work-group // // note: in order to omit conditional statements (out-of-bound check) // the interactions with atoms due to the last work-group // are not calculated for (unsigned int i=get_group_id(0)+1; i=2 j++; locEI+=locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); #endif #if defined(EI_LOOP_UNROLL) && EI_LOOP_UNROLL>=4 j++; locEI+=locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); j++; locEI+=locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); #endif #if defined(EI_LOOP_UNROLL) && EI_LOOP_UNROLL>=8 j++; locEI+=locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); j++; locEI+=locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); j++; locEI+=locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); j++; locEI+=locCharge*sharedCharge[j]/ sqrt(d2(locX,sharedX[j]) + d2(locY,sharedY[j]) + d2(locZ,sharedZ[j])); 37 #endif } } barrier(CLK_LOCAL_MEM_FENCE); } // calculates the interactions with the atoms due to the last work-group if (get_group_id(0)!=groupCount-1){ sharedX[lID] = X[lID + (groupCount-1)*groupSize]; sharedY[lID] = Y[lID + (groupCount-1)*groupSize]; sharedZ[lID] = Z[lID + (groupCount-1)*groupSize]; sharedCharge[lID] = charge[lID + (groupCount-1)*groupSize]; barrier(CLK_LOCAL_MEM_FENCE); if (gID < atomCount) { for(unsigned int j=0; j=2 j++; if ((j+(groupCount-1)*groupSize)=4 j++; if ((j+(groupCount-1)*groupSize)=8 j++; if ((j+(groupCount-1)*groupSize)> 1 = n/2 // optimal access to shared memory on NVIDIA GPUs for(int step = get_local_size(0)>>1; step > 0; step>>=1) { if (lID < step) { shared[lID] += shared[lID+step]; } barrier(CLK_LOCAL_MEM_FENCE); } if (lID == 0) { out[get_group_id(0)] = shared[0]; } } 39 Appendix B C/C++ Algorithm for Potential Energy /*Function computes the electrostatic energy on CPU using brute force, the execution time may be obtained via execTime argument.*/ double computeOnCPU(double * execTime) { #define deltaSqr(a,b) (((a)-(b))*((a)-(b))) // load double precision data for calculation // note: system files are loaded and checked in main() int numOfAtoms = topology.AtomList.GetNumberOfAtoms(); double * dblX = (double *) malloc(sizeof(double)*numOfAtoms); double * dblY = (double *) malloc(sizeof(double)*numOfAtoms); double * dblZ = (double *) malloc(sizeof(double)*numOfAtoms); double * dblCharge = (double *) malloc(sizeof(double)*numOfAtoms); for(unsigned int i = 0; iGetCharge(); } // calculates electrostatic energy of given system clock_t startCPU, endCPU; startCPU = clock(); double all = 0; double ei = 0; double x; double y; double z; double chrg; for (unsigned int i = 0; i