cunuSHT: GPU Accelerated Spherical Harmonic Transforms on Arbitrary Pixelizations

Sebastian Belkner\orcidlinkhttps://orcid.org/0000-0003-2337-2926 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA Université de Genève, Département de Physique Théorique et CAP, 24 Quai Ansermet, CH-1211 Genève 4, Switzerland    Adriaan J. Duivenvoorden \orcidlinkhttps://orcid.org/0000-0003-2856-2382 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA    Julien Carron\orcidlinkhttps://orcid.org/0000-0002-5751-1392 Université de Genève, Département de Physique Théorique et CAP, 24 Quai Ansermet, CH-1211 Genève 4, Switzerland    Nathanael Schaeffer \orcidlinkhttps://orcid.org/0000-0001-5206-3394 ISTerre, Université de Grenoble 1, CNRS, F-38041 Grenoble, France    Martin Reinecke Max-Planck Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany
(20th June 2024)
Abstract

We present cunuSHT \faGithub, a general-purpose Python package that wraps a highly efficient CUDA implementation of the nonuniform spin-00 spherical harmonic transform. The method is applicable to arbitrary pixelization schemes, including schemes constructed from equally-spaced iso-latitude rings as well as completely nonuniform ones. The algorithm has an asymptotic scaling of 𝒪(max3)𝒪superscriptsubscriptmax3\mathcal{O}{(\ell_{\rm max}^{3})}caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ) for maximum multipole maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT and achieves machine precision accuracy. While cunuSHT is developed for applications in cosmology in mind, it is applicable to various other interpolation problems on the sphere. We outperform the fastest available CPU algorithm by a factor of up to 5555 for problems with a nonuniform pixelization and max>4103subscriptmax4superscript103\ell_{\rm max}>4\cdot 10^{3}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT > 4 ⋅ 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT when comparing a single modern GPU to a modern 32323232-core CPU. This performance is achieved by utilizing the double Fourier sphere method in combination with the nonuniform fast Fourier transform and by avoiding transfers between the host and device. For scenarios without GPU availability, cunuSHT wraps existing CPU libraries. cunuSHT is publicly available and includes tests, documentation, and demonstrations.

keywords: Nonuniform Spherical Harmonic Transform, Nonuniform Fast Fourier transform, Cosmic Microwave Background Weak Lensing, CUDA

I Introduction

Spherical harmonic transforms (SHTs) are a key ingredient in signal processing for data sets on the 2-sphere. They are extensively used in active research fields such as cosmology (both in studying the cosmic microwave background (CMB) [1, 2, 3, 4, 5, 6] and large-scale structure of the Universe [7, 8]), gravitational waves [9], meteorology [10], solar physics [11], or solving partial differential equations on the sphere [12].

Modern data sets routinely require the evaluation of the SHT up to max=𝒪(104)subscriptmax𝒪superscript104\ell_{\rm max}=\mathcal{O}(10^{4})roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = caligraphic_O ( 10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT ) for maps with Nmax2similar-to𝑁superscriptsubscriptmax2N\sim\ell_{\rm max}^{2}italic_N ∼ roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT pixels. Here maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT denotes the largest multipole \ellroman_ℓ considered in the SHT. A direct evaluation of the SHT scales as 𝒪(max4)𝒪superscriptsubscriptmax4\mathcal{O}(\ell_{\rm max}^{4})caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT ), which is intractable for large problems. Reducing the computational complexity is thus crucial. One standard optimization is achieved by pixelizing the sphere into rings of constant latitude with equi-angular spaced samples, reducing the problem to 𝒪(max3)𝒪superscriptsubscriptmax3\mathcal{O}(\ell_{\rm max}^{3})caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ). We will refer to this setup as the ring spherical harmonic transform (rSHT).111There exist algorithms for the rSHT that asymptotically scale as 𝒪(max2log2(max))𝒪superscriptsubscriptmax2superscript2subscriptmax\mathcal{O}{(\ell_{\rm max}^{2}\log^{2}(\ell_{\rm max}))}caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_log start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ) ) [13, 14, 15, 16, 17] or 𝒪(max2log2(max)/loglog(max))𝒪superscriptsubscriptmax2superscript2subscriptmaxsubscriptmax\mathcal{O}{(\ell_{\rm max}^{2}\log^{2}(\ell_{\rm max})/\log\log(\ell_{\rm max% })})caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_log start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ) / roman_log roman_log ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ) ) [18]. While there have been improvements over the last few years [19, 20] such methods require high memory usage and significant pre-computations. A general purpose implementation that is competitive with the 𝒪(max3)𝒪superscriptsubscriptmax3\mathcal{O}(\ell_{\rm max}^{3})caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ) counterpart has yet to be published. See [21, 22] for modern implementations of the rSHT. For cases where the transform has to be evaluated on irregularly sampled pixels, ring sampling is not possible. We will refer to this more general setup as the nonuniform spherical harmonic transform (nuSHT). Notable applications of the nuSHT are CMB weak lensing [23], ray tracing [24, 25], or fields where the pixelization changes over time. It should be noted that due to their computational complexity, the rSHT and nuSHT often become the bottleneck in iterative algorithms that repeatedly apply the transforms [26, 27].

In the field of cosmology, nuSHTs are routinely solved via an rSHT and subsequent interpolation to the nonuniform points using bicubic splines [28] or a Taylor series expansion [29]. These methods scale as 𝒪(max3)𝒪superscriptsubscriptmax3\mathcal{O}{(\ell_{\rm max}^{3}})caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ) but only reach relatively low accuracy. A different 𝒪(max3)𝒪superscriptsubscriptmax3\mathcal{O}{(\ell_{\rm max}^{3}})caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ) nuSHT method, proposed by [30], makes use of the double Fourier sphere (DFS) method [31] in combination with the nonuniform fast Fourier transform (nuFFT) to achieve an accurate nuSHT algorithm. Several other implementations of this setup exist [32, 16]. Recently, the lenspyx and DUCC libraries [21, 33] implemented a highly efficient and machine precision accurate implementation of the nuSHT based on the DFS method.

In recent years, there has been a dramatic increase in the availability and usage of graphics processing units (GPUs) dedicated to scientific computing. This development drives the need for GPU based codes and provides an opportunity for increased performance of existing methods. GPUs are optimized for single-instruction multiple-data (SIMD) applications and provide multi-threading well beyond what is achievable with CPUs. Highly parallelizable algorithms can thus greatly benefit from the GPU architecture. In principle, the evaluation of the SHT allows for a large amount of parallel computation, making it a natural target for a GPU implementation. The work presented in [34, 35] was one of the first that explored the use of GPUs for SHTs in the context of cosmology.

Robust and efficient rSHT GPU implementations have been developed in recent years. Notable examples are SHTns [22], supporting the spin-00 and 1111 rSHTs, and S2HAT222https://apc.u-paris.fr/APC_CS/Recherche/Adamis/MIDAS09/software/s2hat/s2hat.html [35], which supports spin-n𝑛nitalic_n transforms. S2FFT [36] is a recent JAX implementation of the spin-n𝑛nitalic_n rSHT that provides differentiable transforms. A first implementation of the nuSHT in a cosmological context on GPUs was presented in [37].

