Mikolaj Szydlarski
email: mikolaj.szydlarski@inria.fr

Pierre Esterie
email: pierre.esterie@lri.fr

Joel Falcou
email: joel.falcou@lri.fr

Laura Grigori
email: laura.grigori@inria.fr

Radek Stompor

Parallel Spherical Harmonic Transforms on heterogeneous architectures (GPUs/multi-core CPUs)

Keywords: Spherical Harmonic Transforms, hybrid architectures, hybrid programming, OpenMP, CUDA, Multi-GPU

Spherical Harmonic Transforms (SHT) are at the heart of many scientific and practical applications ranging from climate modelling to cosmological observations. In many of these areas new, cutting-edge science goals have been recently proposed requiring simulations and analyses of experimental or observational data at very high resolutions and of unprecedented volumes. Both these aspects pose formidable challenge for the currently existing implementations of the transforms.

This paper describes parallel algorithms for computing the SHTs with two variants of intra-node parallelism appropriate for novel supercomputer architectures, multi-core processors and Graphic Processing Units (GPU) and discusses their performance tests, alone and embedded within a top-level, MPI-based parallelization layer ported from the S 2 HAT library, in terms of their accuracy, overall efficiency and scalability. We show that our inverse SHTs with GeForce 400 Series GPUs equipped with latest CUDA architecture ("Fermi") outperforms the state of the art implementation for a multi-core processor executed on a current Intel Core i7-2600K. Furthermore, we show that an MPI/CUDA version of the inverse transform run on a cluster of 128 NVIDIA Tesla S1070 is as much as 3 times faster than the hybrid MPI/OpenMP version executed on the same number of quad-core processors Intel Nahalem for problem sizes motivated by our target applications. For the direct transforms, the performance is however found to be at the best comparable. Here we discuss in detail optimizations of two major steps involved in the transforms calculation, demonstrating how the overall performance efficiency can be obtained, and elucidating the sources of the dichotomy between the direct and the inverse operations.

Parallel Spherical Harmonic Transforms on heterogeneous architectures (GPUs/multi-core CPUs)

Résumé : Pas de résumé Mots-clés : Pas de motclef

Introduction

Spherical harmonic functions define an orthonormal complete basis for signals defined on a 2-dimensional sphere. Spherical harmonic transforms are therefore common in all scientific applications where such signals are encountered. These include a number of diverse areas ranging from weather forecasts and climate modelling, through geophysics and planetology, to various applications in astrophysics and cosmology. In these contexts a direct Spherical Harmonic Transform (SHT) is used to calculate harmonic domain representations of the signals, which often possesses simpler properties and are therefore more amenable to a further investigation. An inverse SHT is then used to synthesize a sky image given its harmonic representation. Both of those are also used as means of facilitating the multiplication of huge matrices, defined on the sphere, and having a property of being diagonal in the harmonic domain. Such matrices play an important role in statistical considerations of signals, which are statistically isotropic and are among key operations involved in Monte Carlo Markov Chain sampling approaches used to study such signals, e.g. [START_REF] Tanner | Tools for statistical inference. Observed data and data augmentation methods[END_REF].

Figure 1: Sky map example synthesised using our MPI/CUDA implementation of the alm2map routine. The units are µK. An overall monopole term, with an amplitude of ∼ 2.7 10 6 µK, has been subtracted to uncover the minute fluctuations.

The specific goal of this work is to assist simulation and analysis efforts related to the Cosmic Microwave Background (CMB) studies. CMB is an electromagnetic radiation left over after the hot and dense initial phase of the evolution of the Universe, which is popularly referred to as the Big Bang. Its observations are one of the primary ways to study the Universe. The CMB photons reaching us from afar carry an image of the Universe at the time when it had just a few percent of its current age of ∼ 13Gyears. This image confirms that the early Universe was nearly homogeneous and only small, 1 part in 10 5 , deviations were present, as displayed in Fig 1 . Those initiated the process of so-called structure formation, which eventually led to the Universe as we observe today. The studies of the statistical properties of these small deviations are one of the major pillars on which the present day cosmology is built. On the observational front the CMB anisotropies are targeted by an entire slew of operating and forthcoming experiments. These include a currently observing European satellite called Planck 1 , which follows in the footsteps of two very successful American missions, COBE 2 and WMAP 3 , as well as a number of balloon-borne and ground-based observatories. Owing to quickly evolving detector technology, the volumes of the data collected by these experiments have been increasing at the Moore's law rate, reaching at the present the values on order of petabytes. These new data aiming at progressively more challenging science goals not only provide, or are expected to do so, images of the sky with an unprecedented resolution, but their scientific exploitation will require intensive, high precision, and voluminous simulations, which in turn will require highly efficient novel approaches to a calculation of SHTs. For instance, the current and forthcoming balloon-borne and ground-based experiments will produce maps of the sky containing as many as O(10 5 -10 6 ) pixels and up to O(10 6 -10 7 ) harmonic modes. Maps from already operating Planck satellite will consist of between O(10 6 ) and O (10 8 ) pixels and harmonic modes. The production of high precision simulations reproducing reliably polarized properties of the CMB signal, in particular generated due to the so called gravitational lensing, requires an effective number of pixels and harmonic modes, as big as O (10 9 ). These latter effects constitute some of the most exciting new lines of research in the CMB area. The SHT implementation, which could successfully address all this needs would not only have to scale well in the range of interest, but also be sufficiently quick to be used as part of massive Monte Carlo simulations or extensive sampling solvers, which are both important steps of the scientific exploitation of the CMB data.

At this time there are a few software packages available, which permit calculations of the SHT. These include, healpix [3], glesp [2], ccsht [1], libpsht [4], and s 2 hat [5], which are commonly used in the CMB research, and some others such as, spharmonickit/s2kit [6] and spherpack [7]. They propose different levels of parallelism and optimization, while implementing, with an exception of two last ones, essentially the same algorithm. Among those the libpsht package [START_REF] Reinecke | Libpsht -algorithms for efficient spherical harmonic transforms[END_REF] offers typically the best performance, due to efficient code optimizations, and is considered as the state of the art implementation on serial and sharedmemory platforms, at least for some of the popular sky pixelizations. Conversely, the implementation offered by the s 2 hat library is best adapted to the distributed memory machines, and fully parallelized and scalable with respect to both memory usage and calculation time. Moreover, it provides a convenient and flexible framework straightforwardly extensible to allow for multiple parallelization levels based on hybrid programming models. Given the specific applications considered in this work, the distributed memory parallelism seems inevitable and this is therefore the s 2 hat package, which we select as a starting point for this work.

The basic algorithm we use hereafter, and which is implemented in both libpsht and s 2 hat libraries as described in the next section, scales as O(R N 2 max ), where R N is the number of rings in the analyzed map and max fixes the resolution and thus determines the number of modes in the harmonic domain, equal to 2 max . For full sky maps, we have usually max ∝ R N ∝ n 1/2 pix , where n pix is the total number of sky pixels, and therefore the typical complexity is O( 3 max ) ∝ O(n 3/2 pix ). We note that further improvements of this overall scaling are possible, for instance, if multiple identical transforms have to, or can, be done simultaneously, what could lower the numerical cost to essentially that of a single transform. In this paper we however focus on the core algorithm and leave an investigation of these potential extensions to future work.

We note that alternative algorithms have been also proposed, some of which display superior complexity. In particular, Driscoll and Healy [START_REF] Driscoll | Computing fourier transforms and convolutions on the 2-sphere[END_REF] proposed a divide-and-conquer algorithm with a theoretical scaling of O(n pix ln 2 n pix ). This approach is limited to special equidistant sphere grids and being inherently numerically unstable requires corrective measures to ensure high precision results. This in turn affects its overall performance. For instance, the software spharmonickit/s2kit, which is the most widely used implementation of this approach, has been found a factor 3 slower than the healpix transforms implementing the standard O(n 3/2 pix ) method at the intermediate resolution of max ∼ 1024 [START_REF] Wiaux | Fast Directional Correlation on the Sphere with Steerable Filters[END_REF]. Other algorithms also exist and typically involve a precomputation step of the same complexity, but often less favorable prefactors, as the standard approach, and an actual calculation of the transforms, which typically exploits either the Fast Multipole Methods, e.g., [START_REF] Suda | A fast spherical harmonics transform algorithm[END_REF][START_REF] Tygert | Fast algorithms for spherical harmonic expansions, ii[END_REF] or matrix compression techniques, e.g., [START_REF] Mohlenkamp | A fast transform for spherical harmonics[END_REF][START_REF] Tygert | Fast algorithms for spherical harmonic expansions, iii[END_REF], to bring down its overall scaling to O(n pix ln n pix ) [START_REF] Suda | A fast spherical harmonics transform algorithm[END_REF], O(n pix ln 2 n pix ) [START_REF] Tygert | Fast algorithms for spherical harmonic expansions, ii[END_REF][START_REF] Tygert | Fast algorithms for spherical harmonic expansions, iii[END_REF], or O(n 5/4 pix ln n pix ) [START_REF] Mohlenkamp | A fast transform for spherical harmonics[END_REF]. The methods of this class typically require significant memory resources needed to store the precomputation products and are advantageous only if a sufficient number of successive transforms of the same type has to be performed in order to compensate for the precomputation costs. The most recent, and arguably satisfactory, implementation of such ideas is a package called wavemoth4 [START_REF] Seljebotn | Wavemoth-Fast Spherical Harmonic Transforms by Butterfly Matrix Compression[END_REF], which achieves a speed-up of the inverse SHT with respect to the libpsht library by a factor ranging from 3 to 6 for max ∼ 4000. In such a case the required extra memory is on order of 40GBytes and it depends strongly, i.e., ∝ 3 max , on the resolution. We also note that in some applications the need to use SHTs can be sometimes by-passed by resorting to approximate but numerically efficient means such as the Fast Fourier Transforms as for instance in the context of convolutions on the sphere as discussed in [START_REF] Elsner | ARKCoS: artifact-suppressed accelerated radial kernel convolution on the sphere[END_REF].

