Spherical Harmonic Transforms

The main target of HealpixMPI.jl, in terms of run time speed up, are the SHTs, which often represent the "bottle neck" of simulation codes or data analysis pipelines.

As mentioned in the introduction, HealpixMPI.jl relies on ducc's state-of-the-art algorithms for performing the spherical harmonic transforms. In particular, its C++ functions are exploited for the computation of Legandre coefficients from alms and maps and vice versa, while the transposition of such coefficients between MPI tasks is entirely coded in Julia.

HealpixMPI.jl implements only the two exact spherical harmonic operators alm2map! (synthesis) and adjoint_alm2map! (adjoint synthesis), leaving to the user the task to implement the two corresponding inverse approximate operators, which can be done through many different approaches, from pixel and ring weights to conjugate gradient solver. Both alm2map! and adjoint_alm2map! are implemented as overloads of the Healpix.jl's functions.

From Alm to Map: synthesis operator

The synthesis SHT (alm2map!) is used to compute a map from a set of $a_{\ell m}$ coefficients. It is generally represented by the matrix operator $\mathrm{Y}$ which is defined through an exact summation as $f(\theta, \phi) = \mathrm{Y} \, a_{\ell m} \quad \text{where} \quad f(\theta, \phi) = \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} a_{\ell m} Y_{\ell m} (\theta, \phi).$

Healpix.alm2map! — Function

alm2map!(d_alm::DAlm{S,N}, d_map::DMap{S,T}, aux_in_leg::StridedArray{Complex{T},3}, aux_out_leg::StridedArray{Complex{T},3}; nthreads::Integer = 0) where {S<:Strategy, N<:Number, T<:Real}
alm2map!(d_alm::DAlm{S,N}, d_map::DMap{S,T}; nthreads::Integer = 0) where {S<:Strategy, N<:Number, T<:Real}

This function performs an MPI-parallel spherical harmonic transform, computing a distributed map from a set of DAlm and places the results in the passed d_map object.

It must be called simultaneously on all the MPI tasks containing the subsets which form exactly the whole map and alm.

It is possible to pass two auxiliary arrays where the Legandre coefficients will be stored during the transform, this avoids allocating extra memory and improves efficiency.

Arguments:

d_alm::DAlm{S,N}: the MPI-distributed spherical harmonic coefficients to transform.
d_map::DMap{S,T}: the MPI-distributed map that will contain the result.

Optionals:

aux_in_leg::StridedArray{Complex{T},3}: (localnm, totnring, ncomp) auxiliary matrix for alm-side Legandre coefficients.
aux_out_leg::StridedArray{Complex{T},3}: (totnm, localnring, ncomp) auxiliary matrix for map-side Legandre coefficients.

Keywords

nthreads::Integer = 0: the number of threads to use for the computation if 0, use as many threads as there are hardware threads available on the system.

source

From Map to Alm: adjoint synthesis operator

The adjoint of the synthesis operator brings us from the map space to the harmonic space, as it is represented by the transpose $\mathrm{Y}^{\mathrm{T}}$. Which is defined through: $\mathrm{Y}^{\mathrm{T}} f(\theta, \phi) \equiv \sum_{i = 1}^{N_{\mathrm{pix}}} Y^*_{\ell m,\, i} \, f_i.$

Note that this does not give directly the $a_{\ell m}$ coefficients, i.e., $\mathrm{Y}^{\mathrm{T}} \mathrm{Y} \neq \mathbf{1}.$ In fact, $\mathrm{Y}^{-1} \simeq \mathrm{W}\, \mathrm{Y}^{\mathrm{T}}.$ Where $\mathrm{W}$ is a diagonal matrix whose non-zero elements are approximately constant and equal to $4 \pi / N_{\mathrm{pix}}$, depending on the map pixelization. $\mathrm{Y}^{-1}$ is in fact an integral operator which must be approximated when implemented numerically.

Healpix.adjoint_alm2map! — Function

adjoint_alm2map!(d_map::DMap{S,T}, d_alm::DAlm{S,N}, aux_in_leg::StridedArray{Complex{T},3}, aux_out_leg::StridedArray{Complex{T},3}; nthreads = 0) where {S<:Strategy, N<:Number, T<:Real}
adjoint_alm2map!(d_map::DMap{S,T}, d_alm::DAlm{S,N}; nthreads::Integer = 0) where {S<:Strategy, N<:Number, T<:Real}

This function performs an MPI-parallel spherical harmonic transform Yᵀ on the distributed map and places the results in the passed d_alm object.