We present cunuSHT, a CUDA accelerated nuSHT algorithm on the GPU. This is, to our knowledge, the first publicly available nuSHT GPU algorithm that reaches machine precision accuracy and achieves significant speed-up compared to the fastest CPU algorithms. We achieve this by carefully combining existing robust and efficient GPU implementations of the rSHT and nuFFT algorithm. cunuSHT does not require memory allocation or calculation on the host, which allows it to be incorporated in GPU-based software.

The remainder of the paper is organized as follows. In Section II we introduce notation and definitions, and present in qualitative terms our implementation. Section III discusses the implementation on the GPU. Section IV shows benchmarks and results. We conclude in Section V. A series of appendices collects further details.

II Nonuniform Spherical Harmonic transform

We introduce our notation and conventions in  II.1, and define the spherical harmonic transform operations that we implement in this paper. We describe the double Fourier sphere method in II.2.

II.1 Definition and properties

The (spin-0) spherical harmonic functions Ym(θ,ϕ)superscriptsubscript𝑌𝑚𝜃italic-ϕY_{\ell}^{m}(\theta,\phi)italic_Y start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT ( italic_θ , italic_ϕ ), with quantum numbers \ellroman_ℓ and m𝑚mitalic_m, with m𝑚-\ell\leq m\leq\ell- roman_ℓ ≤ italic_m ≤ roman_ℓ, are given by,

Ym(θ,ϕ)=Pm(θ)eimϕ,superscriptsubscript𝑌𝑚𝜃italic-ϕsuperscriptsubscript𝑃𝑚𝜃superscript𝑒i𝑚italic-ϕY_{\ell}^{m}(\theta,\phi)=P_{\ell}^{m}(\theta)e^{\mathrm{i}m\phi},italic_Y start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT ( italic_θ , italic_ϕ ) = italic_P start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT ( italic_θ ) italic_e start_POSTSUPERSCRIPT roman_i italic_m italic_ϕ end_POSTSUPERSCRIPT , (1)

where Pm(θ)superscriptsubscript𝑃𝑚𝜃P_{\ell}^{m}(\theta)italic_P start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT ( italic_θ ) are the associated Legendre polynomials. A “general” (or “nonuniform”) spherical harmonic transform (nuSHT) is a linear transformation between a set of spherical harmonic coefficients and field values defined at arbitrary locations on the sphere. We distinguish two types of transforms, with nomenclature inspired by nonuniform Fourier transform literature333The DUCC package uses the names adjoint__\__synthesis__\__general and synthesis__\__general for type 1 and type 2. [38]:

  • Type 1 (also “adjoint nuSHT”, the adjoint operation to type 2 below): given as input a set of N𝑁Nitalic_N values fisubscript𝑓𝑖f_{i}italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT, and N𝑁Nitalic_N locations (θi,ϕi)subscript𝜃𝑖subscriptitalic-ϕ𝑖(\theta_{i},\phi_{i})( italic_θ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_ϕ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ), desired are the coefficients cmsubscript𝑐𝑚c_{\ell m}italic_c start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT defined by,

    cm=i=1NfiYm(θi,ϕi),subscript𝑐𝑚superscriptsubscript𝑖1𝑁subscript𝑓𝑖superscriptsubscript𝑌absent𝑚subscript𝜃𝑖subscriptitalic-ϕ𝑖c_{\ell m}=\sum_{i=1}^{N}f_{i}\>Y_{\ell}^{\dagger m}(\theta_{i},\phi_{i})\,,italic_c start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT = ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_Y start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT † italic_m end_POSTSUPERSCRIPT ( italic_θ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_ϕ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) , (2)

    for \ellroman_ℓ up to some band-limit maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT. We want the result to match a target accuracy ϵitalic-ϵ\epsilonitalic_ϵ requested by the user.

  • Type 2: given as input a set of harmonic coefficients cmsubscript𝑐𝑚c_{\ell m}italic_c start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT up to some band-limit maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT, and a set of N𝑁Nitalic_N locations (θi,ϕi)subscript𝜃𝑖subscriptitalic-ϕ𝑖(\theta_{i},\phi_{i})( italic_θ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_ϕ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ), desired are the field values,

    fi==0maxm=cmYm(θi,ϕi),subscript𝑓𝑖superscriptsubscript0subscriptmaxsuperscriptsubscript𝑚subscript𝑐𝑚subscriptsuperscript𝑌𝑚subscript𝜃𝑖subscriptitalic-ϕ𝑖f_{i}=\sum_{\ell=0}^{\ell_{\rm max}}\sum_{m=-\ell}^{\ell}c_{\ell m}Y^{m}_{\ell% }(\theta_{i},\phi_{i})\,,italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = ∑ start_POSTSUBSCRIPT roman_ℓ = 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT end_POSTSUPERSCRIPT ∑ start_POSTSUBSCRIPT italic_m = - roman_ℓ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_ℓ end_POSTSUPERSCRIPT italic_c start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT italic_Y start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT ( italic_θ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_ϕ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) , (3)

    again respecting a target accuracy ϵitalic-ϵ\epsilonitalic_ϵ as requested by the user.

In matrix notation, type 2 may be written as,

𝐟=𝐘𝐜,𝐟𝐘𝐜\mathbf{f}=\mathbf{Y}\mathbf{c},bold_f = bold_Yc , (4)

where the vector 𝐜𝐜\mathbf{c}bold_c collects the harmonic coefficients, the vector 𝐟𝐟\mathbf{f}bold_f the output field values, and the entries of the matrix 𝐘𝐘\mathbf{Y}bold_Y are the spherical harmonics. Type 1 is represented by the adjoint matrix 𝐘=[𝐘t]superscript𝐘superscriptdelimited-[]superscript𝐘𝑡\mathbf{Y}^{\dagger}=\left[\mathbf{Y}^{t}\right]^{*}bold_Y start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT = [ bold_Y start_POSTSUPERSCRIPT italic_t end_POSTSUPERSCRIPT ] start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT. It is worth noting that type 1 is not the inverse to type 2, except in special cases. We also use the qualifiers type 1 and type 2 for the analogous nonuniform (or uniform) Fourier transforms, where the spherical harmonics and coefficients are replaced by their plane wave counterparts.

In typical applications, the total number of points N𝑁Nitalic_N is comparable to the squared band-limit, max2superscriptsubscriptmax2\ell_{\rm max}^{2}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. In this case the naive computational complexity of these operations is 𝒪(max4)𝒪superscriptsubscriptmax4\mathcal{O}(\ell_{\rm max}^{4})caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT ).

II.2 Double Fourier sphere method

We adopt the approach proposed in [33] and implement type 1 and type 2 transforms using the double Fourier sphere (DFS) method. In this approach, the matrix 𝐘𝐘\mathbf{Y}bold_Y of the type 2 nuSHT is decomposed into 4 matrices,

𝐘=𝐍𝐅𝐃𝐒.𝐘𝐍𝐅𝐃𝐒\mathbf{Y}=\mathbf{NFDS}\,.bold_Y = bold_NFDS . (5)