In many practical applications, as the ones driving this research, SHTs are just one of many processing steps, which need to be performed to accomplish a final task. In such a context, these are the memory high-water mark and algorithm scalability, which frequently emerge as the two most relevant requirements, as the cost of the SHT transforms is often subdominant and their complexity more favorable than that of many other typical operations. From the memory point of view, the standard algorithm, having the smallest memory footprint, remains therefore an attractive option. The important question is then whether its efficient, scalable, parallel implementation is indeed possible. Such a question is particularly pertinent in the context of the heterogeneous architectures and hybrid programming models and this is the context investigated in this work.

The past work on the parallelization of SHTs includes a parallel version of the two step algorithm of Driscoll and Healy introduced by [START_REF] Inda | On the efficient parallel computation of legendre transforms[END_REF], the algorithm of [START_REF] Drake | Algorithm 888: Spherical harmonic transform algorithms[END_REF] based on rephrasing the transforms as matrix operations and therefore making them effectively vectorizable, and a shared memory implementation available in the libpsht package. More recently, algorithms have been developed for NVIDIA General-Purpose Graphics Processing Units (GPGPU) [START_REF] Soman | Accelerating spherical harmonic transforms on the nvidia gpu[END_REF][START_REF] Hupca | Spherical harmonic transform with gpus[END_REF]. This latter work provided a direct motivation for the investigation presented here, which describes what, to the best of our knowledge, is the first hybrid design of parallel SHT algorithms suitable for a cluster of GPU-based architectures and current high-performance multi-core processors, and involving hybrid OpenMP/MPI and MPI/CUDA programming. We find that once carefully optimized and implemented, the algorithm displays nearly perfect scaling in both cases, at least up to 128 MPI processes mapped on the same number of pairs of multi-core CPU-GPU. We find that inverse SHT with GeForce 400 Series equipped with the latest CUDA device ("Fermi") outperforms state of the art implementation for a multi-core processors executed on latest Intel Core i7-2600K, while the direct transforms in both these cases perform comparably.

This paper is organised as follows. In section 2 we introduce the algebraic background of the spherical harmonic transforms. In section 3 we describe a basic, sequential algorithm and list useful assumptions concerning a sphere pixelisation, which facilitate a number of acceleration techniques used in our approach. In following section, we introduce a detailed description of our parallel algorithms along with two variants suitable for clusters of GPUs and clusters of multi-cores. Section 5 presents results for both implementation and finally section 6 concludes the paper.

Algebraic background

Definitions and notations

Just as the Fourier basis facilitates a quick and efficient evaluation of convolutions in one or more dimensional spaces, the spherical harmonic basis does so but for functions defined on the surface of a sphere, i.e., depending on spherical coordinates: θ ∈ (0, π) (colatitude) and φ ∈ [0, 2π) (the angle counterclockwise about positive z-axis from the positive x-axis). We can express a band-limited approximation f to a realvalued spherical function f (θ, φ) f : S 2 → R in terms of the spherical harmonic basis functions as

f (θ, φ) = max =0 m=- a m Y m (θ, φ) , (1)
where a m are spherical harmonic coefficients describing f in the harmonic domain, Y m is a spherical harmonic basis function of degree and order m, and the bar over f indicates that it is merely a bandlimited approximation of f with a bandwidth max . The spherical harmonic coefficients are computed by taking the scalar inner product of f with the corresponding basis function, which can be expressed as the integral

a lm = π θ=0 2π φ=0 dθ dφ f (θ, φ)Y † lm (θ, φ) sin θ, (2)
where † denotes complex conjugation.

In actual applications the spherical harmonic transforms are rather used to project grid point data (a discretized scalar field s n ) on the sphere onto the spectral modes in an analysis step (as displayed in equation ( 3)) and an inverse transform reconstructing the grid point data from the spectral information in a synthesis step (as displayed in equation ( 4)). Depending on the choice of max and the grid geometry the transform s n → a lm may only be approximate, what is indicated by a tilde over a.

ã m = {θn,φn} sn (θn, φn) Y m (θn, φn) , (3)
sn (θn, φn) = max =0 m=- a m Y m (θn, φn) . (4)
In CMB applications, s n is a vector of pixelized data, e.g. the brightness of incoming CMB photons, assigned to n pix locations (θ n , φ n ) on the sky defined as centres of suitable chosen sky pixels. The parameter max defines a maximum order of the Legendre function and thus a band-limit of the field s n , and in practice it is set by the experimental resolution.

The basis functions Y m (θ n , φ n ) are defined in terms of normalised associated Legendre functions P m of degree l and order m,

Y m (θn, φn) ≡ P m (cos θn) e imφn , (5)
where we use the relation Y l,-m = (-1) m Y † lm to compute all spherical harmonics in terms of the associated Legendre functions P lm with m ≥ 0. Consequently, by separating variables following equation (5) we can reduce the computation of the spherical harmonic transform to a regular Fourier transform in the longitudinal coordinate φ n followed by a projection onto the associated Legendre functions,

ã m = {θn}     {φn} s n (θ n , φ n ) e imφn   P m (cos θ n )   . (6)
The associated Legendre functions satisfy the following recurrence (with respect to the multipole number for a fixed value of m) critical for the algorithms developed in this paper,

P +2,m (x) = β +2,m x P +1,m (x) + β +2,m β +1,m P m (x) , (7)
where

β m = 4 2 -1 2 -m 2 . ( 8
)
The recurrence is initialised by the starting values,

P mm (x) = 1 2 m m! (2m + 1)! 4π 1 -x 2 m ≡ µ m 1 -x 2 m , (9)
P m+1,m (x) = β +1,m x P mm (x) . ( 10
)
The recurrence is numerically stable but a special care has to be taken to avoid under-or overflow for large values of max [START_REF] Górski | HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere[END_REF] e.g., for increasing m the P mm values can become extremely small such that they can no longer be represented by the binary floating-point numbers in IEEE 7545 standard. However, since we have freedom to rescale all the values of P m , we can dynamically rescale current values, while performing the recurrence, to prevent a potential overflow on a subsequent step. More precisely on each step the newly computed value of the associated Legendre function is tested and rescaled if found to be too close to the over-or underflow limits. The rescaling coefficients (e.g., in form of their logarithms) are kept track of and used to rescale back all the computed values of P m at the end as required. This scheme is based on two facts. First, the values of the associated Legendre functions calculated via the recurrence change gradually and rather slowly on each step. Second, their actual values as needed by the transforms are representable within the range of the double precision values. The specific implementation of these ideas used in all our algorithms follows that of the s 2 hat software and the HEALPix package. It is based on a use of a precomputed vector of values, sampling the dynamic range of the representable double precision numbers and thus avoids any explicit computation of numericallyexpensive logarithms and exponentials. This scaling vector is used to compare the values of P m computed on each step of the recurrence, and then it is use to rescale them if needed. For sake of simplicity, in this paper we assume that this is an integral part of associated Legendre function evaluation (independent of architecture) and we omit explicit description of algorithm related to rescaling. For more details we refer to [START_REF] Hupca | Spherical harmonic transform with gpus[END_REF] and documentation of the HEALPix package.

Specialized formulation

For a general grid of O(n pix ) points on the sphere, a direct calculation has computational complexity O( 2 max n pix ) which is a serious limitation for problems in range of our interest. For instance, the current and forthcoming balloon-borne and ground-based experiments will produce maps of sky containing as many as n pix ∈ [10 5 , 10 6 ] pixels and up to O(10 6 -10 7 ) harmonic modes. Maps from already operating Planck satellite will consist of between O(10 6 ) and O(10 8 ) pixels and harmonic modes. However, the algorithm we use hereafter and which is described in the next section, scales as O(n pix max ). We obtain such complexity with help of additional geometrical constraints imposed on the permissible pixelizations. These are [START_REF] Górski | HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere[END_REF],

• the map pixels are arranged on R N ≈ √ n pix iso-latitudinal rings, where latitude of each ring is identified by a unique polar angle θ n (for sake of simplicity hereafter we will refer to this angle by the ring index i.e., r n → θ n , where n ∈ [0, R N ])

