Parallel spherical harmonic transforms on heterogeneous architectures (graphics processing units/multi‐core CPUs)

Abstract

SUMMARY 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 SHT with two variants of intra‐node parallelism appropriate for novel supercomputer architectures, multi‐core processors and Graphic Processing Units (GPU). It also discusses their performance, alone and embedded within a top‐level, Message Passing Interface‐based parallelisation layer ported from the S2HAT library, in terms of their accuracy, overall efficiency and scalability. We show that our inverse SHT run on GeForce 400 Series GPUs equipped with latest Compute Unified Device Architecture 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 Message Passing Interface/Compute Unified Device Architecture version of the inverse transform run on a cluster of 128 Nvidia Tesla S1070 is as much as 3 times faster than the hybrid Message Passing Interface/OpenMP version executed on the same number of quad‐core processors Intel Nehalem for problem sizes motivated by our target applications. Performance of the direct transforms is however found to be at the best comparable in these cases. We discuss in detail the algorithmic solutions devised for the major steps involved in the transforms calculation, emphasising those with a major impact on their overall performance and elucidates the sources of the dichotomy between the direct and the inverse operations.Copyright © 2013 John Wiley & Sons, Ltd.

High-precision simulations of the weak lensing effect on cosmic microwave background polarization

“…The major bottleneck of the code performance is due to the need of calculating a single inverse spherical harmonic transform which is required to obtain the overpixelized map of the unlensed signal. This can certainly be alleviated further by using better algorithms and/or numerical implementations, e.g., capitalizing on hardware accelerators such as GPGPU (Hupca et al 2012;Szydlarski et al 2011;Fabbian et al 2012;Reinecke & Seljebotn 2013). We leave these code optimizations for future work.…”

“…It has a nearly perfect memory scalability obtained via a memory distribution of all main pixel and harmonic domain objects (i.e., maps and harmonic coefficients), and ensures very good load balance from the memory and calculation points of view. It is a very flexible tool that allows a simultaneous, multimap analysis of any iso-latitude pixelization, symmetric with respect to the equator, with pixels equally distributed in the azimuthal angle, and provides support for a number of pixelization schemes, including the above mentioned ECP; see Szydlarski et al (2011) for more details. The core of the library is written in F90 with a C interface and it uses the message passing interface (MPI) to institute distributed memory communication, which ensures its portability.…”

“…As can also be seen from the right panel of Fig. 14, the runtime connected to the inverse spherical harmonic transform of the unlensed sky, despite being perfectly load-balanced, tends to flatten due to required internal communication steps and precomputations to initialize the recurrence to compute spherical harmonics. These are per se subdominant steps, but they are expected to play a role for a very fine parallelization (Szydlarski et al 2011). The overall remapping procedure of pixel values requires a number of operation of about O(N pix /n) and is subdominant, since it operates on a lowerresolution map, and perfectly scalable because it does not require any communication.…”

Detecting the tensor-to-scalar ratio with the pure pseudospectrum reconstruction ofB-mode

“…In the case of azimuthally symmetric patches, the numerical computation of such Fisher matrices (either F or F rr ), can be performed in a reasonable time using the expression found in the appendix F of Ref. [22], and by using the s 2 hat package to perform spherical harmonic transforms [41][42][43][44]. (The use of this massively parallel package allows for a rapid computation of the covariance matrix for large sky coverages.).…”

Iterative destriping and photometric calibration forPlanck-HFI, polarized, multi-detector map-making

“…For each noise realization, we constructed the residual map subtracting the signal+white noise map from the destriped map. We then computed full-sky angular power spectra for each residual map using the Xpol (Tristram 2006) code, based on the S 2 HAT library 5 (Hupca et al 2010;Szydlarski et al 2011) and averaged them over the simulated realizations. First we investigate the effect of large signal variations within the pixels on the offset determination.…”