The final operation 𝐍𝐍\mathbf{N}bold_N is a nonuniform Fourier transform of type 2 to the given locations, and the role of 𝐅𝐃𝐒𝐅𝐃𝐒\mathbf{FDS}bold_FDS is to produce the needed Fourier coefficients.

The matrix 𝐒𝐒\mathbf{S}bold_S is an iso-latitude rSHT, that transforms the input harmonic coefficients onto an equi-spaced grid in both ϕitalic-ϕ\phiitalic_ϕ and θ𝜃\thetaitalic_θ, covering the entire sphere.

Refer to caption
Figure 1: Illustration of the doubling step in the DFS method. The upper half shows a spherical map calculated on a rectangular grid and is mirrored along the θ=π𝜃𝜋\theta=\piitalic_θ = italic_π axis. The mirrored image is split in half across the ϕitalic-ϕ\phiitalic_ϕ-direction in the center, and swapped. The result is a function on the torus, with a 2D ordinary Fourier series having the exact same band-limit as the spherical harmonic series of the original map defined on the sphere. This allows the use of efficient nonuniform Fast Fourier Transform for accurate interpolation.

𝐃𝐃\mathbf{D}bold_D is a “doubling” operation, that extends the range of θ𝜃\thetaitalic_θ from [0,π]0𝜋[0,\pi][ 0 , italic_π ] to [0,2π)02𝜋[0,2\pi)[ 0 , 2 italic_π ), see Fig. 1. The doubling is performed by extending the meridians across the south pole back up to the north pole. The essential point is that the resulting map, seen as a map on the doubly-periodic torus, has a standard Fourier series with exactly the same Fourier band-limit444This may be seen for example from the well-known Fourier representation of the Wigner d𝑑ditalic_d-matrices [39, 40], and using the relation Pm(θ)=dm0(θ)(2+1)/4πsubscriptsuperscript𝑃𝑚𝜃subscriptsuperscript𝑑𝑚0𝜃214𝜋P^{m}_{\ell}(\theta)=d^{\ell}_{m0}(\theta)\sqrt{(2\ell+1)/4\pi}italic_P start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT ( italic_θ ) = italic_d start_POSTSUPERSCRIPT roman_ℓ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_m 0 end_POSTSUBSCRIPT ( italic_θ ) square-root start_ARG ( 2 roman_ℓ + 1 ) / 4 italic_π end_ARG. as the spherical harmonic band-limit of the input array cmsubscript𝑐𝑚c_{\ell m}italic_c start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT. Finally, 𝐅𝐅\mathbf{F}bold_F is simply the standard 2D Fourier transform that produces the Fourier coefficients input to 𝐍𝐍\mathbf{N}bold_N from the doubled map.

The adjoint operator 𝐘superscript𝐘\mathbf{Y}^{\dagger}bold_Y start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT, or the type 1 nuSHT is, by definition,

𝐘=𝐒𝐃𝐅𝐍.superscript𝐘superscript𝐒superscript𝐃superscript𝐅superscript𝐍\mathbf{Y^{\dagger}}=\mathbf{S^{\dagger}D^{\dagger}F^{\dagger}N^{\dagger}}\,.bold_Y start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT = bold_S start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT bold_D start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT bold_F start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT bold_N start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT . (6)

𝐍superscript𝐍\mathbf{N}^{\dagger}bold_N start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT is a nuFFT of type 1 that produces Fourier frequencies from the input locations and field values. 𝐅superscript𝐅\mathbf{F}^{\dagger}bold_F start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT produces from these frequencies a 2D map on the torus. 𝐃superscript𝐃\mathbf{D}^{\dagger}bold_D start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT (the adjoint doubling matrix) effectively “folds” this doubled map. The resulting map is 2π2𝜋2\pi2 italic_π-periodic in the ϕitalic-ϕ\phiitalic_ϕ-direction, and the θ𝜃\thetaitalic_θ-direction goes again from 00 to π𝜋\piitalic_π. The map is then transformed to harmonic coefficients with 𝐒superscript𝐒\mathbf{S}^{\dagger}bold_S start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT, an iso-latitude type 1 rSHT.

III Implementation

We discuss the concrete GPU implementation. Readers interested in the CPU equivalent may consult [33].

A GPU is designed to efficiently apply a single instruction on multiple data (SIMD). On the hardware side, it achieves this with Streaming Multiprocessors (SMs) (at the order of 100), that contain a number of simple processors for arithmetic operations (at the order of 100) that execute “warps” of 32 threads in parallel.

On the software side, a GPU accelerated program is executed via a number of threads that are arranged in thread blocks. The GPU is responsible for distributing the thread blocks across the SMs. High throughput is achieved by overloading SMs with many threads as to hide data latency and by ensuring that memory is accessed in multiples of the warp size.

We differentiate between the GPU memory that is “close” to the processor units and can be accessed fast by the device, and host memory, that is managed by the host system of the GPU, and which is generally slow to access by the GPU. We show data transfer benchmarks in Appendix A. Our implementation avoids data transfer and usage of host memory altogether; intermediate results are kept in GPU memory. This is realized by cupy-arrays in combination with a C++-binding nanobind, [41], handily providing a nanobind-cupy interface.

Our implementation of the individual operators 𝐍𝐍\mathbf{N}bold_N, 𝐃𝐃\mathbf{D}bold_D, 𝐅𝐅\mathbf{F}bold_F, and 𝐒𝐒\mathbf{S}bold_S and their adjoints are realised as follows.

For the (adjoint) synthesis (𝐒superscript𝐒\mathbf{S}^{\dagger}bold_S start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT) 𝐒𝐒\mathbf{S}bold_S, we use the highly efficient software package SHTns [22], and calculate the SHTs onto a Clenshaw-Curtis (CC) grid. The GPU implementation requires the sample size to be divisible by 4444.

For iso-latitude rings, in order to achieve best efficiency for the Legendre transform part (the θ𝜃\thetaitalic_θ part of the transform, which is the critical part), modern top-performing CPU and GPU codes like SHTns use on-the-fly calculation of Pm(θ)superscriptsubscript𝑃𝑚𝜃P_{\ell}^{m}(\theta)italic_P start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT ( italic_θ ) using efficient recurrence formulas put forward recently [42]. This allows to keep memory usage low: indeed, only the recurrence coefficients that are independent of θ𝜃\thetaitalic_θ need to be stored, which requires only 𝒪(max2)𝒪superscriptsubscriptmax2\mathcal{O}(\ell_{\rm max}^{2})caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) memory, the same order as the data. It leaves two dimensions along which to parallelize: θ𝜃\thetaitalic_θ and m𝑚mitalic_m, and requires a sequential loop over \ellroman_ℓ to compute the Pm(θ)superscriptsubscript𝑃𝑚𝜃P_{\ell}^{m}(\theta)italic_P start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT ( italic_θ ) recursively. When maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT is larger than 500 to 1000, this leaves enough parallelization opportunities to efficiently use all of the GPU compute units. The computational complexity stays 𝒪(max3)𝒪superscriptsubscriptmax3\mathcal{O}(\ell_{\rm max}^{3})caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ).