• within each ring r n , pixels must be equidistant (∆φ =const), though their number can vary from a ring to another ring.

The requirement of the iso-latitude distribution for all pixels helps to save on the number of required FLOPs as the associated Legendre function need now be calculated only once for each pixel ring. Consequently, nearly all spherical grids, which are currently used in CMB analysis (for instance HEALPix and GLESP), as well as other fields conform with such assumptions.

Taking into account these restrictions and definition of basis function Y m (5) we can rewrite (adapted from [START_REF] Górski | HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere[END_REF]) the synthesis step from equation (4) as

s (rn, φn) = max m=-max e imφ 0 ∆ A m (rn) , (11)
where ∆ A m (r n ) is a set of functions such as

∆ A m (rn) ≡                        max =0 a 0 P 0 (cos rn) , m = 0, max =m a m P m (cos rn) , m > 0, max =|m| a † |m| P |m| (cos rn) , m < 0, (12)
and respectively the analysis step(3) as

ã m = R N n=0 ∆ S m (rn) P m (cos rn) , (13)
where

∆ S m (rn) ≡ {φn} s (rn, φn) e -imφ 0 . ( 14
)
The φ 0 denotes the φ angle of the first pixel in the ring.

From above, we can see that the spectral harmonic transforms can be split into two successive computations: one calculating the associated Legendre functions, Eqs. ( 12) and [START_REF] Górski | HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere[END_REF], and the other performing the series of Fast Fourier Transforms (FFTs), Eqs. [START_REF] Elsner | ARKCoS: artifact-suppressed accelerated radial kernel convolution on the sphere[END_REF] and ( 14) [START_REF] Muciaccia | Fast Spherical Harmonic Analysis: A Quick Algorithm for Generating and/or Inverting Full-Sky, High-Resolution Cosmic Microwave Background Anisotropy Maps[END_REF]. This splitting immensely facilitates implementation, for instance, in order to evaluate Fourier transforms we can use third-party libraries devoted to this problem. However, we have to bear in mind that different libraries offer different robustness with strong dependence on length of FFT. This focus our special attention since the number of samples (pixels) per ring may vary and this strongly affects efficiency of the FFT algorithms. For instance, well-know fftpack library [START_REF] Swarztrauber | Vectorizing the FFTs[END_REF] used in healpix [START_REF] Górski | HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere[END_REF] and libpsht [START_REF] Reinecke | Libpsht -algorithms for efficient spherical harmonic transforms[END_REF] libraries performs well only for a length N F F T , whose prime decomposition only contains the factors 2, 3, and 5, while for prime number lengths its complexity

degrades from O (N F F T log N F F T ) to O(N 2 F F T ).
For this reason we employ by default in all versions of our algorithm the FFTW library [START_REF] Frigo | The design and implementation of FFTW3[END_REF] which uses O (N F F T log N F F T ) algorithms even for prime sizes.

Time consumption

Three parameters typically determine the sizes and numerical complexity characteristic of our problem. These are a total number of pixels, n pix , a number of rings of pixels in the analyzed map, R N , and a maximum order of the associated Legendre functions, max . In particular a number of modes in the harmonic domain is ∼ 2 max . In well defined cases these three parameters are not completely independent. For example, in the case of full sky maps, we have usually max ∝ R N ∝ n 1/2 pix . This imply that recurrence step requires O R N 2 max floating point operations (FLOPS), due to the fact that for each of R N rings we have to calculate a full set, up to max , of the associated Legendre functions, P m , at cost ∝ 2 max . The Fourier transform step, as mentioned above, has then subdominant complexity O (R N max log max ). These theoretical expectations agree with our experimental results shown in Fig. 2 depicting a typical breakdown of the average overall time between the main steps involved in the SHTs computation as obtained for a run performed on quad-core Intel i7-2600K processor. We can observe that indeed the evaluation of the Legendre transform dominates by far over FFT in terms of time of computation. This motivated us to study the possible improvements of this part of the calculation. We explore this issue in the context of the heterogeneous architectures and hybrid programming models. In particular we introduce a parallel algorithm that is suitable for clusters of accelerators, namely multi-core processors and GPUs which we employ for an efficient evaluation of spherical harmonic transforms.

Basic algorithm

We proceed by developing algorithms which agree with assumptions listed in the previous section. Then, in the next section we introduce their parallel versions and we discuss possible improvements specific to the targeted architectures.

Algorithm 1 presents the computation of discrete, direct SHT which is a realisation of equation ( 3) for a given "bandwidth" max and input data vector s obtained from the samples of the original function which we wish to transform. Respectively, algorithm 2 presents the computation of discrete, inverse transform Algorithm 1 Direct transform (Eq. 3)

Require: s ∈ Rn pix values step 1 -∆ S m (r) calculation for every ring r do for every m = 0, ..., mmax do • calculate ∆ S m (r) via FFT Eq. ( 14) end for(m) end for(r) step 2 -Pre-computation for every m = 0, ..., mmax do

• µm pre-computation Eq. (9) end for (m) step 3 -Core calculation for every ring r do for every m = 0, ..., mmax do for every = m, ..., max do • compute P m via Eq. (7)

• update a m , Eq. (13) end for ( ) end for (m) end for (r) return a m ∈ C ×m from equation (3) which takes as input a set of complex numbers a m , interpreted as Fourier coefficients in a spherical harmonic expansion and returns a set of sample values i.e., real data vector s.

Similarities

Both algorithms use a divide and conquer approach i.e., in this approach our problem of computing projection onto associated Legendre functions is decomposed into smaller subproblems of a similar form (rings of pixels). The subproblems are solved recursively by a further subdivision (computation of P m via recurrence with respect to the multipole number for a fixed value of m) and finally their solutions are combined to solve the original problem. Furthermore, both algorithms consist of three main steps:

• the pre-computation of the starting values of the recurrence [START_REF] Drake | Algorithm 888: Spherical harmonic transform algorithms[END_REF] in the loop over all m ∈ [0, m max ] (step 2 in Algorithm 1 and step 1 in Algorithm 2 ),

• the evaluation of the associated Legendre function for each ring of pixels via the recurrence from equation (7). This requires three nested loops, where the loop over is set innermost to facilitate the evaluation of the recurrence with respect to the multipole number for a fixed value of m (step 3 in Algorithm 1 and step 2 in Algorithm 2 ).

• step where via FFTs we either generate the cosine transform representation of the input vector s (11) (step 1 in Algorithm 1) or to complete computation in inverse algorithm an abelian FFT is performed to compute the sums (14) (step 3 in Algorithm 2).

Differences

Both algorithms are almost the same in terms of number and type of arithmetic operations required for the computation (see for instance Figure 2 or a detailed analysis in [START_REF] Melton | An analysis of the spectral transfrom operations in climat and weather models[END_REF]). The main difference consists in the dimension of the partial results which need to be updated within the nested loop where we evaluate the associated Legendre functions (step 3 in Algorithm 1 and step 2 in Algorithm 2). In case of the direct transform, for each ring of pixels on the sphere we need to update all non-zero values of the matrix a m ∈ C max ×mmax by a product between the precomputed inner sums defined as ∆ S m (14) and associated Legendre function for each fixed m. In the case of the inverse transform, where the input is a set of complex numbers a m , the nested loop (step 2 in Algorithm 2) is used to compute the inner sum ∆ A m [START_REF] Frigo | The design and implementation of FFTW3[END_REF], which can be obtained as a matrix-vector product

∆ A m =      ∆ A m (r 0 , m) ∆ A m (r 1 , m) . . . ∆ A m (R N , m)      =     
P mm (cos r 0 ) . . . P max m (cos r 0 ) P mm (cos r 1 ) . . . P m axm (cos r 1 ) . . . . . . . . .

P mm (cos r R N ) . . . P max m (cos r R N )           a mm a m+1m . . . a max m      . ( 15
)
As we will see further in this report, those differences may cause big contrast in terms of evaluation time and memory consumption. This effect can be additionally enhanced by the type of the architecture on which we execute our algorithm.

Parallel Spherical Harmonic Transforms

Top-level distributed-memory parallel framework

The top level parallelism employed in this work is adopted from the S 2 HAT library developed by the last author. Though not new, it has not been described in detail in the literature yet and we present below its first detailed description including its performance model and analysis. The framework is based on a distributed-memory programming model and involves two computational stages separated by a single instance of a global communication involving all the processes performing the calculations. These together with a specialized, but flexible data distribution, are the major hallmarks of the framework. The two stages of the computations can be further parallelized and/or optimized and we discuss two examples of potential intra-node parallelization approaches below. We note that the framework as such does not imply what algorithm should be used by the MPI processes, e.g., [START_REF] Seljebotn | Wavemoth-Fast Spherical Harmonic Transforms by Butterfly Matrix Compression[END_REF], though of course its overall performance and scalability will depend on a specific choice.

