Search

Purpose The current work aims to present a parallel code using the open multi-processing (OpenMP) programming model for an adaptive multi-resolution high-order finite difference scheme for solving 2D conservation laws, comparing efficiencies obtained with a previous message passing interface formulation for the same serial scheme and considering the same type of 2D formulations laws.Design/methodology/approach The serial version of the code is naturally suitable for parallelization because the spatial operator formulation is based on a splitting scheme per direction for which the flux components are numerically computed by a Lax–Friedrichs factorization independently for each row or column. High-order approximations for numerical fluxes are computed by the third-order essentially non-oscillatory (ENO) and fifth-order weighted essentially non-oscillatory (WENO) interpolation schemes, assuming sparse grids in each direction. The grid adaptivity is obtained by a cubic interpolating wavelet transform applied in each space dimension, associated to a threshold operator. Time is evolved by a third order TVD Runge–Kutta method.Findings The parallel formulation is implemented automatically at compiling time by the OpenMP library routines, being virtually transparent to the programmer. This over simplifies any concerns about managing and/or updating the adaptive grid when compared to what is necessary to be done when other parallel approaches are considered. Numerical simulations results and the large speedups obtained for the Euler equations in gas dynamics highlight the efficiency of the OpenMP approach.Research limitations/implications The resulting speedups reflect the effectiveness of the OpenMP approach but are, to a large extension, limited by the hardware used (2 E5-2620 Intel Xeon processors, 6 cores, 2 threads/core, hyper-threading enabled). As the demand for OpenMP threads increases, the code starts to make explicit use of the second logical thread available in each E5-2620 processor core and efficiency drops. The speedup peak is reached near the possible maximum (24) at about 22, 23 threads. This peak reflects the hardware configuration and the true software limit should be located way beyond this value.Practical implications So far no attempts have been made to parallelize other possible code segments (for instance, the ENO|-WENO-TVD code lines that process the different data components which could potentially push the speed up limit to higher values even further. The fact that the speedup peak is located close to the present hardware limit reflects the scalability properties of the OpenMP programming and of the splitting scheme as well. Consequently, it is likely that the speedup peak with the OpenMP approach for this kind of problem formulation will be close to the physical (and/or logical) limit of the hardware used.Social implications This work is the result of a successful collaboration among researchers from two different institutions, one internationally well-known and with a long-term experience in applied mathematics for industrial applications and the other in a starting process of international academic insertion. In this way, this scientific partnership has the potential of promoting further knowledge exchange, involving students and other collaborators.Originality/value The proposed methodology (use of OpenMP programming model for the wavelet adaptive splitting scheme) is original and contributes to a very active research area in the past years, namely, adaptive methods for conservation laws and their parallel formulations, which is of great interest for the entire scientific community. [ABSTRACT FROM AUTHOR]

In modern physics it has become common to find the solution of a problem by solving numerically a set of PDEs. Whether solving them on a finite difference grid or by a finite element approach, the main calculations are often applied to a stencil structure. In the last decade it has become usual to work with so called big data problems where calculations are very heavy and accelerators and modern architectures are widely used. Although CPU and GPU clusters are often used to solve such problems, parallelisation of any calculation ideally starts from a single processor optimisation. Unfortunately, it is impossible to vectorise a stencil structured loop with high level instructions. In this paper we suggest a new approach to rearranging the data structure which makes it possible to apply high level vectorisation instructions to a stencil loop and which results in significant acceleration. The suggested method allows further acceleration if shared memory APIs are used. We show the effectiveness of the method by applying it to an elastic wave propagation problem on a finite difference grid. We have chosen Intel architecture for the test problem and OpenMP (Open Multi-Processing) since they are extensively used in many applications. [ABSTRACT FROM AUTHOR]

To optimize short-range force computations in Molecular Dynamics (MD) simulations, multi-threading and SIMD optimizations are presented in this paper. With respect to multi-threading optimization, a Partition-and-Separate-Calculation (PSC) method is designed to avoid write conflicts caused by using Newton’s third law. Serial bottlenecks are eliminated with no additional memory usage. The method is implemented by using the OpenMP model. Furthermore, the PSC method is employed on Intel Xeon Phi coprocessors in both native and offload models. We also evaluate the performance of the PSC method under different thread affinities on the MIC architecture. In the SIMD execution, we explain the performance influence in the PSC method, considering the “ if-clause ” of the cutoff radius check. The experiment results show that our PSC method is relatively more efficient compared to some traditional methods. In double precision, our 256-bit SIMD implementation is about 3 times faster than the scalar version. [ABSTRACT FROM AUTHOR]

In current computer architectures, data movement (from die to network) is by far the most energy consuming part of an algorithm (≈20 pJ/word on-die to ≈10,000 pJ/word on the network). To increase memory locality at the hardware level and reduce energy consumption related to data movement, future exascale computers tend to use many-core processors on each compute nodes that will have a reduced clock speed to allow for efficient cooling. To compensate for frequency decrease, machine vendors are making use of long SIMD instruction registers that are able to process multiple data with one arithmetic operator in one clock cycle. SIMD register length is expected to double every four years. As a consequence, Particle-In-Cell (PIC) codes will have to achieve good vectorization to fully take advantage of these upcoming architectures. In this paper, we present a new algorithm that allows for efficient and portable SIMD vectorization of current/charge deposition routines that are, along with the field gathering routines, among the most time consuming parts of the PIC algorithm. Our new algorithm uses a particular data structure that takes into account memory alignment constraints and avoids gather/scatter instructions that can significantly affect vectorization performances on current CPUs. The new algorithm was successfully implemented in the 3D skeleton PIC code PICSAR and tested on Haswell Xeon processors (AVX2-256 bits wide data registers). Results show a factor of×2 to×2.5 speed-up in double precision for particle shape factor of orders 1-3. The new algorithm can be applied as is on future KNL (Knights Landing) architectures that will include AVX-512 instruction sets with 512 bits register lengths (8 doubles/16 singles). [ABSTRACT FROM AUTHOR]

We present new versions of the previously published C and CUDA programs for solving the dipolar Gross–Pitaevskii equation in one, two, and three spatial dimensions, which calculate stationary and non-stationary solutions by propagation in imaginary or real time. Presented programs are improved and parallelized versions of previous programs, divided into three packages according to the type of parallelization. First package contains improved and threaded version of sequential C programs using OpenMP. Second package additionally parallelizes three-dimensional variants of the OpenMP programs using MPI, allowing them to be run on distributed-memory systems. Finally, previous three-dimensional CUDA-parallelized programs are further parallelized using MPI, similarly as the OpenMP programs. We also present speedup test results obtained using new versions of programs in comparison with the previous sequential C and parallel CUDA programs. The improvements to the sequential version yield a speedup of 1.1–1.9, depending on the program. OpenMP parallelization yields further speedup of 2–12 on a 16-core workstation, while OpenMP/MPI version demonstrates a speedup of 11.5–16.5 on a computer cluster with 32 nodes used. CUDA/MPI version shows a speedup of 9–10 on a computer cluster with 32 nodes. New version program summary Program Title: DBEC-GP-OMP-CUDA-MPI: (1) DBEC-GP-OMP package: (i) imag1dX-th, (ii) imag1dZ-th, (iii) imag2dXY-th, (iv) imag2dXZ-th, (v) imag3d-th, (vi) real1dX-th, (vii) real1dZ-th, (viii) real2dXY-th, (ix) real2dXZ-th, (x) real3d-th; (2) DBEC-GP-MPI package: (i) imag3d-mpi, (ii) real3d-mpi; (3) DBEC-GP-MPI-CUDA package: (i) imag3d-mpicuda, (ii) real3d-mpicuda. Program Files doi: http://dx.doi.org/10.17632/j3z9z379m8.1 Licensing provisions: Apache License 2.0 Programming language: OpenMP C; CUDA C. Computer: DBEC-GP-OMP runs on any multi-core personal computer or workstation with an OpenMP-capable C compiler and FFTW3 library installed. MPI versions are intended for a computer cluster with a recent MPI implementation installed. Additionally, DBEC-GP-MPI-CUDA requires CUDA-aware MPI implementation installed, as well as that a computer or a cluster has Nvidia GPU with Compute Capability 2.0 or higher, with CUDA toolkit (minimum version 7.5) installed. Number of processors used: All available CPU cores on the executing computer for OpenMP version, all available CPU cores across all cluster nodes used for OpenMP/MPI version, and all available Nvidia GPUs across all cluster nodes used for CUDA/MPI version. Journal reference of previous version: Comput. Phys. Commun. 195 (2015) 117; ibid. 200 (2016) 406. Does the new version supersede the previous version?: Not completely. OpenMP version does supersede previous AEWL_v1_0 version, while MPI versions do not supersede previous versions and are meant for execution on computer clusters and multi-GPU workstations. Nature of problem: These programs are designed to solve the time-dependent nonlinear partial differential Gross–Pitaevskii (GP) equation with contact and dipolar interaction in a harmonic anisotropic trap. The GP equation describes the properties of a dilute trapped Bose–Einstein condensate. OpenMP package contains programs for solving the GP equation in one, two, and three spatial dimensions, while MPI packages contain only three-dimensional programs, which are computationally intensive or memory demanding enough to require such level of parallelization. Solution method: The time-dependent GP equation is solved by the split-step Crank–Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation, in either imaginary or real time, over small time steps. The contribution of the dipolar interaction is evaluated by a Fourier transformation to momentum space using a convolution theorem. MPI parallelization is done using the domain decomposition. The method yields the solution of stationary and/or non-stationary problems. Reasons for the new version: Previously published C and Fortran programs [1] for solving the dipolar GP equation are sequential in nature and do not exploit the multiple cores or CPUs found in typical modern computers. A parallel implementation exists, using Nvidia CUDA [2], and both versions are already used within the ultra-cold atoms community [3]. However, CUDA version requires special hardware, which limits its usability. Furthermore, many researchers have access to high performance computer clusters, which could be used to either further speed up the computation, or to work with problems which cannot fit into a memory of a single computer. In light of these observations, we have parallelized all programs using OpenMP, and then extended the parallelization of three-dimensional programs using MPI to distributed-memory clusters. Since the CUDA implementation uses the same algorithm, and thus has the same structure and flow, we have applied the same data distribution scheme to provide the distributed-memory CUDA/MPI implementation of three-dimensional programs. Summary of revisions: Package DBEC-GP-OMP: Previous serial C programs [1] are here improved and then parallelized using OpenMP (package DBEC-GP-OMP). The main improvement consists of switching to real-to-complex (R2C) Fourier transform, which is possible due to the fact that input of the transform is purely real. In this case the result of the transform has Hermitian symmetry, where one half of the values are complex conjugates of the other half. The fast Fourier transformation (FFT) libraries we use can exploit this to compute the result faster, using half the memory. To parallelize the programs, we have used OpenMP with the same approach as described in [4], and extended the parallelization routines to include the computation of the dipolar term. The FFT, used in computation of the dipolar term, was also parallelized in a straightforward manner, by using the built-in support for OpenMP in FFTW3 library [5]. With the introduction of multiple threads memory usage has increased, driven by the need to have some variables private to each thread. To reduce the memory consumed, we resorted to using techniques similar to the ones used in our CUDA implementation [2], i.e., we have reduced the memory required for FFT by exploiting the aforementioned R2C FFT, and reused the memory with pointer aliases whenever possible. Package DBEC-GP-MPI: Next step in the parallelization (package DBEC-GP-MPI) was to extend the programs to run on distributed-memory systems, i.e., on computer clusters using domain decomposition with MPI programming paradigm. We chose to use the newly-implemented threaded versions of the programs as the starting point. Alternatively, we could have used serial versions, and attempt a pure MPI parallelization, however we have found that OpenMP-parallelized routines better exploit the data locality and thus outperform the pure MPI implementation. Therefore, our OpenMP/MPI-parallelized programs are intended to run one MPI process per cluster node, and each process would spawn the OpenMP threads as needed on its cluster node. Note that this is not a requirement, and users may run more than one MPI process per node, but we advise against it due to performance reasons. With the suggested execution strategy (one MPI process per cluster node, each spawning as many threads as CPU cores available), OpenMP threads perform most of the computation, and MPI is used for data exchanges between processes. There are numerous ways to distribute the data between MPI processes, and we decided to use a simple one-dimensional data distribution, also known as slab decomposition. Data is distributed along the first (slowest changing) dimension, which corresponds to NX spatial dimension in our programs (see Fig. 1 ). Each process is assigned a different portion of the NX dimension, and contains the entire NY and NZ spatial dimensions locally. This allows each process to perform computation on those two dimensions in the same way as before, without any data exchanges. In case the computation requires whole NX dimension to be local to each process, we transpose the data, and after the computation, we transpose the data back. Fig. 1 Illustration of data distribution between MPI processes. On the left, the data are distributed along the NX dimension, while on the right the same data are redistributed along the NY dimension. Transpose routine can be implemented in many ways using MPI, most commonly using MPI_Alltoall function, or using transpose routines from external libraries, like FFTW3 [5] or 2DECOMP&FFT [6]. Since we already rely on FFTW3 library for FFT, we have utilized its dedicated transpose interface to perform the necessary transformations. To speed up transpose operation, we do not perform full transposition of data, but rather leave it locally transposed. That is, we transform from local_NX × NY × NZ, stored in row-major order, to NX × local_NY × NZ in row-major order (where local_NX = NX / number_of_processes, and equivalently for local_NY). This approach has an additional benefit that we do not have to make significant changes in the way array elements are processed, and in most cases we only have to adjust the loop limit of the non-local dimension. Package DBEC-GP-MPI-CUDA: The aforementioned data distribution scheme can be also applied to the CUDA version of programs [2]. However, there is no support for CUDA in FFTW3, and cuFFT (used in CUDA programs for FFT) does not provide equivalent MPI or transpose interface. Instead, we developed our own transpose routines, and used them in FFT computation. One example of manual implementation of transpose routines is shown in Ref. [7], and while we could readily use the same code, we wanted to have the same result as when using FFTW3. To achieve this, we use the same basic principle as in Ref. [7], first we create a custom MPI data type that maps to portions of the data to be exchanged, followed by an all-to-all communication to exchange the data between processes, see Fig. 2 for details. Fig. 2 Example of a transpose routine of a 4×4×4 data between four MPI processes. Initially, all processes have 1 / 4 of the NX dimension, and whole NY and NZ dimensions. After transposing, each process has full NX and NZ dimensions, and 1 / 4 of the NY dimension. The implemented transpose routines are also used to compute a distributed-memory FFT, performed over all MPI processes. To divide the computation of a multidimensional FFT, in our case three-dimensional, we use a well-known row–column algorithm. The basic idea of the algorithm is perhaps best explained on a two-dimensional FFT of N × M data, stored in row-major order, illustrated in Fig. 3 . First the N one-dimensional FFTs of length M are performed (along the row of data), followed by a transpose, after which data are stored as M × N in row-major format. Now M FFTs of length N can be performed along what used to be a column of original data, but are stored as rows after transposing. Finally, an optional transpose can be performed to return the data in their original N × M form. In three dimensions, we can perform a two-dimensional FFT, transpose the data, and perform the FFT along the third dimension. This algorithm can be easily adapted for distributed memory systems. We use advanced cuFFT interface for local computation of FFT, and use our transpose routine to redistribute the data. Note that DBEC-GP-MPI-CUDA programs can be easily modified to work on a single workstation with multiple GPU cards, or a computer cluster with multiple GPU cards per node. In that case, for each GPU card a separate MPI process should be launched and the programs should be modified to assign a separate GPU card for processes on the same cluster node. Fig. 3 Illustration of four stages of row–column FFT algorithm. The last transpose operation may be omitted, and often yields better performance. MPI output format: Given that the distributed memory versions of the programs can be used for much larger grid sizes, the output they produce (i.e., the density profiles) can be much larger and difficult to handle. To alleviate this problem somewhat, we have switched to a binary output instead of the textual. This allowed us to reduce the size of files, while still retaining precision. All MPI processes will write the output to the same file, at the corresponding offset, relieving the user of the task of combining the files. The binary output can be subsequently converted to textual, for example by using hexdump command on UNIX-like systems. We have developed a simple script which converts the output from binary to textual format and included it in the software package. Testing results: We have tested all programs on the PARADOX supercomputing facility at the Scientific Computing Laboratory of the Institute of Physics Belgrade. Nodes used for testing had two Intel Xeon E5-2670 CPUs (with a total of 2 × 8 = 16 CPU cores) with 32 GB of RAM and one Nvidia Tesla M2090 GPU with 6 GB of RAM, each connected by Infiniband QDR interconnect. The presented results are obtained for arbitrary grid sizes, which are not tailored to maximize performance of the programs. We also stress that execution times and speedups reported here are calculated for critical parallelized parts of the programs performing iterations over imaginary or real time steps, and they exclude time spent on initialization (threads initialization, MPI environment, allocation/deallocation of memory, creating/destroying FFTW plans, I/O operations). As a part of its output, each program separately prints initialization time and time spent on iterations for GP propagation. The latter time is used to calculate a speedup, as a speedup obtained this way does not depend on the number of iterations and is more useful for large numbers of iterations. The testing of OpenMP versions of programs DBEC-GP-OMP was performed with the number of threads varying from 1 to 16. Table 1 and Fig. 4 show the obtained absolute wall-clock times, speedups, and scaling efficiencies, as well as comparison with the previous serial version of programs [1]. As we can see from the table, improvements in the FFT routine used already yield a speedup of 1.3 to 1.9 for single-threaded ( T = 1 ) 2d and 3d programs compared to the previous serial programs, and somewhat smaller speedup for 1d programs, 1.1 to 1.3. The use of additional threads brings about further speedup of 2 to 2.5 for 1d programs, and 9 to 12 for 2d and 3d programs. From Fig. 4 we see that for 1d programs, although speedup increases with the number of threads used, the efficiency decreases due to insufficient size of the problem, and one can achieve almost maximal value of speedup already with T = 4 threads, while still keeping the efficiency around 50%. We also see, as expected, that speedup and efficiency of 2d and 3d programs behave quite well as we increase the numbers of threads. In particular, we note that the efficiency is always above 60%, making the use of all available CPU cores worthwhile. Fig. 4 Speedup in the execution time and scaling efficiency of DBEC-GP-OMP programs compared to single-threaded runs: (a) imag1dX-th, (b) real1dX-th, (c) imag2dXY-th, (d) real2dXY-th, (e) imag3d-th, (f) real3d-th. Scaling efficiency is calculated as a fraction of the obtained speedup compared to a theoretical maximum. Grid sizes used for testing are the same as in Table 1 . Speedups and efficiencies of imag1dZ-th, real1dZ-th, imag2dXZ-th, and real2dXZ-th (not reported here) are similar to those of imag1dX-th, real1dX-th, imag2dXY-th, and real2dXY-th, respectively. For testing of MPI versions we have used a similar methodology to measure the strong scaling performance. For OpenMP/MPI programs DBEC-GP-MPI, the obtained wall-clock times are shown in Table 2 , together with the corresponding wall-clock times for the OpenMP programs DBEC-GP-OMP that served as a baseline to calculate speedups. The testing was done for varying number of cluster nodes, from 4 to 32, and the measured speedup ranged from 11 to 16.5. The corresponding graphs of speedups and efficiencies are shown in Fig. 5 , where we can see that the speedup grows linearly with the number of nodes used, while the efficiency remains mostly constant in the range between 40% and 60%, thus making the use of OpenMP/MPI programs highly advantageous for problems with large grid sizes. Fig. 5 Speedup in the execution time and scaling efficiency of DBEC-GP-MPI programs compared to single-node OpenMP runs: (a) imag3d-mpi, (b) real3d-mpi. Scaling efficiency is calculated as a fraction of the obtained speedup compared to a theoretical maximum. Grid size used for testing is the same as in Table 2 . For CUDA/MPI programs DBEC-GP-MPI-CUDA we observe similar behavior in Table 3 and in Fig. 6 . The obtained speedup with N = 32 nodes here ranges from 9 to 10, with the efficiency between 30% and 40%. While the efficiency is slightly lower than in the case of OpenMP/MPI programs, which could be expected due to a more complex memory hierarchy when dealing with the multi-GPU system distributed over many cluster nodes, the speedup still grows linearly and makes CUDA/MPI programs ideal choice for use on GPU-enabled computer clusters. Additional benefit of using these programs is their low CPU usage (up to one CPU core), allowing for the possibility that same cluster nodes are used for other CPU-intensive simulations. Fig. 6 Speedup in the execution time and scaling efficiency of DBEC-GP-MPI-CUDA programs compared to single-node runs of previous CUDA programs [2]: (a) imag3d-mpicuda, (b) real3d-mpicuda. Scaling efficiency is calculated as a fraction of the obtained speedup compared to a theoretical maximum. Grid size used for testing is the same as in Table 3 . The introduction of distributed transposes of data creates some overhead, which negatively impacts scaling efficiency. This is more evident in the CUDA/MPI version, as the transpose algorithm is inferior to the one provided by FFTW3. In our tests, both MPI versions of programs failed to achieve speedup on less than 4 nodes, due to the introduction of the transpose routines. We therefore recommend using MPI versions only on 4 or more cluster nodes. The MPI versions are highly dependent not only on the configuration of the cluster, mainly on the speed of interconnect, but also on the distribution of processes and threads, NUMA configuration, etc. We recommend that users experiment with several different configurations to achieve the best performance. The results presented are obtained without extensive tuning, with the aim to show the base performance. Finally, we note that the best performance can be achieved by evenly distributing the workload among the MPI processes and OpenMP threads, and by using grid sizes which are optimal for FFT. In particular, the programs in DBEC-GP-OMP package have the best performance if NX, NY, and NZ are divisible by the number of OpenMP threads used. Similarly, for DBEC-GP-MPI programs the best performance is achieved if NX and NY are divisible by a product of the number of MPI processes and the number of OpenMP threads used. For DBEC-GP-MPI-CUDA programs, the best performance is achieved if NX and NY are divisible by a product of the number of MPI processes and the number of Streaming Multiprocessors (SM) in the GPU used. For all three packages, the best FFT performance is obtained if NX, NY and NZ can be expressed as 2 a 3 b 5 c 7 d 1 1 e 1 3 f , where e and f are either 0 or 1, and the other exponents are non-negative integer numbers [8]. Additional comments, restrictions, and unusual features: MPI programs require that grid size (controlled by input parameters NX, NY and NZ) can be evenly distributed between the processes, i.e., that NX and NY are divisible by the number of MPI processes. Since the data is never distributed along the NZ dimension, there is no such requirement on NZ. Programs will test if these conditions are met, and inform the user if not (by reporting an error). Additionally, MPI versions of CUDA programs require CUDA-aware MPI implementation. This allows the MPI runtime to directly access GPU memory pointers and avoid having to copy the data to main RAM. List of CUDA-aware MPI implementations can be found in Ref. [9]. Acknowledgments V.L., S.Š. and A.B. acknowledge support by the Ministry of Education, Science, and Technological Development of the Republic of Serbia under projects ON171017 , OI1611005 , and III43007 , as well as SCOPES project IZ74Z0-160453 . L.E. Y.-S. acknowledges support by the FAPESP of Brazil under project 2012/21871-7 and 2014/16363-8 . P.M. acknowledges support by the Science and Engineering Research Board, Department of Science and Technology, Government of India under project No. EMR/2014/000644 . S.K.A. acknowledges support by the CNPq of Brazil under project 303280/2014-0 , and by the FAPESP of Brazil under project 2012/00451-0 . References: [1] R. Kishor Kumar, L. E. Young-S., D. Vudragović, A. Balaž, P. Muruganandam, and S. K. Adhikari, Comput. Phys. Commun. 195 (2015) 117. [2] V. Lončar, A. Balaž, A. Bogojević, S. Škrbić, P. Muruganandam, S. K. Adhikari, Comput. Phys. Commun. 200 (2016) 406. [3] R. Kishor Kumar, P. Muruganandam, and B. A. Malomed, J. Phys. B: At. Mol. Opt. Phys. 46 (2013) 175302; H. Al-Jibbouri, I. Vidanović, A. Balaž, and A. Pelster, J. Phys. B: At. Mol. Opt. Phys. 46 (2013) 065303; R. R. Sakhel, A. R. Sakhel, and H. B. Ghassib, J. Low Temp. Phys. 173 (2013) 177; B. Nikolić, A. Balaž, and A. Pelster, Phys. Rev. A 88 (2013) 013624; X. Antoine and R. Duboscq, Comput. Phys. Commun. 185 (2014) 2969; J. Luo, Commun. Nonlinear Sci. Numer. Simul. 19 (2014) 3591; K.-T. Xi, J. Li, and D.-N. Shi, Physica B 436 (2014) 149; S. K. Adhikari, Phys. Rev. A 90 (2014) 055601; M. C. Raportaru, J. Jovanovski, B. Jakimovski, D. Jakimovski, and A. Mishev, Rom. J. Phys. 59 (2014) 677; A. I. Nicolin, A. Balaž, J. B. Sudharsan, and R. Radha, Rom. J. Phys. 59 (2014) 204; A. Balaž, R. Paun, A. I. Nicolin, S. Balasubramanian, and R. Ramaswamy, Phys. Rev. A 89 (2014) 023609; A. I. Nicolin and I. Rata, High-Performance Computing Infrastructure for South East Europe’s Research Communities: Results of the HP-SEE User Forum 2012, in Springer Series: Modeling and Optimization in Science and Technologies 2 (2014) 15; S. K. Adhikari, J. Phys. B: At. Mol. Opt. Phys. 48 (2015) 165303; S. K. Adhikari, Phys. Rev. E 92 (2015) 042926; T. Khellil and A. Pelster, arXiv:1512.04870 (2015); H. L. C. Couto and W. B. Cardoso, J. Phys. B: At. Mol. Opt. Phys. 48 (2015) 025301; R. R. Sakhel, A. R. Sakhel, and H. B. Ghassib, Physica B 478 (2015) 68; L. Salasnich and S. K. Adhikari, Acta Phys. Pol. A 128 (2015) 979; X. Antoine and R. Duboscq, Lecture Notes Math. 2146 (2015) 49; E. Chiquillo, J. Phys. A: Math. Theor. 48 (2015) 475001; S. Sabari, C. P. Jisha, K. Porsezian, and V. A. Brazhnyi, Phys. Rev. E 92 (2015) 032905; W. Wen, T. K. Shui, Y. F. Shan, and C. P. Zhu, J. Phys. B: At. Mol. Opt. Phys. 48 (2015) 175301; P. Das and P. K. Panigrahi, Laser Phys. 25 (2015) 125501; Y. S. Wang, S. T. Ji, Y. E. Luo, and Z. Y. Li, J. Korean. Phys. Soc. 67 (2015) L1504; A. I. Nicolin, M. C. Raportaru, and A. Balaž, Rom. Rep. Phys. 67 (2015) 143; V. S. Bagnato, D. J. Frantzeskakis, P. G. Kevrekidis, B. A. Malomed, and D. Mihalache, Rom. Rep. Phys. 67 (2015) 5; J. B. Sudharsan, R. Radha, H. Fabrelli, A. Gammal, and B. A. Malomed, Phys. Rev. A 92 (2015) 053601; K.-T. Xi, J. Li, and D.-N. Shi, Physica B 459 (2015) 6; E. J. M. Madarassy and V. T. Toth, Phys. Rev. D 91 (2015) 044041; F. I. Moxley III, T. Byrnes, B. Ma, Y. Yan, and W. Dai, J. Comput. Phys. 282 (2015) 303; D. Novoa, D. Tommasini, J. A. Nóvoa-López, Phys. Rev. E 91 (2015) 012904; Y. H. Wang, A. Kumar, F. Jendrzejewski, R. M. Wilson, M. Edwards, S. Eckel, G. K. Campbell, and C. W. Clark, New J. Phys. 17 (2015) 125012; T. Khellil, A. Balaž, and A. Pelster, New J. Phys. 18 (2016) 063003; T. Khellil and A. Pelster, J. Stat. Mech.-Theory Exp. (2016) 063301; J. Akram and A. Pelster, Phys. Rev. A 93 (2016) 023606; S. K. Adhikari, Laser Phys. Lett. 13 (2016) 035502; J. Akram and A. Pelster, Phys. Rev. A 93 (2016) 033610; J. Akram, B. Girodias, and A. Pelster, J. Phys. B: At. Mol. Opt. Phys. 49 (2016) 075302; S. K. Adhikari and S. Gautam, Phys. Rev. A 93 (2016) 013630; Ž. Marojević, E. Göklü, and C. Lämmerzahl, Comput. Phys. Commun. 202 (2016) 216; A. Paredes and H. Michninel, Phys. Dark Universe 12 (2016) 50; J. Akram and A. Pelster, Laser Phys. 26 (2016) 065501; T. Mithun, K. Porsezian, and B. Dey, Phys. Rev. A 93 (2016) 013620; C.-Y. Lai and C.-C. Chien, Phys. Rev. Appl. 5 (2016) 034001; S. K. Adhikari, Laser Phys. Lett. 13 (2016) 085501; K. Manikandan, P. Muruganandam, M. Senthilvelan, and M. Lakshmanan, Phys. Rev. E 93 (2016) 032212; R. R. Sakhel, A. R. Sakhel, H. B. Ghassib, and A. Balaž, Eur. Phys. J. D 70 (2016) 66; W. Bao, Q. Tang, and Y. Zhang, Commun. Comput. Phys. 19 (2016) 1141; R. Kishor Kumar, T. Sriraman, H. Fabrelli, P. Muruganandam, and A. Gammal, J. Phys. B: At. Mol. Opt. Phys. 49 (2016) 155301; A. Bogojević, A. Balaž, and A. Belić, Phys. Rev. E 72 (2005) 036128; A. Bogojević, I. Vidanović, A. Balaž, and A. Belić, Phys. Lett. A 372 (2008) 3341; I. Vidanović, A. Bogojević, A. Balaž, and A. Belić, Phys. Rev. E 80 (2009) 066706; A. Balaž, A. Bogojević, I. Vidanović, and A. Pelster, Phys. Rev. E 79 (2009) 036701; A. Balaž, I. Vidanović, A. Bogojević, and A. Pelster, Phys. Lett. A 374 (2010) 1539; A. I. Nicolin, Physica A 391 (2012) 1062; I. Vasić and A. Balaž, arXiv:1602.03538 (2016); O. Voronych, A. Buraczewski, M. Matuszewski, and M. Stobińska, arXiv:1603.02570 (2016); A. M. Martin, N. G. Marchant, D. H. J. O’Dell, and N. G. Parker, arXiv:1606.07107 (2016). [4] D. Vudragović, I. Vidanović, A. Balaž, P. Muruganandam, and S. K. Adhikari, Comput. Phys. Commun. 183 (2012) 2021. [5] FFTW3 library, http://www.fftw.org/ (2016). [6] 2DECOMP&FFT library, http://www.2decomp.org/ (2016). [7] B. Satarić, V. Slavnić, A. Belić, A. Balaž, P. Muruganandam, S. K. Adhikari, Comput. Phys. Commun. 200 (2016) 411. [8] Real-data DFTs with FFTW3, http://www.fftw.org/fftw3_doc/Real_002ddata-DFTs.html (2014); Nvidia’s cuFFT accuracy and performance, http://docs.nvidia.com/cuda/cufft/#accuracy-and-performance (2015). [9] Nvidia’s MPI Solutions for GPUs, https://developer.nvidia.com/mpi-solutions-gpus (2016). [ABSTRACT FROM AUTHOR]

The modern scientific research requires careful simulation of experiments, as well as rapid and qualitative processing of a large amount of the program code. Optimization of calculations implies both algorithm improvement and usage of high-speed parallel computing technologies. This article compares various ways of parallelizing programs that are used in the ROOT software package. [ABSTRACT FROM AUTHOR]

A modern Fortran implementation of three Dirac operators (Wilson, Brillouin, Susskind) in lattice QCD is presented, based on OpenMP shared-memory parallelization and SIMD pragmas. The main idea is to apply a Dirac operator to N v vectors simultaneously, to ease the memory bandwidth bottleneck. All index computations are left to the compiler and maximum weight is given to portability and flexibility. The lattice volume, N x N y N z N t , the number of colors, N c , and the number of right-hand sides, N v , are parameters defined at compile time. Several memory layout options are compared. The code performs well on modern many-core architectures (480 Gflop/s, 880 Gflop/s, and 780 Gflop/s with N v = 12 for the three operators in single precision on a 72-core KNL processor, a 2 × 24 -core Skylake node yields similar results). Explicit run-time tests with CG/BiCGstab inverters confirm that the memory layout is relevant for the KNL, but less so for the Skylake architecture. The ancillary code distribution contains all routines, including the single, double, and mixed precision Krylov space solvers, to render it self-contained and ready-to-use. Program title: testknl_main.f90 CPC Library link to program files: https://doi.org/10.17632/7km9wcspn3.1 Licensing provisions: Creative Commons by 4.0 Programming language: fortran 2008 Supplementary material: Detailed description of each module of the code distribution. Nature of problem: Solving the Dirac equation D u = b for three choices of the Dirac operator D relevant to Lattice QCD. The choices are known as Wilson fermions, Brillouin fermions and Susskind or staggered fermions, respectively. Solution method: Standard solution through sparse iterative solvers CG and BiCGstab, albeit with several right-hand vectors packed into a right-hand matrix. The number of colors, N c , the number of right-hand sides, N v , and the box lengths, N x , N y , N z and N t , can be chosen freely. Additional comments including restrictions and unusual features: The code allows the user to compare several memory layout options for the block of right-hand sides, that is several orderings of the slots pertinent to N c , N v and possibly the Dirac spinor, in the array representation of u and b. No machine dependent optimizations are used, only standard OpenMP shared-memory parallelization and SIMD pragmas are used. [ABSTRACT FROM AUTHOR]

We present new version of previously published Fortran and C programs for solving the Gross–Pitaevskii equation for a Bose–Einstein condensate with contact interaction in one, two and three spatial dimensions in imaginary and real time, yielding both stationary and non-stationary solutions. To reduce the execution time on multicore processors, new versions of parallelized programs are developed using Open Multi-Processing (OpenMP) interface. The input in the previous versions of programs was the mathematical quantity nonlinearity for dimensionless form of Gross–Pitaevskii equation, whereas in the present programs the inputs are quantities of experimental interest, such as, number of atoms, scattering length, oscillator length for the trap, etc. New output files for some integrated one- and two-dimensional densities of experimental interest are given. We also present speedup test results for the new programs. New version program summary Program title: BEC-GP-OMP package, consisting of: (i) imag1d, (ii) imag2d, (iii) imag3d, (iv) imagaxi, (v) imagcir, (vi) imagsph, (vii) real1d, (viii) real2d, (ix) real3d, (x) realaxi, (xi) realcir, (xii) realsph. Catalogue identifier: AEDU_v4_0. Program Summary URL: http://cpc.cs.qub.ac.uk/summaries/AEDU_v4_0.html Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland. Licensing provisions: Apache License 2.0 No. of lines in distributed program, including test data, etc.: 130308. No. of bytes in distributed program, including test data, etc.: 929062. Distribution format: tar.gz. Programming language: OpenMP C; OpenMP Fortran. Computer: Any multi-core personal computer or workstation. Operating system: Linux and Windows. RAM: 1 GB. Number of processors used: All available CPU cores on the executing computer. Classification: 2.9, 4.3, 4.12. Catalogue identifier of previous version: AEDU_v1_0, AEDU_v2_0. Journal reference of previous version: Comput. Phys. Commun. 180 (2009) 1888; ibid. 183 (2012) 2021. Does the new version supersede the previous version?: No. It does supersedes versions AEDU_v1_0 and AEDU_v2_0, but not AEDU_v3_0, which is MPI-parallelized version. Nature of problem: The present OpenMP Fortran and C programs solve the time-dependent nonlinear partial differential Gross–Pitaevskii (GP) equation for a Bose–Einstein condensate in one (1D), two (2D), and three (3D) spatial dimensions in a harmonic trap with six different symmetries: axial- and radial-symmetry in 3D, circular-symmetry in 2D, and fully anisotropic in 2D and 3D. Solution method: The time-dependent GP equation is solved by the split-step Crank–Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation, in either imaginary or real time, over small time steps. The method yields the solution of stationary and/or non-stationary problems. Reasons for the new version: Previously published Fortran and C programs [1,2] for solving the GP equation are recently enjoying frequent usage [3] and application to a more complex scenario of dipolar atoms [4] . They are also further extended to make use of general purpose graphics processing units (GPGPU) with Nvidia CUDA [5] , as well as computer clusters using Message Passing Interface (MPI) [6] . However, a vast majority of users use single-computer programs, with which the solution of a realistic dynamical 1D problem, not to mention the more complicated 2D and 3D problems, could be time consuming. Now practically all computers have multicore processors, ranging from 2 up to 18 and more CPU cores. Some computers include motherboards with more than one physical CPU, further increasing the possible number of available CPU cores on a single computer to several tens. The present programs are parallelized using OpenMP over all the CPU cores and can significantly reduce the execution time. Furthermore, in the old version of the programs [1,2] the inputs were based on the mathematical quantity nonlinearity for the dimensionless form of the GP equation. The inputs for the present versions of programs are given in terms of phenomenological variables of experimental interest, as in Refs. [4,5] , i.e., number of atoms, scattering length, harmonic oscillator length of the confining trap, etc. Also, the output files are given names which make identification of their contents easier, as in Refs. [4,5] . In addition, new output files for integrated densities of experimental interest are provided, and all programs were thoroughly revised to eliminate redundancies. Summary of revisions: Previous Fortran [1] and C [2] programs for the solution of time-dependent GP equation in 1D, 2D, and 3D with different trap symmetries have been modified to achieve two goals. First, they are parallelized using OpenMP interface to reduce the execution time in multicore processors. Previous C programs [2] had OpenMP-parallelized versions of 2D and 3D programs, together with the serial versions, while here all programs are OpenMP-parallelized. Secondly, the programs now have input and output files with quantities of phenomenological interest. There are six trap symmetries and both in C and in Fortran there are twelve programs, six for imaginary-time propagation and six for real-time propagation, totaling to 24 programs. In 3D, we consider full radial symmetry, axial symmetry and full anisotropy. In 2D, we consider circular symmetry and full anisotropy. The structure of all programs is similar. For the Fortran programs the input data (number of atoms, scattering length, harmonic oscillator trap length, trap anisotropy, etc.) are conveniently placed at the beginning of each program. For the C programs the input data are placed in separate input files, examples of which can be found in a directory named input. The examples of output files for both Fortran and C programs are placed in the corresponding directories called output. The programs then calculate the dimensionless nonlinearities actually used in the calculation. The provided programs use physical input parameters that give identical nonlinearity values as the previously published programs [1,2] , so that the output files of the old and new programs can be directly compared. The output files are conveniently named so that their contents can be easily identified, following Refs. [4,5] . For example, file named < code > -out.txt, where < code > is a name of the individual program, is the general output file containing input data, time and space steps, nonlinearity, energy and chemical potential, and was named fort.7 in the old Fortran version. The file < code > -den.txt is the output file with the condensate density, which had the names fort.3 and fort.4 in the old Fortran version for imaginary- and real-time propagation, respectively. Other density outputs, such as the initial density, are commented out to have a simpler set of output files. The users can re-introduce those by taking out the comment symbols, if needed. Table 1 Wall-clock execution times (in seconds) for runs with 1, 6 and 20 CPU cores with different programs using the Intel Fortran ifort (F-1, F-6 and F-20, respectively) and Intel C icc (C-1, C-6 and C-20, respectively) compilers using a workstation with two Intel Xeon E5-2650 v3 CPUs, with a total of 20 CPU cores, and obtained speedups (speedup-F=F-1/F-20, -speedupC=C-1/C-20 ) for 20 CPU cores. F-1 F-6 F-20 speedup-F C-1 C-6 C-20 speedup-C imag1d 32 26 26 1.2 45 28 27 1.7 imagcir 15 15 15 1.0 21 15 15 1.4 imagsph 12 12 12 1.0 19 12 10 1.9 real1d 194 84 72 2.7 304 110 98 3.1 realcir 132 62 57 2.3 182 78 64 2.8 realsph 119 68 67 1.8 191 76 61 3.1 imag2d 190 66 52 3.7 394 77 33 11.9 imagaxi 240 74 56 4.3 499 113 55 9.1 real2d 269 70 47 5.7 483 96 35 13.8 realaxi 132 37 25 5.3 237 51 22 10.8 imag3d 1682 472 366 4.6 2490 545 202 12.3 real3d 15,479 3494 2082 7.4 22,228 4558 1438 15.5 Fig. 1 (a) Speedup of the C and Fortran (F) imag3d programs as a function of the number of CPU cores, measured in a workstation with two Intel Xeon E5-2650 v3 CPUs. The speedup for the run with N CPU cores was calculated as the ratio between wall-clock execution times with one and N CPU cores. (b) Wall-clock time of the same runs as a function of the number of CPU cores. Also, some new output files are introduced in this version of programs. The files < code > -rms.txt are the output files with values of root-mean-square (rms) sizes in the multi-variable cases. There are new files with integrated densities, such as imag2d-den1d_x.txt, where the first part (imag2d) denotes that the density was calculated with the 2D program imag2d, and the second part (den1d_x) stands for the 1D density in the x -direction, obtained after integrating out the 2D density | ϕ ( x , y ) | 2 in the x – y plane over y -coordinate, (1) n 1 D ( x ) = ∫ − ∞ ∞ d y | ϕ ( x , y ) | 2 . Similarly, imag3d-den1d_x.txt and real3d-den1d_x.txt represent 1D densities from a 3D calculation obtained after integrating out the 3D density | ϕ ( x , y , z ) | 2 over y - and z -coordinate. The files imag3d-den2d_xy.txt and real3d-den2d_xy.txt are the integrated 2D densities in the x – y plane from a 3D calculation obtained after integrating out the 3D density over the z -coordinate, and similarly for other output files. Again, calculation and saving of these integrated densities is commented out in the programs, and can be activated by the user, if needed. In real-time propagation programs there are additional results for the dynamics saved in files, such as real2d-dyna.txt, where the first column denotes time, the second, third and fourth columns display rms sizes for the x -, y -, and r -coordinate, respectively. The dynamics is generated by multiplying the nonlinearity with a pre-defined factor during the NRUN iterations, and starting with the wave function calculated during the NPAS iterations. Such files were named fort.8 in the old Fortran versions of programs. There are similar files in the 3D real-time programs as well. Often it is needed to get a precise stationary state solution by imaginary-time propagation and then use it in the study of dynamics using real-time propagation. For that purpose, if the integer number NSTP is set to zero in real-time propagation, the density obtained in the imaginary-time simulation is used as initial wave function for real-time propagation, as in Refs. [4,5] . In addition, at the end of output files < code > -out.txt, we have introduced two new outputs, wall-clock execution time and CPU time for each run. We tested our programs on a workstation with two 10-core Intel Xeon E5-2650 v3 CPUs, and present results for all programs compiled with the Intel compiler. In Table 1 we show different wall-clock execution times for runs on 1, 6 and 20 CPU cores for Fortran and C. The corresponding columns “speedup-F” and “speedup-C” give the ratio of wall-clock execution times of runs on 1 and 20 CPU cores, and denote the actual measured speedup for 20 CPU cores. For the programs with effectively one spatial variable, the Fortran programs turn out to be quicker for small number of cores, whereas for larger number of CPU cores and for the programs with three spatial variables the C programs are faster. We also studied the speedup of the programs as a function of the number of available CPU cores. The performance for the imag3d Fortran and C programs is illustrated in Fig. 1 (a) and (b), where we plot the speedup and actual wall-clock time of the imag3d C and Fortran programs as a function of number of CPU cores in a workstation with two Intel Xeon E5-2650 v3 CPUs, with a total of 20 CPU cores. The plot in Fig. 1 (a) shows that the C program parallelizes more efficiently than the Fortran program. However, as the wall-clock time in Fortran for a single CPU core is less than that in C, the wall-clock times in both cases are comparable, viz. Fig. 1 (b). A saturation of the speedup with the increase of the number of CPU cores is expected in all cases. However, the saturation is attained quicker in Fortran than in C programs, and therefore the use of C programs could be recommended for larger number of CPU cores. For a small number of CPU cores the Fortran programs should be preferable. For example, from Table 1 we see that for 6 CPU cores the Fortran programs are faster than the C programs. In Fig. 1 (a) the saturation of the speedup of the Fortran program is achieved for approximately 10 CPU cores, when the wall-clock time of the C program crosses that of the Fortran program. Additional comments: This package consists of 24 programs, see Program title above. For the particular purpose of each program, please see descriptions below. Running time: Example inputs provided with the programs take less than 30 min in a workstation with two Intel Xeon Processors E5-2650 v3, 2 QPI links, 10 CPU cores (25 MB cache, 2.3 GHz). Program summary (i), (v), (vi), (vii), (xi), (xii) Program title: imag1d, imagcir, imagsph, real1d, realcir, realsph. Title of electronic files in C: (imag1d.c and imag1d.h), (imagcir.c and imagcir.h), (imagsph.c and imagsph.h), (real1d.c and real1d.h), (realcir.c and realcir.h), (realsph.c and realsph.h). Title of electronic files in Fortran 90: imag1d.f90, imagcir.f90, imagsph.f90, real1d.f90, realcir.f90, realsph.f90. Maximum RAM memory: 1 GB for the supplied programs. Programming language used: OpenMP C and Fortran 90. Typical running time: Minutes on a modern four-core PC. Nature of physical problem: These programs are designed to solve the time-dependent nonlinear partial differential GP equation in one spatial variable. Method of solution: The time-dependent GP equation is solved by the split-step Crank–Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation in imaginary time over small time steps. The method yields the solution of stationary problems. Program summary (ii), (iv), (viii), (x) Program title: imag2d, imagaxi, real2d, realaxi. Title of electronic files in C: (imag2d.c and imag2d.h), (imagaxi.c and imagaxi.h), (real2d.c and real2d.h), (realaxi.c and realaxi.h). Title of electronic files in Fortran 90: imag2d.f90, imagaxi.f90, real2d.f90, realaxi.f90. Maximum RAM memory: 1 GB for the supplied programs. Programming language used: OpenMP C and Fortran 90. Typical running time: Hour on a modern four-core PC. Nature of physical problem: These programs are designed to solve the time-dependent nonlinear partial differential GP equation in two spatial variables. Method of solution: The time-dependent GP equation is solved by the split-step Crank–Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation in imaginary time over small time steps. The method yields the solution of stationary problems. Program summary (iii), (ix) Program title: imag3d, real3d. Title of electronic files in C: (imag3d.c and imag3d.h), (real3d.c and real3d.h). Title of electronic files in Fortran 90: imag3d.f90, real3d.f90. Maximum RAM memory: 1 GB for the supplied programs. Programming language used: OpenMP C and Fortran 90. Typical running time: Few hours on a modern four-core PC. Nature of physical problem: These programs are designed to solve the time-dependent nonlinear partial differential GP equation in three spatial variables. Method of solution: The time-dependent GP equation is solved by the split-step Crank–Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation in imaginary time over small time steps. The method yields the solution of stationary problems. [ABSTRACT FROM AUTHOR]

We present hybrid OpenMP/MPI (Open Multi-Processing/Message Passing Interface) parallelized versions of earlier published C programs (Vudragović et al. 2012) for calculating both stationary and non-stationary solutions of the time-dependent Gross–Pitaevskii (GP) equation in three spatial dimensions. The GP equation describes the properties of dilute Bose–Einstein condensates at ultra-cold temperatures. Hybrid versions of programs use the same algorithms as the C ones, involving real- and imaginary-time propagation based on a split-step Crank–Nicolson method, but consider only a fully-anisotropic three-dimensional GP equation, where algorithmic complexity for large grid sizes necessitates parallelization in order to reduce execution time and/or memory requirements per node. Since distributed memory approach is required to address the latter, we combine MPI programming paradigm with existing OpenMP codes, thus creating fully flexible parallelism within a combined distributed/shared memory model, suitable for different modern computer architectures. The two presented C/OpenMP/MPI programs for real- and imaginary-time propagation are optimized and accompanied by a customizable makefile. We present typical scalability results for the provided OpenMP/MPI codes and demonstrate almost linear speedup until inter-process communication time starts to dominate over calculation time per iteration. Such a scalability study is necessary for large grid sizes in order to determine optimal number of MPI nodes and OpenMP threads per node. New version program summary Program title: GP-SCL-HYB package, consisting of: (i) imagtime3d-hyb, (ii) realtime3d-hyb. Catalogue identifier: AEDU_v3_0 Program Summary URL: http://cpc.cs.qub.ac.uk/summaries/AEDU_v3_0.html Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland. Licensing provisions: Apache License 2.0 No. of lines in distributed program, including test data, etc.: 26397. No. of bytes in distributed program, including test data, etc.: 161195. Distribution format: tar.gz. Programming language: C/OpenMP/MPI. Computer: Any modern computer with C language, OpenMP- and MPI-capable compiler installed. Operating system: Linux, Unix, Mac OS X, Windows. RAM: Total memory required to run programs with the supplied input files, distributed over the used MPI nodes: (i) 310 MB, (ii) 400 MB. Larger grid sizes require more memory, which scales with Nx*Ny*Nz. Number of processors used: No limit, from one to all available CPU cores can used on all MPI nodes. Number of nodes used: No limit on the number of MPI nodes that can be used. Depending on the grid size of the physical problem and communication overheads, optimal number of MPI nodes and threads per node can be determined by a scalability study for a given hardware platform. Classification: 2.9, 4.3, 4.12. Catalogue identifier of previous version: AEDU_v2_0. Journal reference of previous version: Comput. Phys. Commun. 183 (2012) 2021. Does the new version supersede the previous version?: No. Nature of problem: These programs are designed to solve the time-dependent Gross–Pitaevskii (GP) nonlinear partial differential equation in three spatial dimensions in a fully anisotropic trap using a hybrid OpenMP/MPI parallelization approach. The GP equation describes the properties of a dilute trapped Bose–Einstein condensate. Solution method: The time-dependent GP equation is solved by the split-step Crank–Nicolson method using discretization in space and time. The discretized equation is then solved by propagation, in either imaginary or real time, over small time steps. The method yields solutions of stationary and/or non-stationary problems. Reasons for the new version: Previous C [1] and Fortran [2] programs are widely used within the ultracold atoms and nonlinear optics communities, as well as in various other fields [3]. This new version represents extension of the two previously OpenMP-parallelized programs (imagtime3d-th and realtime3d-th) for propagation in imaginary and real time in three spatial dimensions to a hybrid, fully distributed OpenMP/MPI programs (imagtime3d-hyb and realtime3d-hyb). Hybrid extensions of previous OpenMP codes enable interested researchers to numerically study Bose–Einstein condensates in much greater detail (i.e., with much finer resolution) than with OpenMP codes. In OpenMP (threaded) versions of programs, numbers of discretization points in X , Y , and Z directions are bound by the total amount of available memory on a single computing node where the code is being executed. New, hybrid versions of programs are not limited in this way, as large numbers of grid points in each spatial direction can be evenly distributed among the nodes of a cluster, effectively distributing required memory over many MPI nodes. This is the first reason for development of hybrid versions of 3d codes. The second reason for new versions is speedup in the execution of numerical simulations that can be gained by using multiple computing nodes with OpenMP/MPI codes. Summary of revisions: Two C/OpenMP programs in three spatial dimensions from previous version [1] of the codes (imagtime3d-th and realtime3d-th) are transformed and rewritten into a hybrid OpenMP/MPI programs and named imagtime3d-hyb and realtime3d-hyb. The overall structure of two programs is identical. The directory structure of the GP-SCL-HYB package is extended compared to the previous version and now contains a folder scripts, where examples of scripts that can be used to run the programs on a typical MPI cluster are given. The corresponding readme.txt file contains more details. We have also included a makefile with tested and verified settings for most popular MPI compliers, including OpenMPI (Open Message Passing Interface) [4] and MPICH (Message Passing Interface Chameleon) [5]. Transformation from pure OpenMP to a hybrid OpenMP/MPI approach has required that the array containing condensate wavefunction is distributed among MPI nodes of a computer cluster. Several data distribution models have been considered for this purpose, including block distribution and block cyclic distribution of data in a 2d matrix. Finally, we decided to distribute the wavefunction values across different nodes so that each node contains only one slice of the X -dimension data, while containing the complete corresponding Y - and Z -dimension data, as illustrated in Fig. 1 . This allows central functions of our numerical algorithm, calcluy, calcuz, and calcnu to be executed purely in parallel on different MPI nodes of a cluster, without any overhead or communication, as nodes contain all the information for Y - and Z -dimension data in the given X -sub-domain. However, the problem arises when functions calclux, calcrms, and calcmuen need to be executed, as they also operate on the whole X -dimension data. Thus, the need for additional communication arises during the execution of the function calcrms, while in the case of functions calclux and calcmuen also the transposition of data between X - and Y -dimensions is necessary, while data in Z dimension have to stay contiguous. Transposition provides nodes with all the necessary X -dimension data to execute functions calclux and calcmuen. However, this needs to be done in each iteration of numerical algorithm, thus necessarily increasing communication overhead of the simulation. Transposition algorithms that were considered where the ones that account for greatest common divisor (GCD) between number of nodes in columns (designated by N ) and rows (designated by M) of a cluster configured as 2d mash of nodes [6]. Two of such algorithms have been tested and tried for implementation: the case when GCD=1 and the case when GCD > 1 . The trivial situation N = M = 1 is already covered by the previous, purely OpenMP programs, and therefore, without any loss of generality, we have considered only configurations with number of nodes in X -dimension satisfying N > 1 . Only the former algorithm (GCD=1) was found to be sound in case where data matrix is not a 2d, but a 3d structure. Latter case was found to be too demanding implementation-wise, since MPI functions and data-types are bound to certain limitations. Therefore, the algorithm with M = 1 nodes in Y -dimension was implemented, as depicted by the wavefunction data structure in Fig. 1 . Fig. 1 BEC wavefunction data structure in the employed algorithm. Data are sliced so that the complete Y - and Z -dimension data reside on a single node for a given range of data in X -dimension, while data in X -dimension are distributed over N = gsize nodes . Figures show data transposition and MPI indexed datatype creation parameters for the case of: (a) sending side and (b) receiving side. Implementation of the algorithm relies on a sliced distribution of data among the nodes, as explained in Fig. 2 . This successfully solves the problem of large RAM consumption of 3d codes, which arises even for moderate grid sizes. However, it does not solve the question of data transposition between the nodes. In order to implement the most effective (GCD=1) transposition algorithm according to Ref. [6], we had to carry out block distribution of data within one data slice contained on a single node. This block distribution of data was done implicitly, i.e. data on one node have been put in a single 1d array (psi) of contiguous memory, in which Z -dimension has stride 1, Y -dimension has stride Nz, and X -dimension has stride Ny*Nz. This is different from previous implementation of the programs, where the wavefunction was represented by an explicit 3d array. This change was also introduced in order to more easily form user MPI datatypes, which allow for implicit block distribution of data, and represent 3d blocks of data within 1d data array. These blocks are then swapped between nodes, effectively performing the transposition in X – Y and Y – X directions. Together with transposition of blocks between the nodes, the block data also have to be redistributed. To illustrate how this works, let us consider example shown in Fig. 1 (a), where one data block has size (Nx/gisze)*(Ny/gsize)*Nz. It represents one 3d data block, swapped between two nodes of a cluster (through one non-blocking MPI_Isend and one MPI_Ireceive operation), containing (Nx/gsize)*(Ny/gsize) 1d rods of contiguous Nz data. These rods themselves need to be transposed within the transposed block as well. This means that two levels of transpositions need to be performed. At a single block level, rods have to be transposed (as indicated in upper left corner of Fig. 1 (a) for sending index type and in Fig. 1 (b) for receiving index type). Second level is transposition of blocks between different nodes, which is depicted by blue arrows connecting different blocks in Fig. 1 . The above described transposition is applied whenever needed in the functions calclux and calcmuen, which require calculations to be done on the whole range of data in X -dimension. When performing renormalization of the wavefunction or calculation of its norm, root-mean-square radius, chemical potential, and energy, collective operations MPI_Gather and MPI_Bcast are also used. Figs. 3 and 4 show the scalability results obtained for hybrid versions of programs for small and large grid sizes as a function of number of MPI nodes used. The baseline for calculation of speedups in the execution time for small grid sizes is previous, purely OpenMP programs, while for large grid sizes, which cannot fit onto a single node, the baseline is hybrid programs with minimal configuration runs on 8 nodes. The figures also show efficacies, defined as percentages of measured speedups compared to the ideal ones. We see that an excellent scalability (larger than 80% compared to the ideal one) can be obtained for up to 32 nodes. The tests have been performed on a cluster with nodes containing 2 × 8 -core Sandy Bridge Xeon 2.6 GHz processors with 32 GB of RAM and Infiniband QDR (Quad Data Rate, 40 Gbps) interconnect. We stress that the scalability depends greatly on the ratio between the calculation and communication time per iteration, and has to be studied for a particular type of processors and interconnect technology. Additional comments: This package consists of 2 programs, see Program title above. Both are hybrid, threaded and distributed (OpenMP/MPI parallelized). For the particular purpose of each program, see descriptions below. Fig. 2 Creation of a user-defined MPI datatype indextype with the function MPI_Type_indexed. Here, count represents the number of blocks, blocklengths array contains lengths of each block, and displacements array contains the displacement of each block from the beginning of the corresponding data structure. For example, if an array of double precision numbers (designated as buffer in the figure) is sent by MPI_Send with the datatype set to indextype, it is interpreted as a block-distributed data structure, as specified when indextype was created. Fig. 3 Speedup in the execution time and efficacy curves of imagtime3d-hyb and realtime3d-hyb programs as a function of the number of MPI nodes used for small grid sizes. The results are obtained on a cluster with nodes containing 2 × 8 -core Sandy Bridge Xeon 2.6 GHz processors with 32 GB of RAM and Infiniband QDR interconnect: (a) speedup of imagtime3d-hyb on a 240 × 200 × 160 grid; (b) efficacy of imagtime3d-hyb on a 240 × 200 × 160 grid; (c) speedup of realtime3d-hyb on a 200 × 160 × 120 grid; (d) efficacy of realtime3d-hyb on a 200 × 160 × 120 grid. Shaded areas in graphs (b) and (d) represent high-efficacy regions, where speedup is at least 80% of the ideal one. Running time: All running times given in descriptions below refer to programs compiled with OpenMPI/GCC compiler and executed on 8–32 nodes with 2 × 8 -core Sandy Bridge Xeon 2.6 GHz processors with 32 GB of RAM and Infiniband QDR interconnect. With the supplied input files for small grid sizes, running wallclock times of several minutes are required on 8–10 MPI nodes. Special features: (1) Since the condensate wavefunction data are distributed among the MPI nodes, when writing wavefunction output files each MPI process saves its data into a separate file, to avoid I/O issues. Concatenating the corresponding files from all MPI processes will create the complete wavefunction file. (2) Due to a known bug in OpenMPI up to version 1.8.4, allocation of memory for indexed datatype on a single node for large grids (such as 800 × 640 × 480 ) may fail. The fix for this bug is already in 3c489ea branch and is fixed in OpenMPI as of version 1.8.5. Fig. 4 Speedup in the execution time and efficacy curves of imagtime3d-hyb and realtime3d-hyb programs as a function of the number of MPI nodes used for large grid sizes. The results are obtained on a cluster with nodes containing 2 × 8 -core Sandy Bridge Xeon 2.6 GHz processors with 32 GB of RAM and Infiniband QDR interconnect: (a) speedup of imagtime3d-hyb on a 1920 × 1600 × 1280 grid; (b) efficacy of imagtime3d-hyb on a 1920 × 1600 × 1280 grid; (c) speedup of realtime3d-hyb on a 1600 × 1280 × 960 grid; (d) efficacy of realtime3d-hyb on a 1600 × 1280 × 960 grid. Shaded areas in graphs (b) and (d) represent high-efficacy regions, where speedup is at least 80% of the ideal one. Program summary (i) Program title: imagtime3d-hyb. Title of electronic files: imagtime3d-hyb.c, imagtime3d-hyb.h. Computer: Any modern computer with C language, OpenMP- and MPI-capable compiler installed. RAM memory requirements: 300 MBytes of RAM for a small grid size 240 × 200 × 160 , and scales with Nx*Ny*Nz. This is total amount of memory needed, and is distributed over MPI nodes used for execution. Programming language used: C/OpenMP/MPI. Typical running time: Few minutes with the supplied input files for a small grid size 240 × 200 × 160 on 8 nodes. Up to one hour for a large grid size 1920 × 1600 × 1280 on 32 nodes (1000 iterations). Nature of physical problem: This program is designed to solve the time-dependent GP nonlinear partial differential equation in three space dimensions with an anisotropic trap. The GP equation describes the properties of a dilute trapped Bose–Einstein condensate. Method of solution: The time-dependent GP equation is solved by the split-step Crank–Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation in imaginary time over small time steps. The method yields solutions of stationary problems. Program summary (ii) Program title: realtime3d-hyb. Title of electronic files: realtime3d-hyb.c, realtime3d-hyb.h. Computer: Any modern computer with C language, OpenMP- and MPI-capable compiler installed. RAM memory requirements: 410 MBytes of RAM for a small grid size 200 × 160 × 120 , and scales with Nx*Ny*Nz. This is total amount of memory needed, and is distributed over MPI nodes used for execution. Programming language used: C/OpenMP/MPI. Typical running time: 10–15 min with the supplied input files for a small grid size 200 × 160 × 120 on 10 nodes. Up to one hour for a large grid size 1600 × 1280 × 960 on 32 nodes (1000 iterations). Nature of physical problem: This program is designed to solve the time-dependent GP nonlinear partial differential equation in three space dimensions with an anisotropic trap. The GP equation describes the properties of a dilute trapped Bose–Einstein condensate. Method of solution: The time-dependent GP equation is solved by the split-step Crank–Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation in real time over small time steps. The method yields solutions of stationary and non-stationary problems. Acknowledgments B.S., V.S., A.B., and A.B. acknowledge support by the Ministry of Education, Science, and Technological Development of the Republic of Serbia under projects ON171017, III43007, ON171009, ON174027 and IBEC, and by DAAD - German Academic and Exchange Service under project IBEC. P.M. acknowledges support by the Science and Engineering Research Board, Department of Science and Technology, Government of India under project No. EMR/2014/000644. S.K.A. acknowledges support by the CNPq of Brazil under project 303280/2014-0, and by the FAPESP of Brazil under project 2012/00451-0. Numerical simulations were run on the PARADOX supercomputing facility at the Scientific Computing Laboratory of the Institute of Physics Belgrade, supported in part by the Ministry of Education, Science, and Technological Development of the Republic of Serbia under project ON171017. References [1] D. Vudragović, I. Vidanović, A. Balaž, P. Muruganandam, S. K. Adhikari, C programs for solving the time-dependent Gross–Pitaevskii equation in a fully anisotropic trap, Comput. Phys. Commun. 183 (2012) 2021. [2] P. Muruganandam and S. K. Adhikari, Fortran programs for the time-dependent Gross–Pitaevskii equation in a fully anisotropic trap, Comput. Phys. Commun. 180 (2009) 1888. [3] R. K. Kumar and P. Muruganandam, J. Phys. B: At. Mol. Opt. Phys. 45 (2012) 215301; L. E. Young-S. and S. K. Adhikari, Phys. Rev. A 86 (2012) 063611; S. K. Adhikari, J. Phys. B: At. Mol. Opt. Phys. 45 (2012) 235303; I. Vidanović, N. J. van Druten, and M. Haque, New J. Phys. 15 (2013) 035008; S. Balasubramanian, R. Ramaswamy, and A. I. Nicolin, Rom. Rep. Phys. 65 (2013) 820; L. E. Young-S. and S. K. Adhikari, Phys. Rev. A 87 (2013) 013618; H. Al-Jibbouri, I. Vidanovic, A. Balaz, and A. Pelster, J. Phys. B: At. Mol. Opt. Phys. 46 (2013) 065303; X. Antoine, W. Bao, and C. Besse, Comput. Phys. Commun. 184 (2013) 2621; B. Nikolić, A. Balaž, and A. Pelster, Phys. Rev. A 88 (2013) 013624; H. Al-Jibbouri and A. Pelster, Phys. Rev. A 88 (2013) 033621; S. K. Adhikari, Phys. Rev. A 88 (2013) 043603; J. B. Sudharsan, R. Radha, and P. Muruganandam, J. Phys. B: At. Mol. Opt. Phys. 46 (2013) 155302; R. R. Sakhel, A. R. Sakhel, and H. B. Ghassib, J. Low Temp. Phys. 173 (2013) 177; E. J. M. Madarassy and V. T. Toth, Comput. Phys. Commun. 184 (2013) 1339; R. K. Kumar, P. Muruganandam, and B. A. Malomed, J. Phys. B: At. Mol. Opt. Phys. 46 (2013) 175302; W. Bao, Q. Tang, and Z. Xu, J. Comput. Phys. 235 (2013) 423; A. I. Nicolin, Proc. Rom. Acad. Ser. A-Math. Phys. 14 (2013) 35; R. M. Caplan, Comput. Phys. Commun. 184 (2013) 1250; S. K. Adhikari, J. Phys. B: At. Mol. Opt. Phys. 46 (2013) 115301; Ž. Marojević, E. Göklü, and C. Lämmerzahl, Comput. Phys. Commun. 184 (2013) 1920; X. Antoine and R. Duboscq, Comput. Phys. Commun. 185 (2014) 2969; S. K. Adhikari and L. E. Young-S, J. Phys. B: At. Mol. Opt. Phys. 47 (2014) 015302; K. Manikandan, P. Muruganandam, M. Senthilvelan, and M. Lakshmanan, Phys. Rev. E 90 (2014) 062905; S. K. Adhikari, Phys. Rev. A 90 (2014) 055601; A. Balaž, R. Paun, A. I. Nicolin, S. Balasubramanian, and R. Ramaswamy, Phys. Rev. A 89 (2014) 023609; S. K. Adhikari, Phys. Rev. A 89 (2014) 013630; J. Luo, Commun. Nonlinear Sci. Numer. Simul. 19 (2014) 3591; S. K. Adhikari, Phys. Rev. A 89 (2014) 043609; K.-T. Xi, J. Li, and D.-N. Shi, Physica B 436 (2014) 149; M. C. Raportaru, J. Jovanovski, B. Jakimovski, D. Jakimovski, and A. Mishev, Rom. J. Phys. 59 (2014) 677; S. Gautam and S. K. Adhikari, Phys. Rev. A 90 (2014) 043619; A. I. Nicolin, A. Balaž, J. B. Sudharsan, and R. Radha, Rom. J. Phys. 59 (2014) 204; K. Sakkaravarthi, T. Kanna, M. Vijayajayanthi, and M. Lakshmanan, Phys. Rev. E 90 (2014) 052912; S. K. Adhikari, J. Phys. B: At. Mol. Opt. Phys. 47 (2014) 225304; R. K. Kumar and P. Muruganandam, Numerical studies on vortices in rotating dipolar Bose–Einstein condensates, Proceedings of the 22nd International Laser Physics Workshop, J. Phys. Conf. Ser. 497 (2014) 012036; A. I. Nicolin and I. Rata, Density waves in dipolar Bose–Einstein condensates by means of symbolic computations, High-Performance Computing Infrastructure for South East Europe’s Research Communities: Results of the HP-SEE User Forum 2012, in Springer Series: Modeling and Optimization in Science and Technologies 2 (2014) 15; S. K. Adhikari, Phys. Rev. A 89 (2014) 043615; R. K. Kumar and P. Muruganandam, Eur. Phys. J. D 68 (2014) 289; J. B. Sudharsan, R. Radha, H. Fabrelli, A. Gammal, and B. A. Malomed, Phys. Rev. A 92 (2015) 053601; S. K. Adhikari, J. Phys. B: At. Mol. Opt. Phys. 48 (2015) 165303; F. I. Moxley III, T. Byrnes, B. Ma, Y. Yan, and W. Dai, J. Comput. Phys. 282 (2015) 303; S. K. Adhikari, Phys. Rev. E 92 (2015) 042926; R. R. Sakhel, A. R. Sakhel, and H. B. Ghassib, Physica B 478 (2015) 68; S. Gautam and S. K. Adhikari, Phys. Rev. A 92 (2015) 023616; D. Novoa, D. Tommasini, and J. A. Nóvoa-López, Phys. Rev. E 91 (2015) 012904; S. Gautam and S. K. Adhikari, Laser Phys. Lett. 12 (2015) 045501; K.-T. Xi, J. Li, and D.-N. Shi, Physica B 459 (2015) 6; R. K. Kumar, L. E. Young-S., D. Vudragović, A. Balaž, P. Muruganandam, and S. K. Adhikari, Comput. Phys. Commun. 195 (2015) 117; S. Gautam and S. K. Adhikari, Phys. Rev. A 91 (2015) 013624; A. I. Nicolin, M. C. Raportaru, and A. Balaž, Rom. Rep. Phys. 67 (2015) 143; S. Gautam and S. K. Adhikari, Phys. Rev. A 91 (2015) 063617; E. J. M. Madarassy and V. T. Toth, Phys. Rev. D 91 (2015) 044041. [4] Open Message Passing Interface (OpenMPI), http://www.open-mpi.org/ (2015). [5] Message Passing Interface Chameleon (MPICH), https://www.mpich.org/ (2015). [6] J. Choi, J. J. Dongarra, D. W. Walker, Parallel matrix transpose algorithms on distributed memory concurrent computers, Parallel Comput. 21 (1995) 1387. [ABSTRACT FROM AUTHOR]

A parallel memetic global optimization algorithm suitable for shared memory multicore systems is proposed and analyzed. The considered algorithm combines two well-known and widely used population-based stochastic algorithms, namely Particle Swarm Optimization and Differential Evolution, with two efficient and parallelizable local search procedures. The sequential version of the algorithm was first introduced as MEMPSODE (MEMetic Particle Swarm Optimization and Differential Evolution) and published in the CPC program library. We exploit the inherent and highly irregular parallelism of the memetic global optimization algorithm by means of a dynamic and multilevel approach based on the OpenMP tasking model. In our case, tasks correspond to local optimization procedures or simple function evaluations. Parallelization occurs at each iteration step of the memetic algorithm without affecting its searching efficiency. The proposed implementation, for the same random seed, reaches the same solution irrespectively of being executed sequentially or in parallel. Extensive experimental evaluation has been performed in order to illustrate the speedup achieved on a shared-memory multicore server. Program summary Program title: p-MEMPSODE Catalogue identifier: AEXJ_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEXJ_v1_0.html Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 9950 No. of bytes in distributed program, including test data, etc.: 141503 Distribution format: tar.gz Programming language: ANSI C. Computer: Workstation. Operating system: Developed under the Linux operating system using the GNU compilers v.4.4.3 (or higher). Uses the OpenMP API and runtime system. RAM: The code uses O ( n × N ) internal storage, n being the dimension of the problem and N the maximum population size. The required memory is dynamically allocated. Word size: 64 Classification: 4.9. Nature of problem: Numerical global optimization of real valued functions is an indispensable methodology for solving a multitude of problems in science and engineering. Many problems exhibit a number of local and/or global minimizers, expensive function evaluations or require real-time response. In addition, discontinuities of the objective function, non-smooth and deceitful landscapes constitute challenging obstacles for most optimization algorithms. Solution method: We implement a memetic global optimization algorithm that combines stochastic, population-based methods with deterministic local search procedures. More specifically, the Unified Particle Swarm Optimization and the Differential Evolution algorithms are harnessed with the derivative-free Torczon’s Multi-Directional Search and the gradient-based BFGS method. The produced hybrid algorithms possess inherent parallelism that is exploited efficiently by means of the OpenMP tasking model. Given the same random seed, the proposed implementation reaches the same solution irrespective of being executed sequentially or in parallel. Restrictions: The current version of the software uses only double precision arithmetic. An OpenMP-enabled (version 3.0 or higher) compiler is required. Unusual features: The software requires bound constraints on the optimization variables. Running time: The running time depends on the complexity of the objective function (and its derivatives if used) as well as on the number of available cores. Extensive experimental results demonstrate that the speedup closely approximates ideal values. [ABSTRACT FROM AUTHOR]

Particle tracking in large-scale numerical simulations of turbulent flows presents one of the major bottlenecks in parallel performance and scaling efficiency. Here, we describe a particle tracking algorithm for large-scale parallel pseudo-spectral simulations of turbulence which scales well up to billions of tracer particles on modern high-performance computing architectures. We summarize the standard parallel methods used to solve the fluid equations in our hybrid MPI/OpenMP implementation. As the main focus, we describe the implementation of the particle tracking algorithm and document its computational performance. To address the extensive inter-process communication required by particle tracking, we introduce a task-based approach to overlap point-to-point communications with computations, thereby enabling improved resource utilization. We characterize the computational cost as a function of the number of particles tracked and compare it with the flow field computation, showing that the cost of particle tracking is very small for typical applications. [ABSTRACT FROM AUTHOR]

The computational performance of a smoothed particle hydrodynamics (SPH) simulation is investigated for three types of current shared-memory parallel computer devices: many integrated core (MIC) processors, graphics processing units (GPUs), and multi-core CPUs. We are especially interested in efficient shared-memory allocation methods for each chipset, because the efficient data access patterns differ between compute unified device architecture (CUDA) programming for GPUs and OpenMP programming for MIC processors and multi-core CPUs. We first introduce several parallel implementation techniques for the SPH code, and then examine these on our target computer architectures to determine the most effective algorithms for each processor unit. In addition, we evaluate the effective computing performance and power efficiency of the SPH simulation on each architecture, as these are critical metrics for overall performance in a multi-device environment. In our benchmark test, the GPU is found to produce the best arithmetic performance as a standalone device unit, and gives the most efficient power consumption. The multi-core CPU obtains the most effective computing performance. The computational speed of the MIC processor on Xeon Phi approached that of two Xeon CPUs. This indicates that using MICs is an attractive choice for existing SPH codes on multi-core CPUs parallelized by OpenMP, as it gains computational acceleration without the need for significant changes to the source code. [ABSTRACT FROM AUTHOR]

The prospect of wildly heterogeneous computer systems has led to a renewed discussion of programming approaches in high-performance computing, of which computational fluid dynamics is a major field. The challenge consists in harvesting the performance of all available hardware components while retaining good programmability. In particular the use of graphic cards is an important trend. This is addressed in the present paper by devising a hybrid programming model to create a heterogeneous data-parallel computation with a single source code. The concept is demonstrated for a one-dimensional spectral-element discretization of a fluid dynamics problem. To exploit the additional hardware available when coupling GPGPU-accelerated processes with excess CPU cores, a straight-forward load balancing model for such heterogeneous environments is developed. The paper presents a large number of run time measurements and demonstrates that the achieved performance gains are close to optimal. This provides valuable information for the implementation of fluid dynamics codes on modern heterogeneous hardware. [ABSTRACT FROM AUTHOR]

We present a cross-platform scaling investigation for an OpenMP parallelization of UTBEST3D – a coastal and regional ocean code based on the discontinuous Galerkin finite element method. The study is conducted for a real life application on an unstructured computational mesh of the Northwest Atlantic with realistic topography and well resolved coast line on a broad selection of current computing platforms. Four numerical setups of increasing physical and computational complexity are used for comparison: barotropic with no vertical eddy viscosity, barotropic with an algebraic eddy viscosity parametrization, baroclinic with an algebraic eddy viscosity, and baroclinic with k – ∊ vertical turbulence closure. In addition to Intel Xeon and IBM Power6/PowerPC architectures, we also include Intel’s new MIC processor Xeon Phi in the evaluation. Good scalability is found across all platforms with Intel Xeon CPUs producing the best runtime results and Xeon Phi demonstrating the best parallel efficiency. [ABSTRACT FROM AUTHOR]

The trend of all modern computer architectures, and the path to exascale, is towards increasing numbers of lower power cores, with a decreasing memory to core ratio. This imposes a strong evolutionary pressure on algorithms and software to efficiently utilise all levels of parallelism available on a given platform while minimising data movement. Unstructured finite elements codes have long been effectively parallelised using domain decomposition methods, implemented using libraries such as the Message Passing Interface (MPI). However, there are many optimisation opportunities when threading is used for intra-node parallelisation for the latest multi-core/many-core platforms. The benefits include increased algorithmic freedom, reduced memory requirements, cache sharing, reduced number of partitions, less MPI communication and I/O overhead. In this paper, we report progress in implementing a hybrid OpenMP–MPI version of the unstructured finite element code Fluidity. For matrix assembly kernels, the OpenMP parallel algorithm uses graph colouring to identify independent sets of elements that can be assembled concurrently with no race conditions. In this phase there are no MPI overheads as each MPI process only assembles its own local part of the global matrix. We use an OpenMP threaded fork of PETSc to solve the resulting sparse linear systems of equations. We experiment with a range of preconditioners, including HYPRE which provides the algebraic multigrid preconditioner BoomerAMG where the smoother is also threaded. Since unstructured finite element codes are well known to be memory latency bound, particular attention is paid to ccNUMA architectures where data locality is particularly important to achieve good intra-node scaling characteristics. We also demonstrate that utilising non-blocking algorithms and libraries are critical to mixed-mode application so that it can achieve better parallel performance than the pure MPI version. [ABSTRACT FROM AUTHOR]

Nowadays, an increasing number of researchers have demonstrated a tendency to choose hybrid CPU–GPU hybrid computing as a high performance computing alternative. Entropic Lattice Boltzmann method (ELBM) parallelization, like many parallel algorithms in the field of rapid scientific and engineering computing, has given rise to much attention for applications of computational fluid dynamics. This study aims to present an efficient implementation of ELBM flow simulation for the D3Q19 model in a hybrid CPU–GPU computing environment, which consists of AMD multi-core CPUs with NVIDIA Graphics Processing Units (GPUs). To overcome the GPU memory size limitation and communication overhead, we propose a set of techniques for the development of an efficient ELBM algorithm for hybrid CPU–GPU computation. Considering the contribution of computational capacity for both the CPU and GPU, an efficient load balancing model is built. The efficiency and accuracy of the proposed approach and established model are tested on a hybrid CPU–GPU accelerated system, where the intensive parts of the computation are dealt with the software framework OpenMP and CUDA. Finally, we show the comparison of resulting computational performance using a hybrid CPU–GPU approach against both a single CPU core and a single GPU device. [ABSTRACT FROM AUTHOR]

Presented is the application of the Split Harten, Lax and van Leer (SHLL) technique applied to parallel computation using a hybrid OpenMP/AVX parallelization paradigm for the Shallow Water Equations and Euler Equations. The key behind the ease of parallelization of the SHLL method for both governing equations is the mathematical/vector splitting in each coordinate direction – this splitting results in a high degree of locality, producing a scheme which is embarrassingly parallel and well suited for the vectorization capacities offered by vector-computing architectures. Here we demonstrate this capacity using the SIMD capacities of modern CPUs, namely the Advanced Vector eXtensions (AVX) capability of recent CPUs. The main feature of AVX is the capacity to perform SIMD operations on 8 floating point variables in parallel – an increase from 4 floating point variables as possible using the previous SIMD Streaming Extensions (SSE). Furthermore, since modern CPU’s employ a large number of cores, we further extend the performance by using AVX on each available CPU core using shared memory (OpenMP) parallelization. We present a direction-split higher order extension to both the SHLL method and apply it to AVX through the use of intrinsic functions in the flux computation and state computation modules. High performance is obtained by ensuring that all flux computations are performed using only AVX intrinsic functions – no computations are performed in serial. Through this approach, a single workstation with 2× Xeon CPU’s (16 physical cores) allows a performance increase of over 117 times that of a single core alone in the flux evaluation kernel. [ABSTRACT FROM AUTHOR]

OFF, an open source (free software) code for performing fluid dynamics simulations, is presented. The aim of OFF is to solve, numerically, the unsteady (and steady) compressible Navier–Stokes equations of fluid dynamics by means of finite volume techniques: the research background is mainly focused on high-order (WENO) schemes for multi-fluids, multi-phase flows over complex geometries. To this purpose a highly modular, object-oriented application program interface (API) has been developed. In particular, the concepts of data encapsulation and inheritance available within Fortran language (from standard 2003) have been stressed in order to represent each fluid dynamics “entity” (e.g. the conservative variables of a finite volume, its geometry, etc…) by a single object so that a large variety of computational libraries can be easily (and efficiently) developed upon these objects. The main features of OFF can be summarized as follows: Programming Language OFF is written in standard (compliant) Fortran 2003; its design is highly modular in order to enhance simplicity of use and maintenance without compromising the efficiency; Parallel Frameworks Supported the development of OFF has been also targeted to maximize the computational efficiency: the code is designed to run on shared-memory multi-cores workstations and distributed-memory clusters of shared-memory nodes (supercomputers); the code’s parallelization is based on Open Multiprocessing (OpenMP) and Message Passing Interface (MPI) paradigms; Usability, Maintenance and Enhancement in order to improve the usability, maintenance and enhancement of the code also the documentation has been carefully taken into account; the documentation is built upon comprehensive comments placed directly into the source files (no external documentation files needed): these comments are parsed by means of doxygen free software producing high quality html and latex documentation pages; the distributed versioning system referred as git has been adopted in order to facilitate the collaborative maintenance and improvement of the code; Copyrights OFF is a free software that anyone can use, copy, distribute, study, change and improve under the GNU Public License version 3. The present paper is a manifesto of OFF code and presents the currently implemented features and ongoing developments. This work is focused on the computational techniques adopted and a detailed description of the main API characteristics is reported. OFF capabilities are demonstrated by means of one and two dimensional examples and a three dimensional real application. Program summary: Program title: OFF Catalogue identifier: AESV_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AESV_v1_0.html Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland Licensing provisions: GNU General Public Licence, version 3 No. of lines in distributed program, including test data, etc.: 60466 No. of bytes in distributed program, including test data, etc.: 595575 Distribution format: tar.gz Programming language: Fortran (standard 2003 or newer); developed and tested with Intel Fortran Compiler v. 12.x or newer. Computer: Designed for shared-memory multi-cores workstations and for hybrid distributed/shared-memory supercomputers, but any computer system with a Fortran (2003+) compiler is suited. Operating system: Designed for POSIX architecture and tested on GNU/Linux one. Has the code been vectorized or parallelized?: Hybrid parallelization by means of MPI library and OpenMP paradigm, tested on up to 256 processors. RAM: [1 MB; 1 GB] x core, simulation-dependent Classification: 4.3, 4.10, 12. External routines: The proprietary library [1] must be linked for producing binary outputs in Tecplot Inc. format; the MPI library [2] must be linked for running on distribute-memory parallel systems Nature of problem: Numerical solution of the Compressible Navier–Stokes equations for multi-fluids multi-phase flows in complex geometries Solution method: Fully-conservative Finite Volume scheme based on very high-order WENO Positivity-Preserving reconstruction technique and Strong Stability Preserving high-order Runge–Kutta time integration. Structured multi-block general curvilinear and body-fitted grids constitute the underlining numerical grids. MPI and OpenMP paradigms are used for parallel computations. Pre-processing tool to deal with commercial meshing softwares is provided as well as post-processing one for numerical results visualization Restrictions: At present, OFF is validated for simulating inviscid and single-phase flows (viscous fluxes computation and fully coupled Eulerian /Lagrangian schemes have still to be validated, but they are already implemented); modern Fortran compiler is mandatory Unusual features: OFF is a complex CFD software strongly based on Object-Oriented Programming paradigm by means of modern Fortran standard (2003 or higher) Additional comments: OFF project adopts Git [3], a free and open source distributed version control system. A public repository dedicated to OFF project [4] has been created on github, a web-based hosting service for software development projects using git versioning system. Finally, a comprehensive documentation [5] is provided parsing source code comments by means of doxygen software [6] Running time: The running time depends on the available computation resources, the complexity of the problem and the selected accuracy of the numerical scheme. Simple one dimensional problems require few minutes on personal workstation-like systems whereas complex three dimensional problems needing high accuracy (e.g. DNS) could require weeks on supercomputers. OFF has proven to have a good scalability on small clusters (e.g. CASPUR facilities, namely on Matrix, a GNU/Linux cluster composed by 320 nodes each one constituted by a dual Opteron quadcore at 2.1GHz with 16/32 GB of RAM) References: [[1]] TecIO Library, Tecplot Inc. proprietary library for I/O binary files in Tecplot format, http://www.tecplot.com/downloads/tecio-library/. [[2]] The Message Passing Interface (MPI) standard, a library specification for message-passing, proposed as a standard by a broadly based committee of vendors, implementors, and users, http://www.mcs.anl.gov/research/projects/mpi/. [[3]] Git, a free and open source distributed version control system, http://git-scm.com/. [[4]] Github, a web-based hosting service for software development projects using git versioning system, https://github.com. [[5]] Ocial OFF documentation, http://szaghi.github.com/OFF/index.html. [[6]] Doxygen, a documentation system for many programming languages, http://www.stack.nl/dimitri/doxygen. [ABSTRACT FROM AUTHOR]

Abstract: Multiphase flow implementations of the lattice Boltzmann method (LBM) are widely applied to the study of porous medium systems. In this work, we construct a new variant of the popular “color” LBM for two-phase flow in which a three-dimensional, 19-velocity (D3Q19) lattice is used to compute the momentum transport solution while a three-dimensional, seven velocity (D3Q7) lattice is used to compute the mass transport solution. Based on this formulation, we implement a novel heterogeneous GPU-accelerated algorithm in which the mass transport solution is computed by multiple shared memory CPU cores programmed using OpenMP while a concurrent solution of the momentum transport is performed using a GPU. The heterogeneous solution is demonstrated to provide speedup of as compared to multi-core CPU solution and compared to GPU solution due to concurrent utilization of both CPU and GPU bandwidths. Furthermore, we verify that the proposed formulation provides an accurate physical representation of multiphase flow processes and demonstrate that the approach can be applied to perform heterogeneous simulations of two-phase flow in porous media using a typical GPU-accelerated workstation. [Copyright &y& Elsevier]

Abstract: A Graphics Processing Unit (GPU)-CUDA C and (Multi-core)-OpenMP versions of the Reaction Ensemble Monte Carlo method (REMC) are presented. The REMC algorithm is a powerful tool to investigate the equilibrium behavior of chemically reacting systems in highly non-ideal conditions. Both the GPU and the Multi-core versions of the code are particularly efficient when the total potential energy of the system must be calculated, as in the constant-pressure systems. Results, obtained in the case of Helium plasma at high pressure, show differences between real and ideal cases. [Copyright &y& Elsevier]

Abstract: In this paper, development of a new moving computational domain method and its parallel implementation are presented. The moving computational domain method is a method where the whole of the computational domain including bodies inside moves in the physical space without the limit of region size. The method is developed for simulation of falling spheres in a long pipe. We investigate flows around multiple bodies that move and interact with each other in accordance with six degrees of freedom motion. Movement of the falling multiple spheres in an infinite long square pipe is simulated. The parallel implementation based on the Domino method, which is a kind of block hyper-plane parallelism, is applied to the time-consuming Poisson solver. The parallel computation is performed using a multi-core processor with up to 24 cores, and its performance is tested. [Copyright &y& Elsevier]

Abstract: The work is devoted to the development of efficient parallel algorithms for large-scale simulations of incompressible flows on hybrid supercomputers based on massively-parallel accelerators. The governing equations are discretized using a high-order finite-volume scheme for Cartesian staggered meshes with the only restriction that, at least, one direction is periodic. Its “classical” MPI+OpenMP parallel implementation for CPUs was designed to scale till 100,000 CPU cores. The new hybrid algorithm is developed on a base of a multi-level parallel model that exploits several layers of parallelism of a modern hybrid supercomputer. In this model, MPI and OpenMP are used on the first two levels to couple nodes of a supercomputer and to engage its CPU cores. Then, computing accelerators are further used by means of the hardware independent OpenCL computing standard. In this way, the implementation is adapted to a general computing model with central processors and math co-processors. In this paper the work is focused on adapting the basic operations of the algorithm to architectures of Graphics Processing Units (GPU) without considering the multi-GPU communication scheme. Technology of porting the code to OpenCL is described, certain optimization approaches are presented and relevant performance results obtaining up to 80–90 GFLOPS on a GPU accelerator are demonstrated. Moreover, the experience with different GPU architectures is summarized and a comparison based on the particular application is given for AMD and NVIDIA GPUs as well as for CUDA and OpenCL frameworks. [Copyright &y& Elsevier]

Highlights: [•] A multi-color Gauss–Seidel (MCGS) method is developed for OpenMP parallelization. [•] Using MCGS, an unstructured-grid Navier–Stokes solver is parallelized. [•] An algorithm for painting multi-colors is developed for MCGS. [•] The uniqueness of the solution of MCGS is confirmed with double precision accuracy. [Copyright &y& Elsevier]

Abstract: The NAME Atmospheric Dispersion Model is a Lagrangian particle model used by the Met Office to predict the propagation and spread of pollutants in the atmosphere. The model is routinely used in emergency response applications, where it is important to obtain results as quickly as possible. This requirement for a short runtime and the increase in core number of commonly available CPUs, such as the Intel Xeon series, has motivated the parallelisation of NAME in the OpenMP shared memory framework. In this work we describe the implementation of this parallelisation strategy in NAME and discuss the performance of the model for different setups. Due to the independence of the model particles, the parallelisation of the main compute intensive loops is relatively straightforward. The random number generator for modelling sub-grid scale turbulent motion needs to be adapted to ensure that different particles use independent sets of random numbers. We find that on Intel Xeon X5680 CPUs the model shows very good strong scaling up to 12 cores in a realistic emergency response application for predicting the dispersion of volcanic ash in the North Atlantic airspace. We implemented a mechanism for asynchronous reading of meteorological data from disk and demonstrate how this can reduce the runtime if disk access plays a significant role in a model run. To explore the performance on different chip architectures we also ported the part of the code which is used for calculating the gamma dose from a cloud of radioactive particles to a graphics processing unit (GPU) using CUDA-C. We were able to demonstrate a significant speedup of around one order of magnitude relative to the serial CPU version. [Copyright &y& Elsevier]

Abstract: The consequences that the recent stagnation in clockrates for CPUs has had on large-scale CFD runs are examined. At first sight, the conclusion is that only massive parallelism at the loop or domain decomposition level will lead to higher FLOP counts. However, the significant differences in advances for CPUs/GPUs versus RAM and interprocessor communication bandwidth lead to a so-called ‘limiting domain size’, below which communication dominates execution times and performance degrades drastically. The consequences of this ‘red-shift’ for the future of CFD are manifold: the time to advance the solution one timestep is limited, implying that even with unlimited number of processors/cores, LES, DES and DNS runs for realistic Reynolds-numbers will require days or weeks of execution. [Copyright &y& Elsevier]

Highlights: [•] Unstructured Navier–Stokes solver written in OpenCL. [•] Addressing code and performance portability. [•] A new multi-line vectorization strategy for the SpMV algorithm. [•] Adaptively choosing the optimal settings by running an optimization step. [•] Good speedup for real examples against a fully optimized parallel OpenMP solver. [Copyright &y& Elsevier]

Abstract: In this paper, development of a new Transmission Mesh method and its parallel implementation are presented. The Transmission Mesh method is a method where the body can move freely in the main mesh that covers the whole of the flow field. The method is developed for simulation of flows around multiple bodies that move and interact with each other in accordance with six degrees of freedom motion. An incompressible fluid/rigid sphere interaction is simulated. The parallel implementation based on the Domino method, which is a kind of block hyper-plane parallelism, is applied to the time-consuming Poisson solver. The parallel computation is performed using a multi-core processor with up to 24 cores, and its performance is tested. [Copyright &y& Elsevier]

Abstract: A new high-resolution code for the direct numerical simulation of a zero pressure gradient turbulent boundary layers over a flat plate has been developed. Its purpose is to simulate a wide range of Reynolds numbers from Re θ =300 to 6800 while showing a linear weak scaling up to 32,768 cores in the BG/P architecture. Special attention has been paid to the generation of proper inflow boundary conditions. The results are in good agreement with existing numerical and experimental data sets. [Copyright &y& Elsevier]

Abstract: Much of the current focus in high performance computing (HPC) for computational fluid dynamics (CFD) deals with grid based methods. However, parallel implementations for new meshfree particle methods such as Smoothed Particle Hydrodynamics (SPH) are less studied. In this work, we present optimizations for both central processing units (CPU) and graphics processing units (GPU) focused on a Lagrangian Smoothed Particle Hydrodynamics (SPH) method. In particular, the obtained performance and a comparison between the most efficient implementations for CPU and GPU are shown using the DualSPHysics code. [Copyright &y& Elsevier]

Abstract: The finite-difference time-domain method (FDTD) allows electromagnetic field distribution analysis as a function of time and space. The method is applied to analyze holographic volume gratings (HVGs) for the near-field distribution at optical wavelengths. Usually, this application requires the simulation of wide areas, which implies more memory and time processing. In this work, we propose a specific implementation of the FDTD method including several add-ons for a precise simulation of optical diffractive elements. Values in the near-field region are computed considering the illumination of the grating by means of a plane wave for different angles of incidence and including absorbing boundaries as well. We compare the results obtained by FDTD with those obtained using a matrix method (MM) applied to diffraction gratings. In addition, we have developed two optimized versions of the algorithm, for both CPU and GPU, in order to analyze the improvement of using the new NVIDIA Fermi GPU architecture versus highly tuned multi-core CPU as a function of the size simulation. In particular, the optimized CPU implementation takes advantage of the arithmetic and data transfer streaming SIMD (single instruction multiple data) extensions (SSE) included explicitly in the code and also of multi-threading by means of OpenMP directives. A good agreement between the results obtained using both FDTD and MM methods is obtained, thus validating our methodology. Moreover, the performance of the GPU is compared to the SSE+OpenMP CPU implementation, and it is quantitatively determined that a highly optimized CPU program can be competitive for a wider range of simulation sizes, whereas GPU computing becomes more powerful for large-scale simulations. [Copyright &y& Elsevier]

Abstract: We present C programming language versions of earlier published Fortran programs (Muruganandam and Adhikari (2009) ) for calculating both stationary and non-stationary solutions of the time-dependent Gross–Pitaevskii (GP) equation. The GP equation describes the properties of dilute Bose–Einstein condensates at ultra-cold temperatures. C versions of programs use the same algorithms as the Fortran ones, involving real- and imaginary-time propagation based on a split-step Crank–Nicolson method. In a one-space-variable form of the GP equation, we consider the one-dimensional, two-dimensional, circularly-symmetric, and the three-dimensional spherically-symmetric harmonic-oscillator traps. In the two-space-variable form, we consider the GP equation in two-dimensional anisotropic and three-dimensional axially-symmetric traps. The fully-anisotropic three-dimensional GP equation is also considered. In addition to these twelve programs, for six algorithms that involve two and three space variables, we have also developed threaded (OpenMP parallelized) programs, which allow numerical simulations to use all available CPU cores on a computer. All 18 programs are optimized and accompanied by makefiles for several popular C compilers. We present typical results for scalability of threaded codes and demonstrate almost linear speedup obtained with the new programs, allowing a decrease in execution times by an order of magnitude on modern multi-core computers. New version program summary Program title: GP-SCL package, consisting of: (i) imagtime1d, (ii) imagtime2d, (iii) imagtime2d-th, (iv) imagtimecir, (v) imagtime3d, (vi) imagtime3d-th, (vii) imagtimeaxial, (viii) imagtimeaxial-th, (ix) imagtimesph, (x) realtime1d, (xi) realtime2d, (xii) realtime2d-th, (xiii) realtimecir, (xiv) realtime3d, (xv) realtime3d-th, (xvi) realtimeaxial, (xvii) realtimeaxial-th, (xviii) realtimesph. Catalogue identifier: AEDU_v2_0. Program Summary URL: http://cpc.cs.qub.ac.uk/summaries/AEDU_v2_0.html. Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland. Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html. No. of lines in distributed program, including test data, etc.: 180583. No. of bytes in distributed program, including test data, etc.: 1188688. Distribution format: tar.gz. Programming language: C and C/OpenMP. Computer: Any modern computer with C language compiler installed. Operating system: Linux, Unix, Mac OS, Windows. RAM: Memory used with the supplied input files: 2–4 MB (i, iv, ix, x, xiii, xvi, xvii, xviii), 8 MB (xi, xii), 32 MB (vii, viii), 80 MB (ii, iii), 700 MB (xiv, xv), 1.2 GB (v, vi). Number of processors used: For threaded (OpenMP parallelized) programs, all available CPU cores on the computer. Classification: 2.9, 4.3, 4.12. Catalogue identifier of previous version: AEDU_v1_0. Journal reference of previous version: Comput. Phys. Commun. 180 (2009) 1888. Does the new version supersede the previous version?: No. Nature of problem: These programs are designed to solve the time-dependent Gross–Pitaevskii (GP) nonlinear partial differential equation in one-, two- or three-space dimensions with a harmonic, circularly-symmetric, spherically-symmetric, axially-symmetric or fully anisotropic trap. The GP equation describes the properties of a dilute trapped Bose–Einstein condensate. Solution method: The time-dependent GP equation is solved by the split-step Crank–Nicolson method by discretizing in space and time. The discretized equation is then solved by propagation, in either imaginary or real time, over small time steps. The method yields solutions of stationary and/or non-stationary problems. Reasons for the new version: Previous Fortran programs  are used within the ultra-cold atoms  and nonlinear optics  communities, as well as in various other fields . This new version represents translation of all programs to the C programming language, which will make it accessible to the wider parts of the corresponding communities. It is well known that numerical simulations of the GP equation in highly experimentally relevant geometries with two or three space variables are computationally very demanding, which presents an obstacle in detailed numerical studies of such systems. For this reason, we have developed threaded (OpenMP parallelized) versions of programs imagtime2d, imagtime3d, imagtimeaxial, realtime2d, realtime3d, realtimeaxial, which are named imagtime2d-th, imagtime3d-th, imagtimeaxial-th, realtime2d-th, realtime3d-th, realtimeaxial-th, respectively. shows the scalability results obtained for OpenMP versions of programs realtime2d and realtime3d. As we can see, the speedup is almost linear, and on a computer with the total of 8 CPU cores we observe in (a) a maximal speedup of around 7, or roughly 90% of the ideal speedup, while on a computer with 12 CPU cores we find in (b) that the maximal speedup is around 9.6, or 80% of the ideal speedup. Such a speedup represents significant improvement in the performance. Summary of revisions: All Fortran programs from the previous version  are translated to C and named in the same way. The structure of all programs is identical. We have introduced the use of comprehensive input files, where all parameters are explained in detail and can be set by a user. We have also included makefiles with tested and verified settings for GNU’s gcc compiler, Intel’s icc compiler, IBM’s xlc compiler, PGI’s pgcc compiler, and Oracle’s suncc (former Sun’s) compiler. In addition to this, 6 new threaded (OpenMP parallelized) programs are supplied (imagtime2d-th, imagtime3d-th, imagtimeaxial-th, realtime2d-th, realtime3d-th, realtimeaxial-th) for algorithms involving two or three space variables. They are written by OpenMP-parallelizing the most computationally demanding loops in functions performing time evolution (calcnu, calclux, calcluy, calcluz), normalization (calcnorm), and calculation of physical quantities (calcmuen, calcrms). Since some of the dynamically allocated array variables are used within such loops, they had to be made private for each thread. This was done by allocating matrices instead of arrays, with the first index in all such matrices corresponding to a thread number. Additional comments: This package consists of 18 programs, see Program title above, out of which 12 programs (i, ii, iv, v, vii, ix, x, xi, xiii, xiv, xvi, xviii) are serial, while 6 programs (iii, vi, viii, xii, xv, xvii) are threaded (OpenMP parallelized). For the particular purpose of each program, please see descriptions below. Running time: All running times given in descriptions below refer to programs compiled with gcc on quad-core Intel Xeon X5460 at 3.16 GHz (CPU1), and programs compiled with icc on quad-core Intel Nehalem E5540 at 2.53 GHz (CPU2). With the supplied input files, running times on CPU1 are: 5 min (i, iv, ix, xii, xiii, xvii, xviii), 10 min (viii, xvi), 15 min (iii, x, xi), 30 min (ii, vi, vii), 2 h (v), 4 h (xv), 15 h (xiv). On CPU2, running times are: 5 min (i, iii, iv, viii, ix, xii, xiii, xvi, xvii, xviii), 10 min (vi, x, xi), 20 min (ii, vii), 1 h (v), 2 h (xv), 12 h (xiv). [Copyright &y& Elsevier]

Abstract: Parallel computation was implemented for multiphase flow simulation of droplet collision using a front-tracking method. Various parallel interfaces were tested, i.e., MPI, OpenMP, and a hybrid of them. It was found that when the array size was comparatively small, OpenMP provided a more efficient platform for the simulation because of faster communication through the hardware memories. If the array size became larger than a critical magnitude, however, MPI could facilitate the computation due to outweighing operation of CPU instead of inter-networking communication. As a consequence, the coupled algorithm using both platforms would be effective as the array size was moderately large. [Copyright &y& Elsevier]

Abstract: The CPU/GPU heterogeneous clusters are important platforms for high performance computing applications. However, there are many challenges for efficiently performing the scientific and engineering legacy code on these heterogeneous systems. In this paper, we endeavor to address the programming-model issue by combining the existing models (i.e., MPI, OpenMP and CUDA). First, two hybrid programming patterns are presented, namely the and . Second, three kernels (i.e., EP, CG and MG) of the NAS parallel benchmarks (NPBs), which are abstracted from many legacy computational fluid dynamics applications, are implemented with the above two patterns. Third, these hybrid implementations are executed on the TianHe-1A supercomputer, and the corresponding experimental results show that significant performance improvement can be achieved with the above patterns. Finally, a detailed performance analysis about the two hybrid patterns is performed and some guidelines for porting the legacy code onto large-scale heterogeneous CPU/GPU clusters are also given. [Copyright &y& Elsevier]

Abstract: In this paper, a hybrid programming (OpenMP, MPI and CUDA) approach is used to study the performance of a parallelized Dynamic Discrete Ordinate Method (DDOM) solver . The parallel computation performances were compared under different scenarios. A hybrid parallelism of MPI and OpenMP performs well in terms of parallel efficiency (>90%) on a 64 core CPU cluster without using any load-balancing technique. This hybrid parallelism model is extended to a GPU cluster. By using massive multicore GPUs, the CUDA-accelerated code achieves a speed 250 times faster with a single GPU and over 780 times faster with a Quad-GPU cluster versus the identical process running on a single thread of CPU. Our results demonstrate that DDOM solver provides good scalability on CPU and GPU clusters. [Copyright &y& Elsevier]

Abstract: We describe a parallelised version of the MOLDY molecular dynamics program. This Fortran code is aimed at systems which may be described by short-range potentials and specifically those which may be addressed with the embedded atom method. This includes a wide range of transition metals and alloys. MOLDY provides a range of options in terms of the molecular dynamics ensemble used and the boundary conditions which may be applied. A number of standard potentials are provided, and the modular structure of the code allows new potentials to be added easily. The code is parallelised using OpenMP and can therefore be run on shared memory systems, including modern multicore processors. Particular attention is paid to the updates required in the main force loop, where synchronisation is often required in OpenMP implementations of molecular dynamics. We examine the performance of the parallel code in detail and give some examples of applications to realistic problems, including the dynamic compression of copper and carbon migration in an iron–carbon alloy. Program summary: Program title: MOLDY Catalogue identifier: AEJU_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEJU_v1_0.html Program obtainable from: CPC Program Library, Queenʼs University, Belfast, N. Ireland Licensing provisions: GNU General Public License version 2 No. of lines in distributed program, including test data, etc.: 382 881 No. of bytes in distributed program, including test data, etc.: 6 705 242 Distribution format: tar.gz Programming language: Fortran 95/OpenMP Computer: Any Operating system: Any Has the code been vectorised or parallelized?: Yes. OpenMP is required for parallel execution RAM: 100 MB or more Classification: 7.7 Nature of problem: Moldy addresses the problem of many atoms (of order 106) interacting via a classical interatomic potential on a timescale of microseconds. It is designed for problems where statistics must be gathered over a number of equivalent runs, such as measuring thermodynamic properities, diffusion, radiation damage, fracture, twinning deformation, nucleation and growth of phase transitions, sputtering etc. In the vast majority of materials, the interactions are non-pairwise, and the code must be able to deal with many-body forces. Solution method: Molecular dynamics involves integrating Newtonʼs equations of motion. MOLDY uses verlet (for good energy conservation) or predictor–corrector (for accurate trajectories) algorithms. It is parallelised using open MP. It also includes a static minimisation routine to find the lowest energy structure. Boundary conditions for surfaces, clusters, grain boundaries, thermostat (Nose), barostat (Parrinello–Rahman), and externally applied strain are provided. The initial configuration can be either a repeated unit cell or have all atoms given explictly. Initial velocities are generated internally, but it is also possible to specify the velocity of a particular atom. A wide range of interatomic force models are implemented, including embedded atom, Morse or Lennard-Jones. Thus the program is especially well suited to calculations of metals. Restrictions: The code is designed for short-ranged potentials, and there is no Ewald sum. Thus for long range interactions where all particles interact with all others, the order-N scaling will fail. Different interatomic potential forms require recompilation of the code. Additional comments: There is a set of associated open-source analysis software for postprocessing and visualisation. This includes local crystal structure recognition and identification of topological defects. Running time: A set of test modules for running time are provided. The code scales as order N. The parallelisation shows near-linear scaling with number of processors in a shared memory environment. A typical run of a few tens of nanometers for a few nanoseconds will run on a timescale of days on a multiprocessor desktop. [Copyright &y& Elsevier]

Abstract: Parallel implementation of a three-dimensional direct simulation Monte Carlo (DSMC) code employing complex data structures and dynamic memory allocation is detailed for shared memory systems using Open Multi-Processing (OpenMP). Several techniques to optimize the serial implementation of the DSMC method are first discussed. Specifically for a 3-level Cartesian grid, a Cartesian-based movement technique including particle indexing is demonstrated to result in a modest decrease in overall simulation expense of 34% compared with a ray-tracing technique combined with stored cell-connectivity. Two strategies for data localization leading to optimal usage of cache memory are demonstrated to speed up certain cell-based functions (such as collision computations) by a factor of 3.38–4.36. The shared-memory parallel implementation using OpenMP is then described in detail. Synchronization points and related critical sections are identified as major factors that impact the OpenMP parallel performance. Techniques to remove all such synchronization points in the OpenMP implementation of the DSMC method are outlined. For dual-core and quad-core systems, speedups of 1.99 and 3.74, respectively, are obtained for a (free-stream flow) test simulation with low granularity. Finally, the parallel performance of identical source code employing OpenMP is shown to be strongly correlated to the underlying computer architecture. Both Symmetric Multiprocessor (SMP) and non-uniform memory access (NUMA) systems are studied in order to achieve a better understanding of their impacts on parallel scalability when using OpenMP. [ABSTRACT FROM AUTHOR]

Abstract: Nowadays, NVIDIA''s CUDA is a general purpose scalable parallel programming model for writing highly parallel applications. It provides several key abstractions – a hierarchy of thread blocks, shared memory, and barrier synchronization. This model has proven quite successful at programming multithreaded many core GPUs and scales transparently to hundreds of cores: scientists throughout industry and academia are already using CUDA to achieve dramatic speedups on production and research codes. In this paper, we propose a parallel programming approach using hybrid CUDA OpenMP, and MPI programming, which partition loop iterations according to the number of C1060 GPU nodes in a GPU cluster which consists of one C1060 and one S1070. Loop iterations assigned to one MPI process are processed in parallel by CUDA run by the processor cores in the same computational node. [Copyright &y& Elsevier]

Abstract: The majority of nuclei available for study in solid state Nuclear Magnetic Resonance have half-integer spin , with corresponding electric quadrupole moment. As such, they may couple with a surrounding electric field gradient. This effect introduces anisotropic line broadening to spectra, arising from distinct chemical species within polycrystalline solids. In Multiple Quantum Magic Angle Spinning (MQMAS) experiments, a second frequency dimension is created, devoid of quadrupolar anisotropy. As a result, the center of gravity of peaks in the high resolution dimension is a function of isotropic second order quadrupole and chemical shift alone. However, for complex materials, these parameters take on a stochastic nature due in turn to structural and chemical disorder. Lineshapes may still overlap in the isotropic dimension, complicating the task of assignment and interpretation. A distributed computational approach is presented here which permits simulation of the two-dimensional MQMAS spectrum, generated by random variates from model distributions of isotropic chemical and quadrupole shifts. Owing to the non-convex nature of the residual sum of squares (RSS) function between experimental and simulated spectra, simulated annealing is used to optimize the simulation parameters. In this manner, local chemical environments for disordered materials may be characterized, and via a re-sampling approach, error estimates for parameters produced. Program summary: Program title: mqmasOPT Catalogue identifier: AEEC_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEEC_v1_0.html Program obtainable from: CPC Program Library, Queen''s University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 3650 No. of bytes in distributed program, including test data, etc.: 73 853 Distribution format: tar.gz Programming language: C, OCTAVE Computer: UNIX/Linux Operating system: UNIX/Linux Has the code been vectorised or parallelized?: Yes RAM: Example: (1597 powder angles)×(200 Samples)×(81 F2 frequency pts)×(31 F1 frequency points)=3.5M, SMP AMD opteron Classification: 2.3 External routines: OCTAVE (http://www.gnu.org/software/octave/), GNU Scientific Library (http://www.gnu.org/software/gsl/), OPENMP (http://openmp.org/wp/) Nature of problem: The optimal simulation and modeling of multiple quantum magic angle spinning NMR spectra, for general systems, especially those with mild to significant disorder. The approach outlined and implemented in C and OCTAVE also produces model parameter error estimates. Solution method: A model for each distinct chemical site is first proposed, for the individual contribution of crystallite orientations to the spectrum. This model is averaged over all powder angles [1], as well as the (stochastic) parameters; isotropic chemical shift and quadrupole coupling constant. The latter is accomplished via sampling from a bi-variate Gaussian distribution, using the Box-Muller algorithm to transform Sobol (quasi) random numbers [2]. A simulated annealing optimization is performed, and finally the non-linear jackknife [3] is applied in developing model parameter error estimates. Additional comments: The distribution contains a script, mqmasOpt.m, which runs in the OCTAVE language workspace. Running time: Example: (1597 powder angles)×(200 Samples)×(81 F2 frequency pts)×(31 F1 frequency points)=58.35 seconds, SMP AMD opteron. References: [[1]] S.K. Zaremba, Annali di Matematica Pura ed Applicata 73 (1966) 293. [[2]] H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods, SIAM, 1992. [[3]] T. Fox, D. Hinkley, K. Larntz, Technometrics 22 (1980) 29. [Copyright &y& Elsevier]

Abstract: We present a software library for numerically estimating first and second order partial derivatives of a function by finite differencing. Various truncation schemes are offered resulting in corresponding formulas that are accurate to order , , and , h being the differencing step. The derivatives are calculated via forward, backward and central differences. Care has been taken that only feasible points are used in the case where bound constraints are imposed on the variables. The Hessian may be approximated either from function or from gradient values. There are three versions of the software: a sequential version, an OpenMP version for shared memory architectures and an MPI version for distributed systems (clusters). The parallel versions exploit the multiprocessing capability offered by computer clusters, as well as modern multi-core systems and due to the independent character of the derivative computation, the speedup scales almost linearly with the number of available processors/cores. Program summary: Program title: NDL (Numerical Differentiation Library) Catalogue identifier: AEDG_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEDG_v1_0.html Program obtainable from: CPC Program Library, Queen''s University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 73 030 No. of bytes in distributed program, including test data, etc.: 630 876 Distribution format: tar.gz Programming language: ANSI FORTRAN-77, ANSI C, MPI, OPENMP Computer: Distributed systems (clusters), shared memory systems Operating system: Linux, Solaris Has the code been vectorised or parallelized?: Yes RAM: The library uses internal storage, N being the dimension of the problem Classification: 4.9, 4.14, 6.5 Nature of problem: The numerical estimation of derivatives at several accuracy levels is a common requirement in many computational tasks, such as optimization, solution of nonlinear systems, etc. The parallel implementation that exploits systems with multiple CPUs is very important for large scale and computationally expensive problems. Solution method: Finite differencing is used with carefully chosen step that minimizes the sum of the truncation and round-off errors. The parallel versions employ both OpenMP and MPI libraries. Restrictions: The library uses only double precision arithmetic. Unusual features: The software takes into account bound constraints, in the sense that only feasible points are used to evaluate the derivatives, and given the level of the desired accuracy, the proper formula is automatically employed. Running time: Running time depends on the function''s complexity. The test run took 15 ms for the serial distribution, 0.6 s for the OpenMP and 4.2 s for the MPI parallel distribution on 2 processors. [Copyright &y& Elsevier]

Recently a splitting approach has been presented for the simulation of sonic-boom propagation. Splitting methods allow one to divide complicated partial differential equations into simpler parts that are solved by specifically tailored numerical schemes. The present work proposes a second order exponential integrator for the numerical solution of sonic-boom propagation modelled through a dispersive equation with Burgers' nonlinearity. The linear terms are efficiently solved in frequency space through FFT, while the nonlinear terms are efficiently solved by a WENO scheme. The numerical method is designed to be highly parallelisable and therefore takes full advantage of modern computer hardware. The new approach also improves the accuracy compared to the splitting method and it reduces oscillations. The enclosed numerical results illustrate that parallelisation on a CPU results in a speedup of 22 times faster than the straightforward sequential version. The GPU implementation further accelerates the runtime by a factor 3, which improves to 5 when single precision is used instead of double precision. [ABSTRACT FROM AUTHOR]