The (adjoint) doubling (𝐃superscript𝐃\mathbf{D}^{\dagger}bold_D start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT) 𝐃𝐃\mathbf{D}bold_D is implemented via CUDA, and we write the arrays in a θ𝜃\thetaitalic_θ-contiguous memory layout, as required by SHTns to keep high efficiency for the Legendre transform. The computational complexity is 𝒪(max2)𝒪superscriptsubscriptmax2\mathcal{O}(\ell_{\rm max}^{2})caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ).

For the type-1111 and type-2222 nuFFT in 2222-dimensions (𝐍superscript𝐍\mathbf{N}^{\dagger}bold_N start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT, 𝐍𝐍\mathbf{N}bold_N), we use cufinuFFT [43] in double precision. The nuFFT method works by utilizing the Fourier transform convolution theorem, and interpolation or convolution onto a slightly larger, up-sampled grid. Highly accurate versions use kernels whose error ϵitalic-ϵ\epsilonitalic_ϵ decrease exponentially as a function of the up-sampling factor. The computational complexity (without planning phase) is 𝒪(max2log(max)+max2|log2(ϵ)|)𝒪superscriptsubscriptmax2subscriptmaxsuperscriptsubscriptmax2superscript2italic-ϵ\mathcal{O}(\ell_{\rm max}^{2}\log(\ell_{\rm max})+\ell_{\rm max}^{2}|\log^{2}% (\epsilon)|)caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_log ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ) + roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | roman_log start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ϵ ) | ) in 2 dimensions. It is worth mentioning that we use the guru interface to cufinuFFT to initialize the nuFFT plans. The plans allow for repeated and fast transforms, without re-initialization. However, this planning step is the most time consuming operation and is therefore done before calling the functions.

For the FFT operations 𝐅𝐅\mathbf{F}bold_F (type 2) and 𝐅superscript𝐅\mathbf{F}^{\dagger}bold_F start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT (type 1), we use the package cupyx [44] and its cuFFT555https://developer.nvidia.com/cufft integration therein. For the type 1 Fourier synthesis, we use double precision accuracy. For the type 2 Fourier synthesis, we use single precision accuracy for ϵ106italic-ϵsuperscript106\epsilon\leq 10^{-6}italic_ϵ ≤ 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT, which increases speed at the cost of a negligible decrease in effective accuracy on the final result. The computational complexity is 𝒪(max2log(max))𝒪superscriptsubscriptmax2subscriptmax\mathcal{O}(\ell_{\rm max}^{2}\log(\ell_{\rm max}))caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_log ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ) ).

FFTs become particularly fast if the prime factorization for the sample size gives many small prime numbers, in the following referred to as a good number. If additional constraints are put on the sample size, good numbers may be more difficult to find, see Appendix B for a discussion and concrete definition.

Ignoring the scaling with accuracy, the overall asymptotic computational complexity (for both type 1 and 2) is,

𝒪(max3)+𝒪(max2log(max))+𝒪(max2).𝒪superscriptsubscriptmax3𝒪superscriptsubscriptmax2subscriptmax𝒪superscriptsubscriptmax2\displaystyle\begin{split}\mathcal{O}(\ell_{\rm max}^{3})+\mathcal{O}(\ell_{% \rm max}^{2}\log(\ell_{\rm max}))+\mathcal{O}(\ell_{\rm max}^{2})\end{split}\,.start_ROW start_CELL caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ) + caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_log ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ) ) + caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) end_CELL end_ROW . (7)

Here, we assume that the number of uniform and nonuniform points are about the same.

𝐘𝐘\mathbf{Y}bold_Y (and 𝐘superscript𝐘\mathbf{Y}^{\dagger}bold_Y start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT) could be further optimized: 𝐒𝐒\mathbf{S}bold_S and 𝐅𝐅\mathbf{F}bold_F both contain a Fourier transform in ϕitalic-ϕ\phiitalic_ϕ-direction, effectively cancelling out each other. Avoiding this reduces 𝐅𝐅\mathbf{F}bold_F to a 1111-dimensional Fourier transform and 𝐒𝐒\mathbf{S}bold_S to a Legendre transformation. This optimization is implemented in the CPU implementation in DUCC. We leave this optimization to a future study for the GPU implementation. We only expect a large speed up for the 𝐅𝐅\mathbf{F}bold_F operator in cunuSHT, which takes about 10 to 20% of the total runtime.

For CMB weak lensing applications, we also provide a pointing routine implemented via CUDA, see Appendix C.

IV Benchmark

We present the scaling and execution time of our implementation. We take as a use case an application in CMB weak lensing [23]. CMB weak lensing describes the deflection of primordial CMB photons by mass fluctuations along the line of sight as they travel through the Universe. These small deflections (their root-mean-square is of the order of a few arcminutes) are large enough to be detectable in CMB sky maps [45, 5, 46]. For some applications, it is necessary to simulate these deflections accurately and efficiently [47, 27, 48].

Refer to caption
Figure 2: Setup and simulated data, here for a problem size of max=3095subscriptmax3095\ell_{\rm max}=3095roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 3095, N2107𝑁2superscript107N\approx 2\cdot 10^{7}italic_N ≈ 2 ⋅ 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT. The top left (right) panel shows a typical CMB temperature map (a typical deflection field) in orthographic projection. The bottom left plot shows a 10×10101010\times 1010 × 10 degree detail view of the difference between the undeflected and deflected CMB in Cartesian projection. On the bottom right, we show a 0.5×0.50.50.50.5\times 0.50.5 × 0.5 degree detail view comparison between the uniform grid (blue squares), nonuniform grid (orange squares), and the relation between them as indicated by the black arrows.

Owing to the deflections, the CMB intensity field f~(θ,ϕ)~𝑓𝜃italic-ϕ\tilde{f}(\theta,\phi)over~ start_ARG italic_f end_ARG ( italic_θ , italic_ϕ ) observed at location (θ,ϕ)𝜃italic-ϕ(\theta,\phi)( italic_θ , italic_ϕ ) is the un-deflected field f𝑓fitalic_f at another location,

f~(θ,ϕ)=f(θ,ϕ),~𝑓𝜃italic-ϕ𝑓superscript𝜃superscriptitalic-ϕ\tilde{f}(\mathbf{\theta},\mathbf{\phi})=f(\mathbf{\theta^{\prime}},\mathbf{% \phi^{\prime}})\,,over~ start_ARG italic_f end_ARG ( italic_θ , italic_ϕ ) = italic_f ( italic_θ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT , italic_ϕ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) , (8)