The major motivation behind such a distributed memory parallelization level comes from the volumes of the data, both maps and harmonic coefficients, to be analyzed in the applications targeted here and from the need to maintain flexibility of the routines designed to be used in massively parallel applications, where SHTs typically constitute only one step of the overall process. The latter requirement calls not only for a flexible interface and an ability to use efficiently as many MPI processes as available, but also for a low memory footprint per process.

Data distribution and algorithm

The simple structure of the top-level parallelism outlined above is facilitated by a data layout, which assumes that both the map s and its harmonic coefficients a m , are distributed. The map distribution is ring-based, i.e., we assign to each process a part of the full map consisting of complete rings of pixels. To do so we make use of the assumed equatorial symmetry of the pixelization schemes and first divide all rings of one of the hemispheres (including the equatorial ring if present) into disjoint subsets made of complete, consecutive rings, and then assign to each process part of the map corresponding to one of such ring subsets together with its symmetric, counterpart from the other hemisphere. If we denote the rings mapped on a processor i as R i , then ∪ nprocs i=0 R i = R, where R is the set of all rings and n procs is the number of processes used in the computation. We also have R i ∩ R j = ∅ if i = j. An example ring distribution is depicted in Figure 3(a), where different colours mark pixels on iso-latitudinal rings mapped on two different processes.

In order to distribute harmonic domain object we divide the 2-dimensional array of the coefficients a m with respect to a number, m, and assign to each process a subset of values of m, M i , and a corresponding subset of the hamonic coefficients, {a m : m ∈ M i and ∈ [m, max ]}. As with the ring subsets, the subsets of m values have to be disjoint, M i ∩ M j = ∅ if i = j and include all the values, ∪ nprocs i=0 M i = M. An example of the harmonic domain distribution is depicted in Figure 3(b), where like in the case of the ring distribution we marked with different colours a m coefficients mapped on two different processes.

This data layout allows to avoid any redundancy in terms of calculations and memory and any need for an intermediate inter-process communication, while at the same time is general enough that at least in principle can accommodate an entire slew of possibilities as driven by real applications. Any freedom in selecting a specific layout should be used to ensure a proper load balancing in terms of both memory and work, and it may, and in general will, depend on what algorithm is used by each of the MPI processes. Though perfect load balancing, memory-and work-wise, may not be always possible, in practice we have found that good compromises usually exist. In particular, the map distribution which assigns the same number of rings to each process is typically satisfactory, even if it may in principle lead to an unbalanced work and/or memory distribution. The latter is indeed the case for the healpix pixelization which does not preserve the same number of samples per ring, what introduces differences in the memory usage between the processes but also in the workload. Some fine tuning could be used here to improve on both these aspects, as it is indeed done in s 2 hat, but even without that the differences are not critical. For the distribution of the harmonic objects we propose M i ≡ {m : m = i + k n procs or m = m max -i -k n procs , where k ∈ [0, m max /2]}, so the process i stores the values of m equal to [i, i + n procs , ..., m max -i -n procs , m max -i]. This usually provides a good first try, which appears satisfactory in a number of cases as it will be explained in the next section. We note that this data layout imposes an upper limit on the number of processes which could be used, and which is given by min(m max /2, R N /2). This is typically sufficiently generous not to lead to any problems, as usually we have R N ∝ m max max and the limit on a number of processes is O ( max ). With this data distribution, Algorithms 1 and 2 require only straightforward and minor modifications to become functional parallel codes. As an example, Algorithm 3 details the parallel inverse transform, which is a parallel version of Algorithm 2. The operations listed there are to be performed by each of the • initialise the recurrence: P mm , P m+1,m using precomputed µ m Eqs. ( 9) & (10)

• precompute recurrence coefficients, β m Eq. ( 8) for every = m + 2, ..., max do • compute P m via the 2-point recurrence, given precomputed β m , Eq. (7)

• update ∆ A m (r), Eq. ( 12) end for ( ) end for (r) and (m)

global communication • ∆ A m (r) , m ∈ M i , all r ∈ R MPI_Alltoallv ⇐⇒ ∆ A m (r) , r ∈ R i , all m ∈ M step 3 -s n calculation for every ring r ∈ R i do
• via FFT calculate s (r, φ) for all samples φ of ring r, given pre-computed ∆ A m (r) for all m. end for(r) return s n for all (r ∈ R i ) n proc MPI processes involved in the computation. Like its serial counterpart, the parallel algorithm can be also split into steps. First, we precompute the starting values of the Legendre recurrence µ m (9), but only results for m ∈ M i are preserved in memory. In next step, we calculate the functions ∆ A m using Eq. ( 12) for every ring r of the sphere, but only for m values in M i . Once this is done, a global (with respect to the number of processes involved in the computation) data exchange is performed so that at the end each process has in its memory the values of ∆ A m (r) computed only for rings r ∈ R i but for all m values. Finally, via FFT we synthesize partial results as in equation [START_REF] Elsner | ARKCoS: artifact-suppressed accelerated radial kernel convolution on the sphere[END_REF]. From the point of view of the framework, steps 1 and 2 constitute the first stage of the computation, followed by the global communication, and the second computation involving here a calculation of the FFT. As mentioned earlier, though the communication is an inherent part of the framework, the two computation stages can be adapted as needed. Hereafter, we adhere to the standard SHT algorithm as used in Algorithm 3 .

Performance analysis

Operation count. We note that the scalability of all the major data products can be assured by adopting an appropriate distribution of the harmonic coefficients and the map. This stems from the fact that apart of the input and output objects, i.e., maps and harmonic coefficients themselves, the volumes of which are inversely proportional to a number of the MPI processes by the construction, the largest data objects derived as intermediate products, are arrays ∆ A or ∆ S . Their sizes are in turn given as O(R N m max /n proc ), and again decrease inversely proportionally with n procs . These objects are computed during the first computation stage, subsequently redistributed as the result of the global communication, and used as an input for the second stage of the computation. Their total volume, as integrated over the number of processes, therefore also determines the total volume of the communication, as discussed in the next Section. As these three objects are indeed the biggest memory consumers in the routines the overall memory usage scales as (n pix + 2 m max max + R N m max )/n procs . For typical pixelization schemes, full sky maps, analyzed at its sampling limit, we have n pix ∼ 2 max , max m max , and

R N ∼ n 1/2
pix , and therefore each of these three data objects contribute nearly evenly to the total memory budget of the routines. Consequently, the actual high memory water mark is typically ∼ 50% higher than the volume of the input and output data. We note that if all three parameters, max , m max , R N , are allowed to assume arbitrary values memory re-using and, in particular, in-place implementations are not straightforward, as it may not be easy or possible to fit any of these data objects in the space assigned for others. For this reason, this kind of options are not implemented in the software discussed here.

Given the assumed data layout, also the number of floating point operations (FLOPS) required for the evaluation of FFTs and associated Legendre functions scales inversely proportionally to the number of the processes. The only exception is related to the cost of the initialization of the recurrence (Eq. 9), which is performed redundantly by each MPI process to avoid an extra communication. Consequently, each process receives at least one m value on order of m max and has to perform at least O(m max) operations as part of the pre-computation, equation [START_REF] Drake | Algorithm 888: Spherical harmonic transform algorithms[END_REF]. A summary of the FLOPs required in each step of the parallel SHTs algorithms is presented in Table 1, where we have assumed here that the subsets R i and M i contain R N /n proc and m max /n proc elements, respectively.

FLOPs Communication cost. In order to determine the overall performance of the framework we have to also estimate the communication cost. The data exchange is performed with help of a single collective operation of the type all-to-all in which each process sends data to, and receives from, every other process. We use a simple model to estimate the cost of the collective communication algorithms in terms of latency and bandwidth use i.e., we assume that the time taken to send a message between any two nodes can be modelled as α + β n, where α is the latency (or startup time) per message, independent of message size, β is the transfer time per byte, and n is the number of bytes transferred. However in case of the collective communication routines, different MPI libraries use different algorithms, which moreover may depend on a size of the message. For specificity, we base our analysis on MPICH [START_REF]MPICH -A portable implementation of MPI[END_REF], a widely used MPI library. Similar considerations can be performed for other libraries.