where (θ,ϕ)superscript𝜃superscriptitalic-ϕ(\theta^{\prime},\phi^{\prime})( italic_θ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT , italic_ϕ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) depends on (θ,ϕ)𝜃italic-ϕ(\theta,\phi)( italic_θ , italic_ϕ ) in a smooth way. Details on how these angles relate to each other are given in Appendix C. For the type-2 nuSHT, the inputs are the set of angles (θ,ϕ)superscript𝜃superscriptitalic-ϕ(\theta^{\prime},\phi^{\prime})( italic_θ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT , italic_ϕ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) and the spherical harmonic coefficients of the un-deflected field f𝑓fitalic_f. For type-1, the input are the same set of angles and a set of values of f~~𝑓\tilde{f}over~ start_ARG italic_f end_ARG. The set up of our use case is shown in Fig. 2. The top left panel shows the input CMB map, the deflection field is shown on the top right. Both are shown in orthographic projection. The bottom maps show a detail view of 10×10101010\times 1010 × 10 degree of the difference between the input and deflected map (left panel) and a detail view of approximately 0.5×0.50.50.50.5\times 0.50.5 × 0.5 degree (right panel) of the nonuniform point (orange squares) relative to the uniform point locations (blue squares). We use the same setup for the adjoint operation, in which case the deflected map becomes the input, and the result becomes the adjoint SHT coefficients.

Benchmarks are run on an NVIDIA A-100 GPU with 80808080 GB of memory and an Intel Xeon Gold 8358 Processor with 32 cores. We set the number of threads to 32323232 for the CPU benchmarks. If not stated otherwise, we choose the following parameters (that mostly affect the nuFFT): an up-sampling factor of 1.251.251.251.25 to reduce the memory usage at a small price of increased computation time666This somewhat low up-sampling factor reduces the minimal possible accuracy, in double precision, to 1010superscript101010^{-10}10 start_POSTSUPERSCRIPT - 10 end_POSTSUPERSCRIPT due to its dependence on the size of the up-sampled grid, but can easily be changed by the user, if needed., a gpu_method utilizing the hybrid scheme, called shared memory, and the default kernel evaluation method. The resulting map (or input map in the case of the adjoint) is calculated on a Gauss-Legendre grid.

With this implementation, we can solve problem sizes of up to max9000similar-tosubscriptmax9000\ell_{\rm max}\sim 9000roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ∼ 9000 on an A-100 with 80808080 GB, by using the pre-computed nuFFT plans and keeping all necessary intermediate results in memory.

Refer to caption
Refer to caption
Figure 3: cunuSHT execution time comparison against DUCC for type 2 (type 1) nuSHT shown in the left (right) panel, comparing one A-100 against one Intel Xeon Gold 8358 Processor with 32 cores. A speed up factor >1absent1>1> 1 means that the GPU is faster, and we show the result for different target accuracies. The error bars show the ±1σplus-or-minus1𝜎\pm 1\sigma± 1 italic_σ variance calculated from 5555 runs. The GPU takes over at max3102similar-tosubscriptmax3superscript102\ell_{\rm max}\sim 3\cdot 10^{2}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ∼ 3 ⋅ 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT for type 2, and between max4102similar-tosubscriptmax4superscript102\ell_{\rm max}\sim 4\cdot 10^{2}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ∼ 4 ⋅ 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT for type 1. The speed up increases up to 5555 (3) for large maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT left (right) panel.

Fig. 3 shows the speed up of the GPU algorithm compared to the CPU as a function of maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT for different accuracies. The left panel shows the evaluation of Eq. (5), the right panel shows Eq. (6). The ±1σplus-or-minus1𝜎\pm 1\sigma± 1 italic_σ variance from 5555 runs is indicated by error bars.

For type 2 nuSHT, we reach a speed up between 1 and 5 times for single precision and 1 to 3 times for double precision, with the speed up increasing with increasing maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT. This increase is expected due to better parallelizability for GPUs for higher maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT. The speed up is mostly independent of the target accuracy, with the smallest accuracies tending to perform better. For type 1 nuSHT, we find lower speed up factors, which may indicate further improvements. However, type 1 directions are generally expected to perform worse on the GPU due to the shear amount of threads that have to be written concurrently to the same memory location. Nevertheless, the GPU code either performs almost as good as the CPU code (small maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT), or better by a factor of up to 3333 for large problem sizes. For high maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT, the speed up depends on the accuracy, with lower accuracies performing better. For double precision (ϵ=1010italic-ϵsuperscript1010\epsilon=10^{-10}italic_ϵ = 10 start_POSTSUPERSCRIPT - 10 end_POSTSUPERSCRIPT), the speed up is diminished due to the double precision penalty that we pay in our implementation.

Refer to caption
Refer to caption
Figure 4: Breakdown and total execution time of the GPU implementation cunuSHT. The top panels show the percentage of time spent with the individual operators for type 2 (left column) and type 1 (right column) nuSHT. For each of the problem sizes we show this for different accuracies (1010superscript101010^{-10}10 start_POSTSUPERSCRIPT - 10 end_POSTSUPERSCRIPT on the left, 106superscript10610^{-6}10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT in the center, 102superscript10210^{-2}10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT on the right). The bottom panels show the total execution time, with ±1σplus-or-minus1𝜎\pm 1\sigma± 1 italic_σ variance as shaded area, for different accuracies and the empirically fitted computational complexity model, Eq. (D1).
Refer to caption
Refer to caption
Figure 5: Same as Fig. 4, but for the CPU implementation DUCC. It is important to note that DUCC implements the individual operators more efficiently by avoiding redundant Fourier transforms and by effectively combining the doubling and Fourier transform operations into one. To reflect this, we have grouped the 𝐃𝐃\mathbf{D}bold_D and 𝐅𝐅\mathbf{F}bold_F contributions in the top panel.

We can get a better understanding of the resulting speed up factors by looking at the time spent on each of the operators on the GPU. This breakdown is shown in the top panels of Fig. 4 for the type 2 (left panel) and type 1 (right panel) nuSHT, as a function of maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT for different accuracies. The respective total execution times are shown in the bottom panels together with the empirically fitted computational complexity model, Eq. (D1). At the top panels, for each problem size maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT, each bar represents a benchmark with an accuracy of (from left to right) 1010superscript101010^{-10}10 start_POSTSUPERSCRIPT - 10 end_POSTSUPERSCRIPT, 106superscript10610^{-6}10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT, 102superscript10210^{-2}10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT. Each bar represent the mean over 5 runs. For small problem sizes, doubling dominates and becomes almost negligible for large maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT. 𝐒𝐒\mathbf{S}bold_S only takes about 20% of the execution time for large maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT, even though it has the worst asymptotic computational complexity. This highlights the quality of the rSHT implementation by SHTns.

The choice of accuracy has an impact on the total execution time, as seen in the bottom panels of Fig. 4. The highest accuracy (ϵ=1010italic-ϵsuperscript1010\epsilon=10^{-10}italic_ϵ = 10 start_POSTSUPERSCRIPT - 10 end_POSTSUPERSCRIPT, red line) takes at most twice as long compared to the low accuracy (ϵ=102italic-ϵsuperscript102\epsilon=10^{-2}italic_ϵ = 10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT, purple line). Our results suggests that improvements in 𝐅𝐅\mathbf{F}bold_F might be possible: due to FFT being in principle a bandwidth limited routine, we would expect the execution time of 𝐅𝐅\mathbf{F}bold_F to be comparable to that of 𝐃𝐃\mathbf{D}bold_D, and only a fraction of that of 𝐒𝐒\mathbf{S}bold_S. For type 1 (right column), the total execution time overall takes longer. The breakdown shows that for the low and intermediate accuracy cases, less time is spend in the 𝐅superscript𝐅\mathbf{F^{\dagger}}bold_F start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT call. This is a consequence of our use of single precision arithmetic FFTs for the ϵ106italic-ϵsuperscript106\epsilon\leq 10^{-6}italic_ϵ ≤ 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT case, which, as mentioned before, is possible for the type 1 nuSHT.

It is interesting to compare this breakdown to the CPU implementation of DUCC. While we cannot expect both breakdowns to be exactly the same due to the different nature of the hardware, large differences may imply possible improvements. Fig. 5 shows the computation time of the CPU implementation of DUCC for the type 2 nuSHT (left column) and type 1 nuSHT (right column) as a function of maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT for different accuracies. Looking at the top left panel, we see that 𝐒𝐒\mathbf{S}bold_S increases with increasing maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT, as expected from its 𝒪(max3)𝒪superscriptsubscriptmax3\mathcal{O}(\ell_{\rm max}^{3})caligraphic_O ( roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ) scaling. The bottom panel of the left figure shows that the total execution time of the high-accuracy run is at most twice as long as the low accuracy one.

Many optimizations have gone into DUCC’s and SHTns’s SHT routines (𝐒𝐒\mathbf{S}bold_S, 𝐒superscript𝐒\mathbf{S}^{\dagger}bold_S start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT) over the years and it is safe to assume that they are close to optimal. Comparing the breakdown of the CPU and GPU implementation for 𝐒𝐒\mathbf{S}bold_S (green bars) shows that DUCC spends more than twice as long with this operator. This hints to potential sub-optimalities in the implementations of the GPU operators. This also becomes apparent when we look at the operator 𝐍𝐍\mathbf{N}bold_N (and 𝐍superscript𝐍\mathbf{N}^{\dagger}bold_N start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT). The scaling with the problem size is much more pronounced for the GPU code.

Fig. 6 shows the effective accuracy as a function of target accuracy for both CPU and GPU. The effective accuracy ϵeffsubscriptitalic-ϵeff\epsilon_{\rm eff}italic_ϵ start_POSTSUBSCRIPT roman_eff end_POSTSUBSCRIPT is calculated by solving Eq. (3) in a brute-force manner (giving us the values ftruesuperscript𝑓truef^{\rm true}italic_f start_POSTSUPERSCRIPT roman_true end_POSTSUPERSCRIPT) and comparing it against Eq. (5) (giving us the values f^estsuperscript^𝑓est\hat{f}^{\rm est}over^ start_ARG italic_f end_ARG start_POSTSUPERSCRIPT roman_est end_POSTSUPERSCRIPT),