Pre-computation Step {2,1} O (m max ) Recurrence Step {3,2} O (R N max m max n proc FFTs Step {1,3} O R N n proc m max log m max
The all-to-all communication in MPICH [START_REF] Thakur | Optimization of collective communication operations in mpich[END_REF] for short messages (≤ 256 kilobytes per message) uses the index algorithm by Bruck et al. [START_REF] Bruck | Efficient algorithms for all-to-all communications in multiport message-passing systems[END_REF]. It is a store-and-forward algorithm that takes log n steps at the expense of some extra data communication ( β n log n instead of β n). For long messages and even number of processes, MPICH uses a pairwise-exchange algorithm, which takes n proc -1 steps. On each step, each process calculates its target process and exchanges data directly with that process i.e., data is directly communicated from source to destination, with no intermediate steps.

Note that in case of even distribution of rings r and m values, the size of message to be exchanged between each pair of processes can be defined as

S msg := R N m max n proc n C , (16)
where n C is the size in bytes of a complex number representation (usually n C = 16).

Finally, merging all the assumptions listed above we can derive a following formula for the total time require for a communication in our parallel SHTs algorithm. In Figure 4 we illustrate theoretical runtimes depending on the number of the employed processes (leftmost panel) and on the size of the problem (middle panel). The transform parameters have been set assuming a standard full sky, fully resolved case, and the HEALPix pixelization, i.e., max

= m max = 2 N side , R N = 4 N side -1, n pix = 12 N 2 side .
Here, N side is the HEALPix resolution parameter (for more details see beginning of the section 5 and [START_REF] Górski | HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere[END_REF]). For the hardware parameters, in the estimation of the communication cost we have used following [START_REF] Rugina | Predicting the running times of parallel programs by simulation[END_REF][START_REF] Sottile | Performance analysis of parallel programs via message-passing graph traversal[END_REF] α = 10 -5 , β = 10 -9 , while for the cost of calculations we have assumed that each MPI process attains the effective performance of 10Giga operations per second, resulting in a time of γ = 10 -10 seconds per FLOP. The choice we have made is meant to be only indicative, yet at its face value the latter number is more typical of the MPI processes being nodes of multi-cores than just single cores, in particular given the usual efficiency of such calculations, which is usually on order of 10 -20%. If less efficient MPI processes are used then the solid curves in Fig. 4 corresponding to the computation times should be shifted upwards by a relevant factor.

The dominant calculation is the recurrence needed to compute the associated Legendre functions in agreement with our earlier numerical result, Fig. 2. The computation scales nearly perfectly with the number of processed used decaying as ∝ 1/n proc and it increases with growing size of the problem, ∝ N 2 side . At the same time the communication cost is seen to be independent on a number of processes. This is because for the considered here numbers of processes the size of the message exchanged between any two processes is always large and the communication time is thus described by the second term of the second equation of Eqs. 17. Given that the single message size, S msgs , decreases linearly with a number of processes, the total communication time does not depend on it. The immediate consequence of these two observations is that both these times will be found comparable if sufficiently large number of processes is used. A number of processes at which such a transition takes place will depend on the constants entering in the calculation of Eqs. 17 and the size of the problem. This dependence of the critical number of processes on the size of the problem is shown in the rightmost panel of Fig. 4 with thick solid contour labelled 1.0. Clearly, from the perspective of the transform runtimes there is no extra gain, which could be obtained from using more processes than given by their critical value. In practice, a number of processes may need to be fixed by memory rather than computational efficiency considerations. As the memory consumption depends on the problem size, i.e., n pix , in the same way as the communication volume, by increasing the number of processes in unison with the size of the problem we can ensure that the transform runs efficiently, i.e., the fraction of the total execution time spent on communicating will not change, and that it will always fit within the machines memory.

We note that some improvement in the overall cost of the communication could be expected by employing a non-blocking communication pattern applied to smaller data objects, i.e., single rows of the arrays, ∆ A and ∆ S corresponding to a single ring or m value, successively and immediately after they become available. Though such an approach could extend the parameter space in which the computation time is dominant and the overall scaling favorable, it will most likely have only a limited impact given that the same, large, total volume of the data, which needs to be communicated. The communication volume could be decreased, and thus again the parameter space with the computation dominance extended, if redundant data distributions and redundant calculations are allowed for in the spirit of the currently popular communication avoiding algorithms, e.g., [START_REF] Grigori | Communication avoiding gaussian elimination[END_REF]. Though these kinds of approaches are clearly of interest they are outside of the scope of this work and are left here for future research.

We also not that the conclusions may disfavour algorithms which require on the one hand big memory overhead and/or need to communicate more. Such algorithms will tend to require a larger number of processes to be used for a given size of the problem and thus will quickly get into bad scaling regime. Their application domain may be therefore rather limited. We emphasize that this general conclusion should be reevaluated case-by-case for any specific algorithms. Likewise, it is clear that the regime of the scalability of the standard SHT algorithm could be extended by decreasing the communication volume. Given that the communicated objects are dense, the latter would typically require some lossy compression techniques, we do not consider them here. Instead we focus our attention on accelerating intra-node calculations. This could improve the transforms runtimes for a small number of processes, however, if successful unavoidably will also result in lowering the critical values for the processes numbers. leading to losing the scalability earlier. Nevertheless, the overall performance of these improved transforms would never be inferior to that of the standard ones whatever a number of processes.

Exploiting intra-node parallelism

In this Section we consider two extensions of the basic parallel algorithm, Algorithm 3, each going beyond its MPI-only paradigm. Our focus hereafter is therefore on the computational tasks as performed by each of the MPI processes within the top-level framework outlined earlier. Guided by the theoretical models from the previous Section we expect that any gain in performance of a single MPI process will translate directly into analogous gain in the overall performance of the transforms at least as long as the communication time is subdominant, i.e., in a relatively broad regime of cases of practical importance. Hereafter we therefore develop implementations suitable for multi-core processors and GPUs and discuss the benefits and trade-offs they imply.

As highlighted by Fig. 4, of the three computational steps, this is the recurrence used to compute the associated Legendre functions which takes by far the dominant fraction of the overall runtime. Consequently, we will pay special attention to accelerating this part of the algorithm, while we will keep on using standard libraries to compute the FFTs. The latter will be shown to deliver a sufficient performance given our goals here. The associated Legendre function recurrence involves three nested loops and we will particularly consider the benefits and downside of their different ordering in order to attain the best possible performance.

Hereafter, for simplicity we will refer to the algorithms for direct and inverse spherical transforms by the core names of their corresponding routines i.e., map2alm and alm2map respectively.

Multithreaded version for multi-core processors

Since the largest and fastest computers in the world today employ both shared and distributed memory architectures, software which use only a distributed memory approaches such as MPI may not fully exploit the shared memory underlying architecture and thus may loose some efficiency. Shared memory approaches, based on explicit multithreading or language-generated multithreading (e.g. OpenMP) are more accurate on shared memory architectures. Thus, many hybrid approaches, typically mixing MPI and OpenMP, have been proposed to better exploit clusters of of multi-core processors. In this spirit, we introduce the second level of parallelism based on multithreaded approach in the computation of SHTs.

As discussed shortly at the beginning of this section, the order of the three nested loops can be changed. For instance, the iterations over m in the outermost loop allow precomputing β m Eq. ( 8), P mm Eq. ( 9) and P m+1,m Eq. ( 10) only once for each ring with extra storage of an order O ( max ). But in this case we also need to store intermediate products which require additional O (R N max ) storage comparable to the size of the input and output. However modern super computers are equipped with memory banks of size in tens of GB. This fostered us for developing multithreaded versions of Algorithm 1 and 2 in which we place loop over m outermost. That is, we introduced the multithread layer via OpenMP for m loop i.e., a set of m ∈ M i values (already assigned to i th MPI process) is evenly divided into a number of subsets equal to the number of threads mapped onto each physical core available on a given multi-core processor. We denote this subset of m values as M T i where subindex i is an i th thread. In such a way each physical core executes one thread concurrently and calculates intermediate results of SHTs for its local (in respect to shared memory) subset of m. The Algorithm 4 describes the step 1 of Algorithm 3 in its multithreaded version. 9) for every ring r do for every = m + 2, ..., max do • precompute β m Eq. ( 8)

Algorithm 4 ∆ A m calculation per thread Require: all m ∈ M T i ⇒ Nth i=1 M T i = M i for all m ∈ M T i do • initialise P mm Eqs. (
• compute P m via β m , Eq. ( 7)

• update ∆ A m (r) in shared memory, Eq. ( 12) end for ( ) end for(r) end for(m)

Load balance between active threads

Selecting a loop over m as an outermost one is certainly the best approach from the point of view of speed optimization (see previous subsection), but such configuration leads to poor load balance. Notice that due to assumptions listed in section 2.2 the -recursion starts from = m (see for example Eq. 12). Therefore, the number of flops required for evaluating the recurrence from equation ( 9) differs for different m. For that reason, we have to carefully map subsets M T i onto active threads to ensure that the total time required for the evaluation of the problem on a given node will be equally divided by the number of available cores.

Because the lenght of the loop over decreases with growing number of m, we create subsets M T i by taking values from the set M i in "max-min" fashion i.e., iteratively we take from M i the maximum value and the minimal one and we assign them to a chosen thread. In the same moment we try to equilibrate the number of such pairs within each subset M T i . Example of a such a mapping is depicted in Figure 5. In next subsection, we introduce a similar multithreaded approach for accelerating the projection of associated Legendre functions using GPU.

SHTs on GPGPU

A CUDA (Compute Unified Device Architecture) device is a chip consisting of hundreds of simple cores, called Streaming Processors (SP) which execute instructions sequentially. Groups of such processors are encapsulated in a Streaming Multiprocessor (SM) which executes instruction in a SIMD (Single Instruction, Multiple Data) fashion. Therefore each SM manages hundreds of active threads in a cyclic pipeline in order to hide memory latency i.e., a number of threads are grouped together and scheduled as a Warp, executing the same instructions in parallel. Therefore CUDA threads are grouped together into a series of thread blocks. All threads in the thread block execute on the same SM, and can synchronise and exchange data using very fast but small (up to 48KB in size for the latest NVIDIA carts) shared memory. The shared Figure 5: Example of two "min-max" pairs of m values assigned to Thread X on Process 1 (left) and 0 (right) respectively. The read arrow represents the direction and the lenght of the -recurrences (Eq. 10) evaluated by a thread for a given "min-max" pair. memory can be used as a user-managed cache enabling higher bandwidth than direct data copy from the slow device global memory.

The restrictions related to memory (both size and speed) and the lack of synchronisation between the blocks of threads drive the structure of our algorithms devised for GPUs. For instance, placing the loop over m as outermost (like in algorithms dedicated for multi-core CPUs) is not appropriate for CUDA applications due to the memory size limitation. Consequently, the vector β m cannot be entirely stored in shared memory. Its values would need to be recomputed for each m and accessed sequentially in the -loop. A possible implementation would require re-computing the β m coefficients for each pass through the -loop, which would significantly affect the performance. Therefore, to agree with the philosophy of programming for GPUs, which consists in using fine grained parallelism and launching a very large number of threads in order to use all the available cores and hide memory latency, we set iterations over r in the outermost loop and parallelize it in such a way that a number of rings (preferably only one) is assigned to one thread. In this model each thread computes the 2-point recurrence for all available m-values (but we assume that GPU is a part of distributed memory system therefore) in contiguous fashion as depicted on Figure 6. This makes it straightforward to plan the computation of β m via Eq. ( 9) in batches, which can be cached and shared inside a thread block. Note that in case of distributed memory systems the number of m values in set M i may vary for different i, while the total number of rings R N stays constant. Therefore, it is easier to plan good load balance for GPU when we parallize our algorithms over r loop. Setting loop over rings outermost in our nested loop triggers two completely different approaches for memory management in direct and inverse SHTs. Therefore we dedicate the following two subsections to describe in detail the differences between map2alm and alm2map algorithms for GPUs.

map2alm

Mapping rings onto big number of threads, indeed facilitates the parallelization of SHTs algorithm on GPUs i.e., in such set-up each active thread perform exactly the same number and type of operations but with different values. However, the direct transform requires that for each ring we update all non-zero values of the matrix a ,m ∈ C mmax× max . This leads to all-reduce operation6 between active threads. This is to say, the partial results generated by all active threads need to be summed together and written in the memory space from which we will read the final solution after all threads will be executed. Note that on CUDA, active threads can communicate only in thread block, which makes this task more difficult.

The Algorithm 5 outlines the Step 3 of algorithm 1 adapted to CUDA. We assume however that we are working on a memory distributed system, therefore the range of m values is limited to M i . To ensure workaround for the high latency device memory and take advantage of the fast shared memory, we first calculate the values of β m vectors in segments, which allows us to fit them to the shared memory i.e., all threads within a block in a "collaborative" way compute in advance β m values and store them in shared memory. Next, we use those precomputed values to evaluate associated Legendre function P m (r). After, each active thread fetches from global memory a precomputed sum ∆ S m (r) and performs a final multiplication Eq. ( 13). Once it's done, we sum all those partial results within a given block of threads using adapted7 version of parallel reduction primitive from CUDA SDK. For such locally reduced result one thread in block updates corresponding entry in global memory. This implies that j = NBT (number of blocks of threads) threads will try to increase the value of a variable at the same space of global memory -which is referred to as a race condition problem, common in multithread applications. Fortunately, race conditions are easy to avoid in CUDA by using atomic operations. An atomic operation is capable of reading, modifying, and writing a value back to memory without the interference of any other threads, which guarantees that a race condition will not occur. Therefore, all final updates of final results in global memory are performed via CUDA atomicAdd function. In case when atomic functions are not supported on a given GPU we can easily avoid race conditions by increasing granularity of recurrence step. In this scenario we loop over m values on CPU and execute M i "smaller " CUDA kernels which for a fixed values m and µ m , evaluate NBT reduced (per block of threads) results a m , which are next copy from device to host memory and sum together using CPU. This variant of our algorithm corresponds to the state HYBRID_EVALUATION = true, in listing 5.

alm2map

The inverse SHT also requires the reduction of partial results, but in contrary to the direct transform, we need to sum them for all as displayed in Eq. [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. With rings assigned to single threads, this leads to rather straightforward implementation. Algorithm 6 outlines the alm2map algorithm devised for GPUs and it is based on the version detailed in [START_REF] Hupca | Spherical harmonic transform with gpus[END_REF], which conforms with all assumptions stated in this paper.

First, we partition the vector of the a m coefficients in the global memory into small tiles. In this way a block of threads can fetch resulting segments in sequence and fit them in the shared memory. Next, we calculate β m values in the same way as in the map2alm algorithm. With all the necessary data in the fast shared memory, in following steps, each thread evaluates in parallel an associated Legendre function and updates the partial result ∆ A m . In contrast to the direct transform, those updates are performed without a race condition problem i.e., each thread stores its temporary values in small space of the shared memory to which it has an exclusive and very fast access. Once, the loop over all is finished, the final partial result is written to the global memory.

Exploitation of hybrid architectures

Many new supercomputers are hybrid clusters i.e., they are a combination of computational nodes interconnected (usually via PCI-Express bus) with multiply GPUs. In such a way, one or more multi-core processors have access to one GPU. Facing these heterogeneities, the front end user of the spherical harmonic package will have to decide which type of architecture to choose for calculation. The three-stage structure of our algorithms allows us to freely map different steps of the algorithms onto different accelerators in hybrid system. However, the steps have to be evaluated sequentially in unchanged order. The intention of this flexibility is to facilitate selection of the most efficient library for FFTs e.g. NVIDIA provides its own library for FFT on GPUs i.e., CUFFT [21] which efficiency strongly depends on the version of CUDA device and the version of the routine. Therefore, on some supercomputers with highly tuned version of FFTW3 library, FFTs perform better on multi-core CPUs. Note that our algorithms require the evaluation of FFTs on complex data, which length may vary from ring to ring.

In practice, on heterogeneous systems we assign each step of the algorithm to the architecture on which it performs faster. The default assignment for both transformation is depicted in Figure 7

A m ∈ C
to their sequential nature, we omit steps in which we evaluate the starting values of the recurrence Eq. ( 9) (step 2 in Algorithm 1 and step 1 in Algorithm 2), which are always computed on CPU.

Figure 7: Default assignment of the major steps of the SHTs algorithm to one of the two architectures, multi-core processor or GPU, in a GPU/CPU heterogeneous system. Step 3 in map2alm algorithm covers two architectures due to optional use of CPU for reducing partial results between the execution of the CUDA kernels.

Numerical experiments

In order to evaluate the accuracy, efficiency and scalability of our new algorithm, we perform a series of tests with different geometry of the 2-sphere grid and different values of max and m max . To validate our algorithms, in each experiment we perform a pair of backward and forward SHTs starting out with a backward transform on a set of random a init m coefficients which are uniformly distributed random numbers in the range (-1, 1). To measure the accuracy of the algorithm we analyse the resulting map and we calculate the error as the difference between the final output a out m and the initial set a init m via the following formula,

D err = m a init m -a out m 2 m a init m 2 . ( 19
)
The pairs of transforms were tested on HEALPix grids with a different resolution parameter N side . Accordingly to this parameter n pix = 12N 2 side , R N = 4N side -1 where maximum lenght of the ring (equator) is equal to 4N side . Finally we also choose max = m max = 2 × N side , except for some specific experiments.

Validation against other implementation

In this experiment, a set of a init m coefficients with a different value of max was generated, then a pair of backward and forward transforms were performed to converted a init m to a HEALPIx map with a fixed N side = 1024 and then back to a out m . Starting with the same input array we used our both implementation (MPI/OpenMP & MPI/CUDA) and the libpsht library [START_REF] Reinecke | Libpsht -algorithms for efficient spherical harmonic transforms[END_REF] to evaluate transforms. We measured the total time and the final error D err Eq. [START_REF]MPICH -A portable implementation of MPI[END_REF].

Figure 8 depicts D err errors for s2hat_{mt,cuda} and libpsht transforms pairs on HEALPIx grids with various max and N side values. Every data point is the maximum error value encountered in the pairs of transforms. As we can observe, all implementations provide almost identical quality of solutions with respect to different sizes of the problem.

Scaling with growing resolution and performance comparison with libpsht

To study the time scaling of the algorithm with growing resolution of the problem defined via N side parameter, we perform a number of experiments with various N side = max /2 parameter. Like in previous experiment we compare our results with libpsht package.

We choose Martin Reinecke's libpsht because it is to our knowledge, the fastest implementation of SHTs suitable for an astrophysical application and it is essentially a realisation of the same algorithm as the one used in our parallel algorithm. This highly optimized implementation is based on explicit multithreading and SIMD extensions (SSE and SSE2) but its functionality is limited to the shared memory systems only. Consequently, we perform all comparison experiments with only one pair of multi-core processor and GPU on three different systems listed in Table 2. In a similar fashion as in our multithreaded version, all CPU-intensive parts of libpsht's transforms algorithm are instrumented for parallel execution on multi-core computers using OpenMP directives. Therefore, in all our tests we set-up for all runs a number of OpenMP threads equal to the number of physical cores (OMP_NUM_THREADS = 4).

The range of pixelization sizes was limit by the size of the memory of a single GPU. Therefore we limit our tests to N side = 4096. We also limit timing measurement in this experiment to the recurrence part only (step 3 in Algorithm 1 and step 2 in Algorithm 2), however in case of CUDA routines, those measurements includes also the time spend on any data movements between GPU and CPU.

The results of this experiment are shown in Figures 9 and10. The figures depict the measured times in seconds and the number of Gflops per second in logarithmic scale. The number or flops include all the operations necessary for the purpose of the numerical stability like rescaling in case of s 2 hat routines (see paragraph 2).

As expected, our multithreaded version is always slower in comparison to libpsht by a constant factor. Mainly, due to lack of SSE intrinsics in our code, which allows on some hardware to perform two or more arithmetic operations of the same kind with a single machine instruction. Especially on latest intel i7-2600K processor, SSE optimised libpsht implementation takes advantage of the extended width of the SIMD register 8 and Advanced Vector Extensions (AVX) which significantly increase parallelism and throughput in floating point SIMD calculations.

The other important observation is that forward and backward SHTs on CPU, with identical parameters, require almost identical time. Furthermore, performance analysis allows several deductions about performance of GPU version of SHTs algorithms.

• It is fairly evident that forward transform on CUDA with slower, not only with respect to its CPU version, but also in comparison with backward transform on the same architecture. Since the implementation of the CUDA kernel for realisation of projection onto associated Legendre polynomials is essentially the same for both directions of the transform, the performance lost is mainly due to very 2 for details) pairs of Intel multicore processor with GPU on a HEALPIx grid. The top pane shows timings for projection on Legendre polynomials, the libpsht code (orange circles) is compared with hybrid s 2 hat codes (blue&green lines). In each case we use a HEALPIx grid with N side = max /2. In case of CUDA routines, timings includes time of data movements between host and device. Bottom pane shows number of GFLOPS per second. Note that GeForce 400 Series equipped with latest CUDA architecture (codename "Fermi") outperform libpsht on the latest, Intel Core i7-2600K. Bear in mind that GeForce family is not design for professional use therefore their double-precision throughput has been reduced to 25% of the full design.

expensive reduction of partial results. First in shared memory and then on slow global memory or on host memory with use of CPU.

• Thanks to the latest CUDA architecture (codename "Fermi") GeForce 400 Series of GPUs by far dominate over old generation Tesla S1070. Note that, unlike Tesla series, GeForce family is not designed for professional use therefore their double-precision throughput has been reduced to 25% of the full design. Nevertheless, those two GPUs (GTX 460 & GTX 480) outperform libpsht on the latest intel i7 processor in the evaluation inverse transform.

• CUDA version seems to favour cases with big N side which indicates usage of many threads (as we explained in §4.2.2, each ring of the sphere is assigned to one CUDA thread) i.e., both GPU routines gain performance with growing number of sky-pixels (thus, also rings).

Scaling with number of multi-core processors and GPUs

The main purpose of using parallel SHTs is to exhibit data parallelism i.e., to allow us perform SHTs in big resolutions by distributing our input and output data across different computing nodes to be processed in parallel. A good parallel algorithm should scale with a size of the problem and a number of computing nodes. To measure those two quantities we perform a series of experiments with different values of N side and a number of MPI processes spanned on hybrid pairs of multi-core CPU-GPU. For all experiments in this section, we set max = 2 × N side . All experiments in this section were executed on CCRT's super computer, one of the first hybrid cluster built, based on Intel and Nvidia processors. It has 1068 nodes dedicated for computation, each node is composed of 2 Intel Nehalem quad-cores (2.93 GHz) and 24 GB of memory. One part of the system has a Multi-GPU support with 48 Nvidia Tesla S1070 servers, 4 Tesla cards with 4GB of memory and 960 processing units are available for each server. Two computation nodes are connected to one Tesla server by the PCI-Express bus, therefore 4 processors handle 4 Tesla cards. The high performance interconnect is based on Infiniband DDR. This system provides a peak performance of 100 Tflops for the Intel part and 200 Tflops for the GPU part. Intel compiler version 12.0.3 and NVIDIA Cuda compiler 4.0 were used for compilation. Furthermore, we span MPI processes in our experiment in such a way that exactly one MPI is assigned to one pair of multi-core processor and GPU (4 MPI processes per one Tesla server). To increase intra-node parallelism we set-up for all runs a number of OpenMP threads equal to the number of physical cores of Nehalem CPU (OMP_NUM_THREADS = 4).

Figure 11 shows maximum overall time measured on one of 128 CPU/GPUs involves in computation of pair of SHTs in different size of the map. As we can observe, our both implementations scale well with respect to the size of the problem. As expected, only the communication cost grows slowly with number of processes. Note, that according to our prediction (see §2.3), the evaluation of the projection onto associated Legendre polynomials dominates by far over the remaining steps and it de facto defines the overall time. On contrary, the time spend on memory transfers between host and the GPU device never exceed 1s, and deviation from perfect scalability visible on plots can be explained by random errors of MPI_Wtime() function used in our experiments for time measurement (due to its limited resolution, MPI_Wtime() rerun trustable results only for elapsed times greater then 1s). In both implementations we assign the evaluation of FFTs to the CPU due to its slightly better performance in comparison to CUFFT also available on Titane (see Fig. 14).

SHTs accuracy on CUDA

Both x86 and CUDA architecture used in our experiments respect the IEEE Standard for Binary Floating-Point Arithmetic (IEEE 754-1985). IEEE 754 standardizes how arithmetic results should be approximated in floating point. Whenever working with inexact results, programming decisions can affect accuracy. It is important to consider many aspects of floating point behaviour in order to achieve the highest performance with the precision required for any specific application. This is especially true in a heterogeneous computing environment where operations are perform on different types of hardware. The IEEE 754 standard defines a handful of basic arithmetic operations like add, subtract, multiply, divide and square root. The results of these operations are guaranteed to be the same for all implementations of the standard, for a given format and rounding mode. However, IEEE standard do not define accuracy of math functions like sinus, cosines or logarithm, which depends on implementation. For instance CUDA math library includes all the math functions listed in the C99 standard plus some additional useful functions [START_REF]NVIDIA CUDA Programming Guide[END_REF]. These functions have been tuned for a reasonable compromise between performance and accuracy which may differ from their equivalent on x86 architectures. For that reason, we performed series of experiment in which we converted the same random set a init m (see beginning of this section) to a HEALPIx map on two different architecture, namely NVIDIA GPU and x86 processor, and we calculate an absolute difference of those two maps (see the example result in Figure 16). In all our tests the maximum difference between two sky-pixels never exceeded value of 1.0 • 10 -11 , which indicates a very good agreement between those two architectures.

Conclusion and future work

This paper describes a parallel algorithm for computing the SHTs with two variants of intra-node parallelism specific for given architectures i.e., multi-core processors and GPUs. We tested both versions in terms of accuracy, overall efficiency and scalability. We showed that our inverse SHTs with GeForce 400 Series INRIA equipped with latest CUDA architecture (codename "Fermi") outperform state of the art implementation for multi-core processor executed on latest Intel Core i7-2600K.

For solving big problems, a cluster of NVIDIA Tesla S1070 can accelerate overall execution of inverse transform by as much as 3 times with respect to the hybrid MPI/OpenMP version executed on the same number of quad-core processors Intel Nahalem. We also showed that this nontrivial algorithm exhibits some difficulties in full adaptation to CUDA architecture i.e., a direct transform is not as well adapted to GPUs as an inverse transform, which impacts significantly overall performance. Nevertheless, direct transform with NVIDIA GTX 480 outperformed libpsht on Core i7-960K.

The current implementation still leaves some space for future optimizations. Therefore, as a future work, we may add SSE intrinsics ( or/and Advanced Vector Extensions (AVX) for latest x86 architectures) in order to introduce low level vectorization on CPUs. From the algorithm side, it would be crucial to enhance the set of our routines to support transforms with spin quantities, which are important part of a CMB science.

Figure 2 :

 2 Figure 2: The overall time breakdown between the main steps of the SHTs algorithm as computed with the identical band limit ( max = 4096) for quad-core Intel i7-2600K processor.

Algorithm 2

 2 Inverse Transform (Eq. 4) Require: a m ∈ C ×m values step 1 -Pre-computation for every m = 1, ..., mmax do • µm precomputation Eq. (9) end for (m) step 2 -∆ A m calculation for every ring r do for every m = 0, ..., mmax do for every = m, ..., max do • compute P m via Eq. (7) • update ∆ A m (r), Eq. (14) end for ( ) end for (m) end for (r) step 3 -s calculation for every ring r do for every m = 0, ..., mmax do • calculate s via FFT and given ∆ A m (r), Eq. (11) end for(m) end for (r) return s ∈ Rn pix

  INRIA (a)The iso-latitudinal rings distribution.(b) Distribution of a m coefficients.

Figure 3 :

 3 Figure 3: Example of a distribution of the iso-latitudinal rings and the a m coefficients among two processes. Each colour marks elements assigned to one process. For rings visualisation we used an August's Conformal Projection of the sphere on a two-cusped epicycloid.

Algorithm 3

 3 General parallel alm2map algorithm ( Code executed by each MPI process ) Require: a m ∈ C ×m for m ∈ M i and all Require: indices of rings r ∈ R i step 1 -Pre-computation for every m ∈ M i do • µ m pre-computation Eq. (9) end for (m) step 2 -∆ A m calculation for every ring r ∈ R and every m ∈ M i do

17 )Figure 4 :

 174 Figure 4: The two left most plots show theoretical run-times (logarithmic scale) as function of a number of processes and a fixed size of the problem, left most panel, and as a function of a resolution/size of the problem for the fixed number of processes, middle. The sizes of the problem shown in the leftmost panel with thick, solid lines, correspond to 2 10 8 pixels and corresponding harmonic modes, as given the HEALPix parameter value N side = 4096. Thin, dashed lines show cases with 4 times fewer (more) pixels, lower and upper lines respectively. In the middle panel the number of processes is assumed to be n proc = 512, thick, solid line, or 256 (1024), as shown by thin dashed lines. The contours in the rightmost panel display a ratio of the total computation to communication times shown as a function of a number of process and size of the problem.