ϵeff=1f¯truei=1Nt(fitruef^iest)2,subscriptitalic-ϵeff1superscript¯𝑓truesuperscriptsubscript𝑖1subscript𝑁𝑡superscriptsuperscriptsubscript𝑓𝑖truesuperscriptsubscript^𝑓𝑖est2\epsilon_{\rm eff}=\frac{1}{\bar{f}^{\rm true}}\sqrt{\sum_{i=1}^{N_{t}}\left(f% _{i}^{\text{true}}-\hat{f}_{i}^{\rm est}\right)^{2}}\,,italic_ϵ start_POSTSUBSCRIPT roman_eff end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG over¯ start_ARG italic_f end_ARG start_POSTSUPERSCRIPT roman_true end_POSTSUPERSCRIPT end_ARG square-root start_ARG ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_POSTSUPERSCRIPT ( italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT true end_POSTSUPERSCRIPT - over^ start_ARG italic_f end_ARG start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_est end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , (9)

with

f¯true=i=1Nt(fitrue)2.superscript¯𝑓truesuperscriptsubscript𝑖1subscript𝑁𝑡superscriptsuperscriptsubscript𝑓𝑖true2\bar{f}^{\rm true}=\sqrt{\sum_{i=1}^{N_{t}}\left(f_{i}^{\text{true}}\right)^{2% }}.over¯ start_ARG italic_f end_ARG start_POSTSUPERSCRIPT roman_true end_POSTSUPERSCRIPT = square-root start_ARG ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_POSTSUPERSCRIPT ( italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT true end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG . (10)

We choose Nt=106subscript𝑁𝑡superscript106N_{t}=10^{6}italic_N start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT random points on the sphere, giving us good sampling across the full sphere, and a sufficiently low variance on ϵeffsubscriptitalic-ϵeff\epsilon_{\rm eff}italic_ϵ start_POSTSUBSCRIPT roman_eff end_POSTSUBSCRIPT.

Refer to caption
Figure 6: Effective accuracy as a function of target accuracy for the CPU (blue), and GPU (orange) implementation.

Both algorithms achieve good effective accuracies, with the CPU code being more conservative. For the GPU implementation, we note that only low effective accuracies, ϵeff>104subscriptitalic-ϵeffsuperscript104\epsilon_{\rm eff}>10^{-4}italic_ϵ start_POSTSUBSCRIPT roman_eff end_POSTSUBSCRIPT > 10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT, can be achieved when nuFFT is executed in single precision. When an accuracy of ϵ=106italic-ϵsuperscript106\epsilon=10^{-6}italic_ϵ = 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT is desired, a single precision nuFFT evaluation is thus insufficient. We find that executing 𝐍𝐍\mathbf{N}bold_N and 𝐍superscript𝐍\mathbf{N}^{\dagger}bold_N start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT in double precision solves this. There is an associated penalty in efficiency due to the double precision calculations being approximately twice as slow. The nuFFT implementation in DUCC circumvents this by allowing for double precision accuracy on the pointing and intermediate results while using single precision accuracy on the Fourier coefficients, hereby reducing the execution time for single precision accuracies.

V Conclusion

We presented cunuSHT, a GPU accelerated implementation of the spherical harmonic transform on arbitrary pixelization that is, to our knowledge, the first of its kind to achieve faster execution when compared against CPU based algorithms. cunuSHT achieves machine precision accuracy by transforming the problem of interpolating on the sphere into a problem of computing a nonuniform fast Fourier transform on the torus. Comparing our implementation executed on an A-100 to the fastest available CPU implementations to date running on a single Intel Xeon Gold 8358 Processor with 32323232 cores, we find that our implementation is up to 5555 times faster. We used highly efficient, publicly available packages that are well tested and robust: SHTns for rSHTs, cufinuFFT for nuFFTs. We found that, although it has the highest asymptotic complexity, the high-quality rSHT implemented in SHTns is not the bottleneck. Our code does not require intermediate transfers between host and device, allowing it to be incorporated within larger GPU-based algorithms. Many applications in cosmology that we have in mind typically require spin-1111 to 3333 transforms. We thus plan to extend this package to spin-n𝑛nitalic_n transforms in the near future. There are, in principle, no obstacles to the generalization to spin-n𝑛nitalic_n by implementing the corresponding Wigner-d𝑑ditalic_d transform.

cunuSHT is a general purpose package distributed via pypi, and also works on standard pixelization schemes such as HEALPix, and can also perform rSHTs. Demonstrations of our package on GitHub present how to integrate it into exisiting pipelines.

Acknowledgment

The authors thank Alex Barnett for helpful discussions about nuFFT, Lehman Garrison for help on the cupy-nanobind implementation, and Libin Lu for general discussion and upgrades to cufinuFFT. This work is supported by the Research Analyst grant from the Simons Foundation and the computing resources of the Flatiron Institute. The Flatiron Institute is supported by the Simons Foundation. SB and JC acknowledges support from a SNSF Eccellenza Professorial Fellowship (No. 186879).

References

Appendix A Data Transfer

All benchmarks are done without measuring the time to transfer the data to device and back. This can be a large part of the overall computation and should be avoided, as is shown in Fig. A1 for different problem sizes. We note here that cunuSHT provides means to keep everything on the GPU without having to transfer the data. Thus, these transfer times can in principle be avoided.

Refer to caption
Figure A1: type 2 nuSHT transfer times “Host-to-Device” (H2D) (green) and “Device-to-Host” (D2H) (yellow) as a function of maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT. The shaded areas show the ±1σplus-or-minus1𝜎\pm 1\sigma± 1 italic_σ variance in transfer time, calculated from ten runs. Note that these transfer times are not generally part of the routine and can completely be avoided. H2D transfer contains the SHT coefficients, while D2H contains the much larger map. For type 1 nuSHT, the transfer times are equally long, but the H2D and D2H are swapped due to the input and output of the function.

Appendix B Good numbers

For fast Fourier transforms, execution time depends strongly on the largest prime factor of the transform length; the smaller it is, the better. Luckily enough, many good numbers for CPU algorithms exist. Assuming Clenshaw-Curtis quadrature, the constraint on the number of rings Nrsubscript𝑁𝑟N_{r}italic_N start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT to sufficiently sample a map with band limit maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT is Nrmax+1subscript𝑁𝑟subscriptmax1N_{r}\geq\ell_{\rm max}+1italic_N start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ≥ roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT + 1. Therefore, after applying the double Fourier sphere method, we must have at least 2Nr+22subscript𝑁𝑟22N_{r}+22 italic_N start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT + 2 samples, and we are free to choose an arbitrarily larger number which happens to be a good FFT size. For the GPU algorithm, however, there is a caveat: due to the additional constraint on Nrsubscript𝑁𝑟N_{r}italic_N start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT to be divisible by 4444 for the rSHT operation, the selection of good numbers is reduced. While there are still good numbers for FFT satisfying both constraints, they are sparser. The full list of Nrsubscript𝑁𝑟N_{r}italic_N start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT that are multiples of 4 up to Nr10000subscript𝑁𝑟10000N_{r}\leq 10000italic_N start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ≤ 10000 for which max=Nr1subscriptmaxsubscript𝑁𝑟1\ell_{\rm max}=N_{r}-1roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = italic_N start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT - 1 can be factored into primes up to 11 is 4, 8, 12, 16, 28, 36, 56, 64, 76, 100, 136, 148, 176, 232, 244, 276, 316, 344, 364, 376, 496, 540, 568, 676, 736, 848, 876, 892, 1156, 1216, 1324, 1332, 1376, 1576, 1716, 1816, 1876, 2080, 2188, 2476, 2696, 2836, 3088, 3268, 3376, 3676, 4236, 4376, 4456, 4852, 5104, 5776, 6076, 6616, 6656, 6876, 7204, 7624, 7876, 8020, 8576, 9076, 9376. The numbers with the highest prime factor being 11 are highlighted with an underline as efficient calculation of cunuSHT only applies to the CPU code here.

One way to enlarge this list could be for example by using a Fejér grid as intermediate grid, as opposed to the currently used Clenshaw-Curtis grid. This would effectively reduce the number of required samples to 2Nr2subscript𝑁𝑟2N_{r}2 italic_N start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT after doubling, resulting in a much wider choice of sample sizes that are also good FFT lengths. Hence restrictions on the number of good numbers are expected to go away in the future.

Appendix C CMB weak lensing pointing

Let e^θ,e^ϕsubscript^𝑒𝜃subscript^𝑒italic-ϕ\hat{e}_{\theta},\hat{e}_{\phi}over^ start_ARG italic_e end_ARG start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT , over^ start_ARG italic_e end_ARG start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT and n^^𝑛\hat{n}over^ start_ARG italic_n end_ARG form the right-handed unit basis vectors at the point on the sphere parametrized by θ,ϕ𝜃italic-ϕ\theta,\phiitalic_θ , italic_ϕ (so, the components of n^^𝑛\hat{n}over^ start_ARG italic_n end_ARG are sinθcosϕ,sinθsinϕ,cosθ𝜃italic-ϕ𝜃italic-ϕ𝜃\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\thetaroman_sin italic_θ roman_cos italic_ϕ , roman_sin italic_θ roman_sin italic_ϕ , roman_cos italic_θ). In our benchmark CMB lensing application, the deflected positions n^superscript^𝑛\hat{n}^{\prime}over^ start_ARG italic_n end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT that define the angles θsuperscript𝜃\theta^{\prime}italic_θ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT and ϕsuperscriptitalic-ϕ\phi^{\prime}italic_ϕ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT at which the CMB field must be evaluated are given by,

n^=cos(α)n^+sin(α)α(αθe^θ+αϕe^ϕ),superscript^𝑛𝛼^𝑛𝛼𝛼subscript𝛼𝜃subscript^𝑒𝜃subscript𝛼italic-ϕsubscript^𝑒italic-ϕ\hat{n}^{\prime}=\cos(\alpha)\hat{n}+\frac{\sin(\alpha)}{\alpha}\left(\alpha_{% \theta}\hat{e}_{\theta}+\alpha_{\phi}\hat{e}_{\phi}\right)\,,over^ start_ARG italic_n end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = roman_cos ( italic_α ) over^ start_ARG italic_n end_ARG + divide start_ARG roman_sin ( italic_α ) end_ARG start_ARG italic_α end_ARG ( italic_α start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT over^ start_ARG italic_e end_ARG start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT + italic_α start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT over^ start_ARG italic_e end_ARG start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT ) , (C1)

where α=αθ2+αϕ2𝛼superscriptsubscript𝛼𝜃2superscriptsubscript𝛼italic-ϕ2\alpha=\sqrt{\alpha_{\theta}^{2}+\alpha_{\phi}^{2}}italic_α = square-root start_ARG italic_α start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_α start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG, and αθsubscript𝛼𝜃\alpha_{\theta}italic_α start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT and αϕsubscript𝛼italic-ϕ\alpha_{\phi}italic_α start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT form the gradient of the lensing potential ΦΦ\Phiroman_Φ,

αθ(n^)+iαϕ(n^)=(θ+isinθϕ)Φ(n^).subscript𝛼𝜃^𝑛isubscript𝛼italic-ϕ^𝑛𝜃i𝜃italic-ϕΦ^𝑛\alpha_{\theta}(\hat{n})+\mathrm{i}\alpha_{\phi}(\hat{n})=\left(\frac{\partial% }{\partial\theta}+\frac{\mathrm{i}}{\sin\theta}\frac{\partial}{\partial\phi}% \right)\Phi(\hat{n}).italic_α start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG ) + roman_i italic_α start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG ) = ( divide start_ARG ∂ end_ARG start_ARG ∂ italic_θ end_ARG + divide start_ARG roman_i end_ARG start_ARG roman_sin italic_θ end_ARG divide start_ARG ∂ end_ARG start_ARG ∂ italic_ϕ end_ARG ) roman_Φ ( over^ start_ARG italic_n end_ARG ) . (C2)