Figure 6 :

 6 Figure 6: Schematic view of -recurrence iterations (eq. 10) evaluated by each thread on CUDA device assigned to Process 1 and 0 respectively.

Figure 8 :

 8 Figure 8: D err for s2hat_{mt,cuda} and libpsht transforms pairs on HEALPIx grids with various max and N side . Every data point is the maximum error value encountered in transforms pairs. Loss of precision seen in the left panel for max > ∼ 2048 and shared by both implementations is due to aliasing inherent to the adopted pixelization scheme. CPU clock speed GPU Compilators line style Core i7-960K 3.20 GHz GeForce GTX 480 gcc 4.4.3 / nvcc 4.0 solid line Core i7-2600K 3.40 GHz GeForce GTX 460 gcc 4.4.5 / nvcc 4.0 dashed line Xeon X5570 2.93 GHz Tesla S1070 icc 12.1.0 / nvcc 4.0 dot line

Figure 9 :

 9 Figure 9: Benchmarks for direct SHT performed on different (see Table2for details) pairs of Intel multicore processor with GPU on a HEALPIx grid. The top plot shows timings for projection on Legendre polynomials, the libpsht code (orange circles) is compared with hybrid s 2 hat codes (blue&green lines). In each case we use a HEALPIx grid with N side = max /2. In case of CUDA routines, timings include time of data movements between CPU and GPU. Bottom pane shows number of GFLOPS per second. Note that CUDA routines loose performance in comparison to the SSE optimised lipsht due to very costly reduction of partial results. First in shared memory and then on host memory with use of CPU. Nevertheless, NVIDIA GTX 480 outperforms libpsht on Core i7-960K.