cunuSHT first obtains αθsubscript𝛼𝜃\alpha_{\theta}italic_α start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT and αϕsubscript𝛼italic-ϕ\alpha_{\phi}italic_α start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT from this equation using SHTns. Then, threading across each ring, we solve for θ,ϕsuperscript𝜃superscriptitalic-ϕ\theta^{\prime},\phi^{\prime}italic_θ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT , italic_ϕ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT in Eq. (C1) on the fly. A benchmark as a function of maxsubscriptmax\ell_{\rm max}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT is shown in Fig. C1.

Refer to caption
Figure C1: Execution time of the pointing routine on the GPU as a function of the problem size. The shaded area shows the ±1σplus-or-minus1𝜎\pm 1\sigma± 1 italic_σ uncertainty calculated from 5 runs.

Appendix D Computational complexity

Table D1 shows the results for the empirical fit of the type 2 and type 1 nuSHT computational complexity models,

C()=α3norm3+β2log()norm2log(norm),𝐶𝛼superscript3superscriptsubscriptnorm3𝛽superscript2superscriptsubscriptnorm2subscriptnormC(\ell)=\alpha\frac{\ell^{3}}{\ell_{\rm norm}^{3}}+\beta\frac{\ell^{2}\log(% \ell)}{\ell_{\rm norm}^{2}\log(\ell_{\rm norm})}\,,italic_C ( roman_ℓ ) = italic_α divide start_ARG roman_ℓ start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG start_ARG roman_ℓ start_POSTSUBSCRIPT roman_norm end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG + italic_β divide start_ARG roman_ℓ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_log ( roman_ℓ ) end_ARG start_ARG roman_ℓ start_POSTSUBSCRIPT roman_norm end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_log ( roman_ℓ start_POSTSUBSCRIPT roman_norm end_POSTSUBSCRIPT ) end_ARG , (D1)

with norm=6067subscriptnorm6067\ell_{\rm norm}=6067roman_ℓ start_POSTSUBSCRIPT roman_norm end_POSTSUBSCRIPT = 6067 the normalization for the unknown prefactors. Compared to Eq. (7) we only account for the terms with the largest complexities. The fits are shown in Fig. 3, Fig. 4, and Fig. 5.

Table D1: Empirical fits of Eq. (7) for type 2 and type 1 nuSHT, for different accuracies for both CPU and GPU.
type 2 Target acc.        α𝛼\alphaitalic_α        β𝛽\betaitalic_β
CPU 1010superscript101010^{-10}10 start_POSTSUPERSCRIPT - 10 end_POSTSUPERSCRIPT 1.521011.52superscript1011.52\cdot 10^{-1}1.52 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT 6.961016.96superscript1016.96\cdot 10^{-1}6.96 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
106superscript10610^{-6}10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT 1.411011.41superscript1011.41\cdot 10^{-1}1.41 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT 6.161016.16superscript1016.16\cdot 10^{-1}6.16 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
102superscript10210^{-2}10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT 1.671011.67superscript1011.67\cdot 10^{-1}1.67 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT 4.111014.11superscript1014.11\cdot 10^{-1}4.11 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
GPU 1010superscript101010^{-10}10 start_POSTSUPERSCRIPT - 10 end_POSTSUPERSCRIPT 1.041061.04superscript1061.04\cdot 10^{-6}1.04 ⋅ 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT 2.721012.72superscript1012.72\cdot 10^{-1}2.72 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
106superscript10610^{-6}10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT 1.061021.06superscript1021.06\cdot 10^{-2}1.06 ⋅ 10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT 1.601011.60superscript1011.60\cdot 10^{-1}1.60 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
102superscript10210^{-2}10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT 1.391041.39superscript1041.39\cdot 10^{-4}1.39 ⋅ 10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT 1.191011.19superscript1011.19\cdot 10^{-1}1.19 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
type 1 Target acc. α𝛼\alphaitalic_α β𝛽\betaitalic_β
CPU 1010superscript101010^{-10}10 start_POSTSUPERSCRIPT - 10 end_POSTSUPERSCRIPT 1.241011.24superscript1011.24\cdot 10^{-1}1.24 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT 7.341017.34superscript1017.34\cdot 10^{-1}7.34 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
106superscript10610^{-6}10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT 1.061011.06superscript1011.06\cdot 10^{-1}1.06 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT 6.671016.67superscript1016.67\cdot 10^{-1}6.67 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
102superscript10210^{-2}10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT 1.681011.68superscript1011.68\cdot 10^{-1}1.68 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT 3.721013.72superscript1013.72\cdot 10^{-1}3.72 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
GPU 1010superscript101010^{-10}10 start_POSTSUPERSCRIPT - 10 end_POSTSUPERSCRIPT 1.311081.31superscript1081.31\cdot 10^{-8}1.31 ⋅ 10 start_POSTSUPERSCRIPT - 8 end_POSTSUPERSCRIPT 4.331014.33superscript1014.33\cdot 10^{-1}4.33 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
106superscript10610^{-6}10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT 1.071061.07superscript1061.07\cdot 10^{-6}1.07 ⋅ 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT 2.841012.84superscript1012.84\cdot 10^{-1}2.84 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
102superscript10210^{-2}10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT 2.161052.16superscript1052.16\cdot 10^{-5}2.16 ⋅ 10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT 1.641011.64superscript1011.64\cdot 10^{-1}1.64 ⋅ 10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT

Appendix E Code examples

The following code block shows a minimum working example for calculating a type 2 and type 1 nuSHT on the GPU for SHT coefficients alm that are deflected by a deflection field dlm_scaled, for an accuracy of epsilon and band limit lmax. Both routines will set up the plans, so that the actual nusht2dX() call may be called repeatedly. 1import cunusht as cu 2 3lenjob_geominfo = (’gl’,{’lmax’: lmax}) 4kwargs = { 5 ’geominfo_deflection’: lenjob_geominfo, 6 ’epsilon’: epsilon, 7 ’nuFFTtype’: 2, 8} 9t = cu.get_transformer(backend=’GPU’)(**kwargs) 10ptg = t.dlm2pointing(dlm_scaled) 11lenmap = t.nusht2d2(alm, ptg, lmax, lenmap) 1import cunusht as cu 2 3lenjob_geominfo = (’gl’,{’lmax’: lmax}) 4kwargs = { 5 ’geominfo_deflection’: lenjob_geominfo, 6 ’epsilon’: epsilon, 7 ’nuFFTtype’: 1, 8} 9t = cu.get_transformer(backend=’GPU’)(**kwargs) 10ptg = t.dlm2pointing(dlm_scaled) 11alm = t.nusht2d1(alm, ptg, lmax, lenmap)