Figure 10 :

 10 Figure 10: Benchmarks for inverse SHT performed on different (see Table2for details) pairs of Intel multicore processor with GPU on a HEALPIx grid. The top pane shows timings for projection on Legendre polynomials, the libpsht code (orange circles) is compared with hybrid s 2 hat codes (blue&green lines). In each case we use a HEALPIx grid with N side = max /2. In case of CUDA routines, timings includes time of data movements between host and device. Bottom pane shows number of GFLOPS per second. Note that GeForce 400 Series equipped with latest CUDA architecture (codename "Fermi") outperform libpsht on the latest, Intel Core i7-2600K. Bear in mind that GeForce family is not design for professional use therefore their double-precision throughput has been reduced to 25% of the full design.

Figure 11 :

 11 Figure 11: Maximum elapsed wall clock time for a pair of transforms with various N side and max = 2×N side . Experiment evaluated with 128 MPI processes spanned on the same number of pairs Intel Nehalem quadcores and Nvidia Tesla S1070.

Figures 12

 12 Figures 12 and 13 depict the runtime distribution over main steps of the algorithm multiplied by the number of MPI processes. In both versions, all steps scale almost perfectly (at least up to 128 MPI processes). As expected, only the communication cost grows slowly with number of processes. Note, that according to our prediction (see §2.3), the evaluation of the projection onto associated Legendre polynomials dominates

Figure 12 :

 12 Figure 12: Maximum elapsed wall clock time distribution over main steps of the algorithm for the hybrid MPI/OpenMP version of the alm2map_mt and map2alm_mt routines. Results for experiment with N side = 16384 and max = 2 × N side .

Figure 13 :

 13 Figure 13: Maximum elapsed wall clock time distribution over steps of the algorithm for the hybrid MPI/CUDA version of the alm2map_cuda and map2alm_cuda routines. Results for experiment with N side = 16384 and max = 2 × N side .

Figure 14 :

 14 Figure 14: The overall difference between the FFTs evaluation on Titane supercomputer via cuFFT on GPU and Intel MKL Fourier Transform wrapped to FFTW3 interface. Values correspond to experiment with N side = 16384 and max = 2 × N side .Finally, if we compare overall time between CPU and GPU version, we find out that implementation of alm2map for CUDA is in all our big test cases, in average ×3.3 faster then multithreaded version. In contrary, map2alm due to time consuming reduction operations remains slower then its counterpart dedicated for CPUs. Average speed-up of hybrid implementation MPI/CUDA with respect to hybrid MPI/OpenMP version is depicted in Figure15

Figure 15 :

 15 Figure 15: Average speed-up of hybrid implementation MPI/CUDA with respect to hybrid MPI/OpenMP version for different number of MPI processes.

  (a) Reference map evaluated on CPU (b) Map generated on Nvidia Tesla S1070 GPU (c) Absolute difference between maps generated on CPU and Nvidia GPU.

Figure 16 :

 16 Figure 16: CMB maps generated on two different architectures and their absolute difference.

Table 1 :

 1 Number of floating point operations required in the parallel SHTs algorithm.

  . Note that due Algorithm 5 a m calculation on CPU&GPU Require: ∆ S m & µ m vectors in GPU global memory for m ∈ M i for every ring r assign to one GPU thread do for every m ∈ M i do for every = m + 2, ..., max do • use precomputed or, if needed, precompute in parallel a segment of β m (8); • copy partial results a GP U m [. . .] from device to host memory and sum together using CPU for i = 0 → NBT do a m = a m + a GPU m [i] end for else • update a m in global memory via atomic addition: a m = a m + a m ,m ∈ C mmax× max matrix in GPU global memory for every ring r assign to one GPU thread do for m ∈ M i do for every = m + 2, ..., max do • use fetched or, if needed, fetch in parallel a segment of a m data;

	• compute P m via Eq. 7;
	• evaluate: a r m = ∆ S m (r) P m (r) • sum all partial results a r m within a block of threads
	NTPB
	a m =	(a r m ) i	(18)
	i=0
	if HYBRID_EVALUATION then
	• for each block of threads with different id = BID, write a m in global device memory
	a GP U m [BID] = a m
	end if
	end for ( )
	end for(m)
	end for(r)
	return a ,m ∈ C mmax× max
	Algorithm 6 ∆ A m calculation on GPU
	Require: a • use precomputed or, if needed, precompute in parallel a segment of β m (8);
	• compute P m via Eq. 7;
	• update ∆ A m (r) for a given prefetched a m and computed P m (12);
	end for ( )
	• save final ∆ A m (r) in global memory
	end for(m)
	end for(r)
	return vector of partial results: ∆

planck: http://sci.esa.int/science-e/www/area/index.cfm? fareaid =17

cobe: http://lambda.gsfc.nasa.gov/product/cobe/

wmap: http://map.gsfc.nasa.gov/ RR n°7635

Wavemoth: https://github.com/wavemoth/wavemoth INRIA

RR n°7635

The IEEE has standardised the computer representation for binary floating-point numbers in IEEE 754. This standard is followed by almost all modern machines.

Reduce-type operation combines the distributed elements using predefined operation (e.g., summation, multiplication etc.) and returns the combined value in the output buffer.

We omitted a part of the code which writes partial result to global memory RR n°7635

The width of the SIMD register in file on i7-2600K is increased from 128 bits to 256 bits in respect with previous generation of i7 processors

Acknowledgments

This work has been supported in part by French National Research Agency (ANR) through COSINUS program (project MIDAS no. ANR-09-COSI-009). The HPC resources were provided by GENCI-[CCRT] (grants 2011-066647 and 2012-066647) in France and by the NERSC in the US, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Some of the results in this paper have been derived using the HEALPix package [13].