Caroline Kulcsar

M John Romein

M Charles

Nicolas Monnier

Nicolas Gac

Cyril Tasse

Erwan Raffin

David Guibert

ExaSKA : Parallelization on a High Performance Computing server for the exascale radiotelescope SKA

Keywords: Parallélisation sur serveur Calcul Parallèle, Radio Astronomie, Problème Inverse, Interpolation, Exascale, GPU Parallel Computing, Radio Astronomy, Inverse Problem, Interpolation, Exascale, GPU . . . . . 78 Radioastronomy, Inverse Problem, Interpolation, GPU, Gridding Contents

La radioastronomie permet de faire des observations du ciel à des longueurs d'ondes invisible à l'oeil nu et ainsi obtenir de nombreuses informations sur l'univers qui nous entoure. Les réseaux d'antennes qui composent ces instruments d'observation génèrent un flux continu de données à des débits qui rendent leur traitement et leur stockage presque impossible. En plus du très grand nombre de données à traiter, la reconstruction d'images du ciel nécessite des algorithmes itératifs coûteux en temps et en puissance de calcul rendant les plus gros supercalculateurs insuffisant pour le faire en temps-réel. Cette thèse vise à explorer les différentes manières de réduire le temps de calcul pour la reconstruction d'une image en se concentrant à la fois sur l'aspect algorithmique et implémentation. Une méthode algorithmique permettant de réduire le coût calculatoire des opérateurs de gridding et degridding est proposée. Une implémentation sur plateforme accélératrice ainsi qu'à grande échelle sur plusieurs noeuds de calcul d'un supercalculateur est également présentée.

Abstract:

Radio astronomy allows us to make observations of the sky at wavelengths invisible to the naked eye and thus obtain a lot of information about the universe that surrounds us. The antenna arrays that make up these observation instruments generate a continuous data flow at rates that make their processing and storage almost impossible. In addition to the very large amount of data to process, the reconstruction of sky images requires iterative algorithms that are costly in terms of time and computational power, making the largest supercomputers insufficient to do it in real-time. This thesis explores different ways to reduce the computation time for image reconstruction by focusing on the algorithmic aspect and the implementation aspect. We have proposed an algorithmic method to reduce the computational cost of expensive operators. We have also proposed an implementation on an accelerator platform and a large-scale implementation on several computational nodes of a supercomputer. 2.4 A basic two-element interferometer observing an arbitrary celestial source. The signal from the source to the antennas is corrupted along its path. The signal is converted into independent voltage v p and v q after reaching two antennas p and q. These voltages are fed into a correlator, which produces four pairs or correlations defined by the visibility matrix V pq . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.5

The (lmn) sky plane is tangent to the celestial sphere at the phase center and parallel to the plane in which the measurement are made (the uv-plane). The w axis points towards the phase center (l = 0, m = 0) [Bhatnagar 2001

Mathematical symbols

• R : real numbers.

• C : complex numbers.

• a : scalar quantity.

• a : column vector.

• A : matrix.

• (.) * : complex conjugate operator.

• (.) T : transpose operator.

• (.) † : adjoint operator.

• (.) H : Hermitian transpose operator.

• ⊗ : Kronecker product between two matrices.

• * : convolution product.

• ∥a∥ 2 : l 2 norm.

• < • > : Average over small time and frequency.

Preamble

This thesis, realized in laboratoire des Signaux et Système at Centralesupelec, was co-financed by Région Île-de-France in its project of support of scientific employment in connection with the needs of the companies, and by ATOS.

Motivation of this work

The visible light emitted by celestial objects does not tell us everything about their current state, history, and future evolution. Planets, stars, and galaxies naturally emit electromagnetic waves in the radio frequency range. These "invisible" observations allow us to trace hydrogen, the most abundant element in the Universe, in various states and thus about the universe's formation or the structuration of galaxies due to hydrogen, the most abundant element. Radio astronomers, therefore, study the sky at many wavelengths to understand better what surrounds us in the universe. One means at their disposal is to make images, or hyperspectral images, from the radio information of all these celestial structures. The sky is observed using a radio interferometer, an antenna array that simulates a single telescope by exploiting interferometry techniques and aperture synthesis. Unfortunately, the observation of the sky is not straightforward, and many steps of data processing are mandatory. The signals received by the antennas are corrupted 2 Chapter 1: Introduction by many effects, from the emission of electromagnetic waves to reception. Moreover, the samples, named visibilities, correspond to Fourier coefficients of the observed sky image. Therefore, reconstructing the sky image is a problem of Fourier synthesis where the data are corrupted and must be corrected.

Modern radio telescopes such as LOw-Frequency ARray (LOFAR) or Square Kilometre Array (SKA) offer a large collection area, a very high sensitivity, and a high resolution. However, due to the unprecedented data rate generated by the antennas, reaching the optimal theoretical performance for image reconstruction is a real challenge. Indeed, reconstructing an image from irregular and sparse Fourier coefficients requires iterative algorithms by adding a priori information on the objects to reconstruct.

In this thesis, we propose specific methods and implementations to accelerate image reconstruction algorithms in order to be ready for the new radio telescopes.

Contributions

The ways to reduce the computation time for reconstructing an image with a given algorithm can be divided into two points. The first one is based on the algorithmic part. We can modify some parts without changing the algorithm to speed up the process while keeping the same final result. The second one is the software implementation on dedicated hardware. Modern hardware encourages parallelization because of its design in order to speed up the calculations. This parallelization can be done at fine grain by optimizing a specific computational kernel on a dedicated platform. It can also be done at a large scale by starting from the algorithm and splitting the work into several computational resources.

Our contributions will focus on these two points. First, we will propose a method for interpolation on a Fourier grid to decrease this operation's computational cost. Indeed, interpolation is used to map visibilities on a uniform grid (or vice versa) in order to use an FFT. The computational cost of this step when using iterative algorithms on a large number of data is quickly prohibitive. Our method, named Grid to Grid (G2G), makes the computational cost decrease keeping the same accuracy as standard methods. We will also propose a parallel implementation on multi-core CPU and GPU of this method.

Then, we will propose a multi-core multi-node parallelization of the DDFacet imaging framework, named McMn DDFacet, in order to reduce the execution time of a sky image reconstruction. The amount of data to be processed with modern radio telescopes and the complexity of the corrections of the effects corrupting the data make the computation time prohibitive. DDFacet is at the forefront of correcting direction-dependent effects for a wide field of view observations. Our method decreases the computation time and gives hints to find a balance between the hardware resources used and the speedup.

Thesis Outline 3

Thesis Outline

• Chapter 2 introduces the problem setup and the state-of-the-art. All notations and mathematical problems related to radio interferometry imaging are introduced. The different image reconstruction methods are also presented in more detail, including the most commonly used by the community. We also review the challenges related to imaging with recent radio telescopes to underline the interest in the work of this thesis.

• Chapter 3 presents algorithmic interpolation methods, gridding and degridding, used in image reconstruction algorithms. We propose a method named Grid to Grid (G2G), allowing us to reduce the computation time and the memory footprint. Moreover, we show the implementation of this method on CPU and GPU accelerator platform in order to compare it with state-of-the-art.

• Chapter 4 presents the distribution of a DDFacet imaging framework on an HPC server. We present this implementation in detail ranging from multicore parallelization with shared memory to multi-node parallelization with distributed memory.

Chapter 2

Problem setup and state-of-the-art Radio astronomy imaging aims to reconstruct sky images from large-scale data acquired by radio telescopes. The model of these data is not trivial because of the effects disturbing the measured signal. Consequently, sky reconstruction is complex and computationally expensive because these effects must be corrected for processing the huge amount of data. Recent radio telescopes, such as the SKA, have such a high data rate to process that the algorithms and infrastructure to process the data need to scale up. Indeed, with SKA, dataflows with sizes of several TB are directly sent in the pipeline for quasi-real-time processing without any storage capabilities. We will focus on the Science Data Processor stage, which is responsible for the reconstruction through iterative algorithms deployed on large 6 Chapter 2: Problem setup and state-of-the-art HPC servers. Therefore, working on the algorithm optimization and the HPC scalability deployment is required to meet the real-time requirement.

In this chapter, we will first see how a radio telescope works and how to reconstruct an image from the acquired data. We will then discuss the different algorithmic methods of sky reconstruction. Finally, we will see the computational challenges related to the large amount of data to process to reconstruct images as quickly as possible.

Context of study

This section is about the study context of radio astronomy. Radio telescopes generate a large amount of data, increasing with the most recent instruments. This data must be processed, which is very expensive because of the irregular sampling.

Radio Interferometer

Radio astronomy aims to observe radio-wave emissions from structures of the Universe to obtain more information on this wide variety of structures, such as galaxies, pulsars, or quasars. We have expected astronomical structures to emit electromagnetic radio waves since the 1860s and Maxwell's equation. In this chapter, we will first see how a radio telescope works and how to reconstruct an image from the acquired data. We will then discuss the different algorithmic methods of sky reconstruction. Finally, we will see the computational challenges related to the large amount of data to process to reconstruct images as quickly as possible. However, it was only during the 1930s that K.Jansky accidentally received a signal from the center of the Milky Way while working on a large directional antenna at the Bell Telephone Laboratories.

The real beginning of radio astronomy as a science allowing observation of celestial structures started after the Second World War, in the 1950s, when the cost of radar equipment had drastically dropped. At first, mainly single-antenna radio telescopes were used. For such a telescope, the angular resolution is given by

θ ∼ λ D , (2.1)
where λ is the wavelength, and D is the diameter of the instrument. In radio, the resolution is commonly insufficient. In order to obtain a better one, as well as a higher sensitivity, scientists have built antenna receivers with larger collecting areas.

To illustrate this work, Fig. 2.1 compares the 9 meters of Ruber's dish, built in 1937, and the famous Arecibo Observatory in Puerto Rico, built in 1963, with a diameter of 305 meters [Goldsmith 1996]. More recently, the radiotelescope Fast in China, fully operational since 2020, has diameter of 500 meters, reaching an resolution of 2.9 arcmins at the 21 cm wavelength [Nan 2006]. However, these radio telescopes have limits. The main one is the difficulty of building colossal structures that require a high financial cost as well as several workers to provide efficient maintenance. Moreover, the limit of building a large dish implies a limitation in improving the angular resolution. In order to overcome this problem, scientists had to find innovative techniques such as interferometry and aperture synthesis.

The first radio observation using two dipoles, before the concept of aperture synthesis theorized by Ryle in 1960 [START_REF] Ryle | [END_REF]], was made in 1946 by Ryle and Vanberg. The principle is to simulate a large antenna by collecting signals from several antennas instead of one. For a pair of antennas physically separated by a distance named baseline, each antenna's measured signal is correlated with the signal of the other antenna to measure the delay between these two signals. Since this simulated telescope can be composed of an array of several antennas, its virtual diameter, also called aperture, is the largest distance between two antennas in the array. The angular resolution of a radio interferometer is then given by

θ ∼ λ B max , (2.2)
where B max is the longest baseline of the instrument. Large radio telescope projects composed of several dishes are widespread because they offer the advantage of the high sensitivity brought by the dishes in addition to a better angular resolution. The Very Large Telescope (VLA), built in the 1980s in the USA and modernized many times since its construction, consists of 27 antennas of 25 meters in diameter, covering a frequency range from 74 MHz to 50 GHz [START_REF] Perley | [END_REF]]. The Atacama Large Millimeter/submillimeter Array (ALMA), a network of 66 dishes of 12 meters in diameter, is a radio telescope built in Chili in 2003 covering the frequency band from 31 GHz to 1000 GHz [Webber 2013]. Besides this lack of flexibility, radio interferometers built with an array of dishes are not very flexible because it is necessary to mechanically move the antennas to observe in the direction of interest in the sky. Moreover, the construction cost is still relatively high because of the large number of antennas to build.

The new generation of radio interferometers, also called Software Telescope (ST), is now composed of many small, fixed, and omnidirectional antennas. Unlike the dishes arrays, to change this new generation's observation direction, we apply a phase shift to the data at a software level, equivalent to simulating the observation's phase center in another direction. After being correlated, a beam is steered toward the observation and applied to the signal measured by the antennas. If we count n a antennas in our radio telescope, the number of baselines is

N bl = n a (n a -1) 2 (2.3)
The number of baselines is the instantaneous number of correlation measurements. As these radio telescopes usually have many more antennas than those with dishes, the effect on the final result is a better sensitivity, a better resolution, and thus a better image quality. Great Britain, Ireland, Latvia, and Sweden. This telescope provides a large field of view in the 10-240MHz frequency band. In addition, several observations can be made in parallel, with up to 488 beams being formed. Another example is the SKA telescope, which is still under construction. It is the biggest radio telescope project ever made. This radio telescope will be built on two sites. The first one, called SKA-mid, will be mainly in South Africa, with several stations in other countries of Africa. This station will have 197 dishes to cover the 350 MHz to 15 GHz frequency range. The second site, called SKA-low, will be in Australia and will have more than 130.000 antennas gathered in 256 stations, similar to LOFAR, to cover the 50 MHz to 350 MHz frequency range. The objectives of this telescope are very ambitious because of the high sensitivity as well as the high resolution that the instrument will provide. The main science cases will be the observation of the epoch of reionization, some tests of general relativity through the study of pulsars, the search for extraterrestrial life, and many others.

Radio Interferometric Measurement Equation

Let us consider a sky consisting of a single point source which can be described as a single source of quasi-monochromatic signal. At a specific time, the signal lives in a 3-dimensional orthonormal coordinate system (x, y, z) space where z stands for the direction of propagation of the signal. The amplitude of the signal e can be described by a column vector of two complex numbers e x and e y :

e = e x e y ∈ C 2 (2.4)
The main assumption is that all transformations along the signal path w.r.t. e are linear [Hamaker 1996, Smirnov 2011].

When the antennas receive the signal, the signal is converted into complex voltages by the antenna feeds. Assuming two feeds, a and b (two linear dipoles), the corresponding voltage v a and v b are then linear w.r.t e. Thus, we define this complex voltage for the antenna p as

v p (θ) = v pa v pb e = J p (θ)e (2.5)
where θ is an unknown parametrization vector and J p is known as the Jones matrix [START_REF] Jones | [END_REF]]. The Jones matrix is a 2 × 2 complex matrix that describes all the perturbations (considered linear) that occur along the signal path from the emitting source to the receiving antenna p.

The array's antennas, p and q, measure two independent voltages v p and v q . As described in Fig. 2.4, in an interferometer, these voltages are fed into a correlator to produce four pairs of correlations

< v pa v * qa >, < v pa v * qb >, < v pb v * qa >, < v pb v * qb >. v *
qb is the complex conjugate of v qb , and brackets denote averaging over small time and frequency bin. The 2 × 2 resulting correlation matrix is then defined by [Smirnov 2011] as

V pq (θ) = E{v p (θ)v H q (θ)} = J p (θ)XJ H q (θ) (2.6)
where

X = E{ee H } = I + Q U + iV U -iV I -Q is the brightness matrix.
In this notation, (I, Q, U, V ) are the stokes parameters describing the polarization state of the signal [START_REF] Stokes | On the Composition and Resolution of Streams of Polarized Light from different Sources[END_REF]]. These parameters respectively represent the total intensity, two intensities of linear polarization, and the intensity of circular polarization. As the Jones matrices, in Eq.2.6, are considered constant over the integration time interval, they can be taken out of the expectation matrix. Eq 2.6 is the first form of the Radio Interferometric Measurement Equation (RIME) for a single source. Jones matrices J p and J q describe all the perturbations that corrupt the original signal. With the assumption of linearity of all transformations along the path w.r.t. e, J p and J q can be decomposed into the multiplication of several 2 × 2 matrix where each one is describing a specific perturbation. Thus, we defined the multiple matrix multiplication as the Jones chain, such as

J p = J p,n ...J p,2 J p,1 .
(2.7)

The order of each element of the Jones chain describes the physical order of signal perturbation. For example, J p,1 is the last perturbation and corresponds to the electronic gain of the antenna. Moreover, the order cannot be changed as matrix multiplication is non-commutative unless diagonal. These Jones matrices describe different kinds of corruption and can be represented as a scalar that will affect both signal components equally, a diagonal matrix that will affect each component differently, or a rotation matrix. A specific Jones chain has been built over the years in the literature [Hamaker 1996, Noordam 1996]. We will use a simplified version of this Jones chain, such as

J p = G p D p K p , (2.8)
with G p the diagonal matrix representing the electronic complex gain of the p th antenna, D p the matrix representing all direction-dependent effects like the perturbation of the ionosphere or the primary beam, and K p a scalar depending on the geometry of the antenna array. K p describes the phase difference with respect to the phase center due to the difference in the propagation path of the signal. Its expression is defined by

K p = exp -2iπ(u p l + v p m + w p (n -1)) ,
(2.9)

Figure 2.4: A basic two-element interferometer observing an arbitrary celestial source. The signal from the source to the antennas is corrupted along its path.

The signal is converted into independent voltage v p and v q after reaching two antennas p and q. These voltages are fed into a correlator, which produces four pairs or correlations defined by the visibility matrix V pq .

with (u p , v p , w p ) the position of the p th antenna in unit of wavelength such as the phase center is (u 0 = 0, v 0 = 0, w 0 = 0), and the coordinate (l, m, n = √ 1 -l 2 -m 2 ) is the direction-cosine of the arriving signal, see Fig. 2.5. Considering a signal from a single source and taking into account only the effects related to the geometry of our radiotelescope, Eq 2.6 becomes

V pq = K p XK H q = Xe -2iπ upql+vpqm+wpq(n-1) ,
(2.10) with u pq = u p -u q , and V pq expressing the visibility as a function of the baseline of the uvw coordinate system.

In real life, there is not a single source or multiple discrete sources but rather a continuous distribution of the sky brightness X(l, m), with each antenna having a specific Jones matrix J p (l, m) for each direction. Thus, the full sky RIME, taking Eq 2.10 into account becomes

V pq = G p lm D p (l, m)X(l, m)D H q (l, m)e -2iπ(ul +vm+w(n-1)) dl dm n G H q .
(2.11)

Figure 2.5: The (lmn) sky plane is tangent to the celestial sphere at the phase center and parallel to the plane in which the measurement are made (the uv-plane).

The w axis points towards the phase center (l = 0, m = 0) [Bhatnagar 2001].

The RIME formalism using Jones matrix is helpful for understanding and considering all the effects disturbing the signal easily. However, another formalism is widely used: the RIME using 4 × 4 Mueller matrices [Thompson 2001, Rau 2009]. We then define G pq = G p ⊗ G H q and J pq = D p ⊗ D H q as the direction-independent and direction-dependent 4 × 4 Mueller matrices for the baseline pq at a specific time and frenquency, and ⊗ the Kronecker product. Taking the vectorized version of x = vec(X) and v pq = vec(V pq ), the visibility equation becomes v pq = G pq lm J pq x(l, m)e -2iπ(ul +vm+w(n-1)) dl dm n .

(2.12)

From now, we'll consider the visibilities corrected for the electronic complex gain G pq , the interferometer array is coplanar, and we observe an infinitesimal source with a small Field of View (FoV). We can then apply Van Cittert-Zernike theorem, and Eq 2.12 becomes

v pq = lm x(l, m)e -2iπ(ul +vm) dl dm (2.13)
This simplification highlights the nature of the data. Therefore, the sampling of the Fourier transform of the intensity of the observed source provides the visibilities.

In that case, the problem of imaging is a problem of Fourier synthesis.

Synthesis imaging

The Earth's rotation causes the baselines to rotate as well, so the position of the visibilities changes over time. The total number of visibilities, M , depends on the number of baselines N bl , the number of frequencies N chan , the number of polarization N pol , the duration of the observation T , and the integration time for a sample ∆t, such as

M = N bl × N chan × N pol × T ∆t .
(2.14)

The uv-coordinates of the visibilities will therefore change over time. Fig. 2.6 shows the uv-coverage evolution regarding the VLA's observation time. v i e 2iπ(u i l+v i m) .

(2.15)

The dirty image reconstruction quality depends on the uv-coverage. The greater the coverage, the better the image quality; therefore, reconstructing the real sky will be easier. On the other hand, this also comes with a higher computational cost to make an image. Fig. 2.7 shows the evolution of the image quality obtained as a function of the uv-coverage.

The most straightforward approach to compute the dirty image is to use the Direct Fourier Transform method. In this case, it involves computing the 2D DFT of an N p × N p = P pixels image for the M visibilities of the observation. The computational complexity is well known to be

C DF T = O(P M ).
(2.16)

However, the size of the images (> 10e8 pixels) and the number of visibilities generated by recent radiotelescopes make this method prohibitive. Indeed, LOFAR has 50 stations, and it is planned to build 512 stations for SKA-Low. The number of baselines and, thus, the number of visibilities generated will increase exponentially, so this method will not be able to scale up.

Using the FFT would avoid this greedy DFT algorithm but requires that the visibilities' uv position are sampled uniformly, which is not the case. The most common solution is to interpolate visibilities with non-uniform uv-coverage into a uniform uv-coverage, an interpolation known as gridding. More technical details will be given in Sec. 3. The main idea is to convolve the data with an appropriate convolution kernel whose result is sampled on the uniform grid. The Fourier grid g at coordinate (u j , v j ) is obtained as

g(u j , v j ) = M i=0 C † (u i -u j , v i -v -j)v i ,
(2.17)

where C † is the gridding convolution kernel. Once all the visibilities are gridded on the grid g, an inverse FFT is applicable to obtain the dirty image y. In this case, the computational complexity becomes

C grid = O(P log 2 P + C 2 supp M ) (2.18)
where C 2 supp is the size of the convolution kernel in 2D.

Calibration

The incoming radio electromagnetic waves are affected by several effects along the path from the source to the antenna array. The perturbations effects, involving propagation delays, amplitude drops, or polarization disturbances, are characterized into two types. The first type is the instrumental effect. Indeed, each antenna has a unique electrical gain and an anisotropic beam pattern. The antennas are, therefore, sensitive to direction due to the sidelobes of the beam. The second is the environmental effect. Indeed, the troposphere and the ionosphere, two layers of the atmosphere, affect the incoming signal with diffraction and refraction at specific frequencies. The results of all these perturbations are artifacts of the reconstructed image.

The purpose of calibration is to estimate and correct these perturbations. The RIME allows modeling these perturbations in the form of Jones matrices. For example, G p is the Jones matrix to be determined, which characterizes the electrical gain of the antenna p. D p , which in our model corresponds to the Direction-Dependent Effects (DDEs), can be expanded to characterize each specific perturbation, such as ionospheric effects or antenna beam. The calibration methods depend on the architecture of the radio telescope. Several regimes are presented in the literature [Wijnholds 2010, Ollier 2018], and each regime depends on the size of the FoV and the distance between the antennas. Some regimes are illustrated in Fig. 2.8.

Beyond that, two main strategies exist for calibration. The first one is the external calibration strategy. This method estimates the unknown gain from sources whose structure, intensity, and position are well known. However, this method is limited to sources with a calibration close to the FoV. Moreover, when we have a large FoV, the estimations are correct only in a small area of the FoV. The second method is called self-calibration. In this method, both the sky and the perturbation parameters are estimated. The sky is estimated from perturbation parameters that are considered known. Then, the perturbations are estimated from the previously estimated sky. Both are then estimated iteratively as a form of alternating optimization.

Radio Interferometric Imaging

Radio interferometric imaging aims to reconstruct an estimate of the true x sky from the visibilities v. Variations of the CLEAN algorithm are among the most widely used methods in the radio astronomy community, although many other methods exist.

Inverse Problem

By stacking all the visibilities acquired at each time, we get a discrete version of the forward problem. Thus, the equation Eq 2.13 can be written

v = SF x + n (2.19)
where v ∈ C M is the visibility measurement as a vector of M samples, the true sky can be represented by N p × N p = P pixels grid as x ∈ R P . F and S are linear operators such as F ∈ C P ×P is the normalized Fourier transform, and S ∈ C M ×P the degridding operator maps P Fourier coefficients to M visibilities. In this model, we consider an additive independent and identically distributed Gaussian noise n ∈ C M such that n ∼ N (0, Σ). With the a priori information that the noise is not correlated, the covariance is a diagonal matrix

Σ =        σ 2 1 0 . . . 0 0 σ 2 2 . . . . . . . . . . . . 0 . . . . . . σ 2 M       
where σ i is the variance of the visibility v i . Radio interferometric imaging aims to reconstruct the sky x from the visibilities v by inverting the forward model. The formulation of Eq. 2.19, taking H = SF , is similar to many problems and applications in signal processing where the original signal must be recovered thanks to an inversion. In most cases, the objects or quantities we are interested in cannot be measured directly but through more or less sophisticated measurement tools. We can, for example, quote medical imaging (MRI, tomography, ...), geology (earthquake, ...), satellite imaging, or as in our case, astronomy.

A problem is said to be well-posed in the sense of Hadamard [Hadamard 1902] when it verifies the following three conditions:

• a solution exists,

• the solution is unique,

• the solution's behavior changes continuously with the initial conditions.

A problem is said to be ill-posed if one of the three previous conditions is not reached. A radio interferometer samples the Fourier transform of the observed sky. The geometry of the radio interferometer, the number of antennas in the array, and the observation time increase the number of samples making better coverage, see Fig. 2.6. However, the sampling will remain incomplete. As we have no information on the unsampled Fourier coefficients, an infinite number of solutions can complete these missing samples. Therefore, there is an infinite number of solutions to the imaging problem. Thus, this problem does not respect Hadamard's conditions. Reconstructing the observed sky from the visibilities is then an ill-posed inverse problem.

In the case of radio interferometric imaging, we consider knowing the probability distribution of the noise. Therefore, the density probability function of the likelihood is

p(v|x, Σ) = 1 (2π) M 2 |Σ| 1 2 × exp - 1 2 (v -SF x) † Σ -1 (v -SF x) , (2.20)
where † denotes the conjugate transpose. Note that the probability distribution in Eq. 2.20 is only valid in our case where the noise is decorrelated, and the real and imaginary parts of the visibilities have the same variance.

The maximum likelihood estimator is a popular estimator that maximizes p(v|x, Σ) when the data are known, and we want to estimate the source signal x. In radio astronomy imaging, these conditions are reached. Thus, the estimator becomes

x = arg max x p(v|x, Σ) = arg min x -log p(v|x, Σ) = arg min x 1 2 (v -SF x) † Σ -1 (v -SF x) = arg min x 1 2 ∥v -SF x∥ 2 Σ = arg min x J(x)
Consequently, we fall back on the weighted least squares solution. The argument of the minimum is explicit and can be obtained where the gradient of the objective function J(x) is zero. The gradient of J(x) is then

∇J(x) = ∇( 1 2 v † Σ -1 v -x T F † S † Σ -1 v + 1 2 x T F † S † Σ -1 SF x) = -F † S † Σ -1 v + F † S † Σ -1 SF x,
where the conjugate transpose of x ∈ R P is x † = x T . By setting the gradient of the objective function to 0, we fall back on the normal equation

(F † S † Σ -1 SF )x = F † S † Σ -1 v.
(2.21)

In the discrete case, the operators live in finite-dimensionnal space. Therefore, at least one minimizer x exists and dependent continuously on the measurement v [Idier 2001]. However, since the problem is ill-posed and has an infinite number of solutions, the term in parenthesis to the left of the equation cannot be inverted, and the problem cannot be solved. Therefore, the selection of a unique solution requires the addition of prior information on the search solution to be taken into account.

The addition of a priori information on the object of interest is characterized by an additional constraint in the optimization problem. Thus, the problem becomes (2.22) where in addition to the data fidelity criterion J(x), R(x) is the regularization function, and α is the regularization parameter that is set to find the best compromise between the data fidelity term and the regularization term. The regularization function is set to give a priori information on what is observed. In radio astronomy, it can be a constraint of positivity, a constraint on the presence of punctual or extended objects, or more. Therefore, choosing the right regulation function(s) is essential to ensure a good reconstruction [Thiébaut 2017]. Let us take an example where R(x) is quadratic, and therefore discontinuities will be excessively penalized :

x = arg min x J(x) + αR(x),
x = arg min x J(x) + α∥x∥ 2 .
(2.23) Thus, the regularized weighted least squares solution is

(F † S † Σ -1 SF + αI)x = F † S † Σ -1 v. (2.24)
As the number of visibilities and the number of pixels in the image to be reconstructed are very large, we consider imaging as a large-scale problem.

Let us define operator H = F † S † Σ -1 SF ∈ R P ×P as a large-scale matrix with a very high memory cost. Unless α! = 0 , the left-hand side of Eq.2.24 is invertible using iterative numerical methods [START_REF] Liu | [END_REF]]. However, even with these methods, the computation time of these algorithms for matrices with an order like H is prohibitive. Therefore, the explicit solution of Eq. 2.24 is not reasonably achievable. In order to minimize the regularized criterion of Eq. 2.22 or Eq. 2.23, it is necessary to use iterative algorithms, of which a partial list is described in the following sections.

One remark is that, even without considering DDEs, S † Σ -1 S is not diagonal. Therefore, F † S † Σ -1 SF is not a circulant Toeplitz matrix and, thus, cannot be considered as a convolution operator. Indeed, because of the convolutional interpolation operators S † and S, the matrix product S † Σ -1 S is not diagonal. Fig. 2.9 shows an example of an S † S matrix built with random visibilities for a 15 × 15 pixels image with Σ set as the identity matrix I. We can note that the matrix is rather a band diagonal with a Toeplitz aspect. The values outside the band diagonal are very low and depend on the convolution function used for the gridding and degridding. Therefore, there will be a small difference between the PSF taken for each pixel of the image.

CLEAN -the historical method

This section presents the general ideas behind the existing imaging algorithms. The CLEAN algorithm is the most widely used type of algorithm in the radio astronomy imaging community. This algorithm is separated into two families: the historical Single Scale CLEAN and the more recent multi-frequency multi-scale CLEAN. We will describe the idea behind these two methods before reviewing the other main imaging methods.

Single scale CLEAN

Several significant modifications have characterized the evolution of the Single-scale CLEAN algorithm. The first one was the most basic form of this algorithm [Högbom 1974]. A few years after, a new version improved the computational performance of CLEAN [Clark 1980]. Finally, a more sophisticated implementation, which is still used today, mitigates some approximations introduced by the first implementations [Schwab 1984].

The philosophy of CLEAN comes from the intuition that if we observe only spaced source points in the sky, then the dirty image is the convolution of the sky with the Point Spread Function (PSF) at the sources positions. With this assumption and Σ approximated as the identity matrix I, F † S † SF becomes a circulant matrix, such as

F † S † SF x ≈         a 1 a n .. a 2 a 2 a 1 .. a 3 .. .. .. .. a n-1 a n-2 .. a 0 a n a n-1 .. a 1                 x 0 x 1 .. x n-1 x n         (2.25)
Therefore, due to the diagonalization of the circulant matrices in Fourier space, we have

S † S ≈           c 1 0 .. 0 0 0 c 2 .. 0 0 .. .. .. .. .. 0 0 .. 0 0 0 0 .. c P -1 0 0 0 .. 0 c P           (2.26)
From this convolution approximation, we define the PSF y P SF as F † S † SF applied to a sky with only the central pixel set to 1.

y P SF = F † S † SF 0 1 (2.27)
where 0 1 is a vector containing zeros everywhere except the central pixel. This approximation is close to the truth if we consider the Van Cittert-Zernike theorem with a coplanar array of antennas. However, there are some limitations due to the convolution kernel used for the gridding and degridding operators, see Chap 3. From Eq 2.25 and Eq 2.27, we can write

F † S † SF x ≈ y P SF * x (2.28)
Thereby, we can approximate the gradient of the data fidelity term of the minimization problem [START_REF] Arras | [END_REF], Rau 2009], such as

∇J(x) = F † S † (v -SF x) = F † S † v -F † S † SF x ≈ y -y P SF * x (2.29)
The solution to the problem is to set the gradient to 0, which brings us back to the well-known notion that the Dirty image is obtained by convolving the real sky by the PSF.

y = y P SF * x (2.30)
The practical solution, assimilated to a sparse deconvolution and formalized by Högbom [Högbom 1974] is simplified in Alg. 1. The first step is to start with an empty sky model. Then, iteratively, we center the PSF on the brightest pixel of the dirty image. The PSF is multiplied by the flux value of the pixel, and a fraction of it is subtracted from the dirty image y. At the same time, the same fraction that has been removed from the brightest pixel is added to the model. This action is then repeated iteratively until the residual of the dirty image (also called residual image δy) is at the same level as the noise or at another tunable criterion. The PSF must be twice as large as the image for cases where the pixel is located on the edge of the image. This type of algorithm is similar to what we found under the name of Matching Pursuit [Mallat 1993a]. */

subtracting the model convolved by the PSF from the dirty image.

δy (k) = y -y P SF * x (k) (2.31)
In practice, the convolution is done in the Fourier domain to make it faster. An evaluation of the residual image is called a major cycle.

Finaly, the Cotton & Schwab variant [Schwab 1984], simplified in Alg. 3, is a generalization of Clark using the exact evaluation of the gradient, Eq. 2.29, rather than the convolution approximation. This has the effect of correcting the data when the antenna array is not coplanar and allows incorporating the w-term through wprojection [Cornwell 2008b]. On the other hand, it greatly impacts the performance and the time to evaluate a gradient. For this reason, it is desirable to make as few evaluations as possible. In the end, a mix of these three variants is used. A mix of Högbom and Clark is used to add components to the model. This step is also called the minor cycle. Furthermore, from time to time, the residual image δy is updated using the Cotton & Schwab variant.

Multi-scale and multi-frequencies CLEAN

These two algorithms are extensions of the CLEAN algorithm, which took into account only source points, allowing to better solve the problem with extended sources and multi-frequency images.

Multi-scale Clean [Cornwell 2008a, Offringa 2017] works in a similar way to CLEAN. Instead of using only one delta function to model a source point, there is a function dictionary to differentiate the cases. The implementations follow Schwab's idea of major and minor cycles, with only the minor cycle that is different. There are many variations in the implementation of multi-scale CLEAN in different imagers. The differences between these implementations come from how to choose and which kernels (Gaussian, tapered quadratic function, etc.) are used in the minor cycle.

Mutli-frequency CLEAN [Rau 2011[START_REF] Arras | [END_REF] allows to Deconvolve a multifrequency image. This algorithm allows adding the uv-coverages of several frequencies on the same grid. The advantage is to have better uv-coverage to decrease the problem's difficulty. In addition, this also has a gain in memory cost since it reduces the number of slices of the image cube. The principle is to consider the sky as a collection of source points with a Taylor Polynomial Spectrum.

Other methods

Even though CLEAN and all its variants are still widely used, much work and literature are done on alternative methods. Indeed, the advantage of CLEAN is that it is a well-known algorithm of which biases are known and controlled. On the other hand, there are also many weaknesses that these other methods are intended to atenuate. See below for a short overview of the main alternative methods.

Bayesian methods

The forward problem, Eq. 2.19 can be formulated as a Bayesian Statistical Inference problem. This formulation allows us to derive the estimate of x and provides tools to analyze and quantify the uncertainty in the solution. Considering a Gaussian noise, the likelihood function of the statistical model p(v|x) associated with the forward problem can be expressed as

p(v|x) ∝ exp -∥v -SF x∥ 2 2 2σ 2 .
(2.32)

However, as said previously, the imaging problem is ill-posed. Bayesian methods face this difficulty by adding a priori information about the object to be estimated. We thus have a prior distribution p(x) to regularize the problem, and, eventually, some additionally parameters ω.. Therefore, the prior and the observation are combined using Bayes' theorem to obtain the a posteriori distribution

p(x|v, ω) = p(v|x, ω)p(x, ω) p(v, ω) . (2.33)
The additionnal parameter ω represents any element,excluding data, that affect the expression of the probabilities, such as the noise. In Bayesian methods, the maximum a posteriori (MAP) solution gives the best image from the data [Thiébaut 2017], such as

x M AP = arg max x p(x|v, ω).
(2.34)

Several convex optimization methods can be used for efficient MAP estimation [Green 2015, Cai 2018b]. However, convex optimization algorithms do not provide uncertainty quantification. In radio interferometry, the problem being severely ill-posed, the uncertainty quantification becomes essential information. Traditional Markov chain Monte Carlo (MCMC) sampling methods provide uncertainty information. However, these methods are computationally expensive and cannot meet current radio telescope computation requirements. Recent literature proposes MCMC methods to reduce the computation time, as well as post-processing MAP estimation to add uncertainty quantification [Cai 2018a, Cai 2018b].

Compressed sensing methods

The Compressed Sensing (CS) theory acknowledges that a large variety of signals in nature are sparse. Indeed, signals measured are assumed to be sparse or compressible in some sparse basis. In the case of radio interferometric imaging, the main CS assumption is that the underlying signal has a sparse representation x = Ψα, with α ∈ C D contains only a few non-zero elements, and Ψ ∈ C P ×D a dictionary (e.g. a collection of wavelet bases).

[ McEwen 2011] reviewed that it results in an optimization problem called Basis Pursuit (BP). This problem aims to minimize the L1 norm of the coefficient α of the signal represented is a sparsity basis Ψ under a constraint of the L2 norm of the noise min

α ∥α∥ 1 s.t. ∥v -SF Ψα∥ 2 ≤ ε, (2.35)
where ε is a bound on th L2 norm of the noise n. The synthesis problem finds the image representation α with the final image obtained from x = Ψα [START_REF] Onose | [END_REF]].

Unlike synthesis-based problems, analysis-based problems recover the sky itself, solving min

x ∥Ψ † x∥ 1 s.t. ∥v -SF x∥ 2 ≤ ε, (2.36)
Recently, CS and convex optimization methods have been applied to deconvolution in radio interferometry and have shown encouraging results regarding the reconstruction quality [START_REF] Wiaux | [END_REF], Carrillo 2014, Carrillo 2012, Dabbech 2015].

The relationship between CLEAN and CS was first studied by [Marsh 1987]. Moreover, as pointed out by [START_REF] Wiaux | [END_REF]], CLEAN can be formulated in terms of a Matching Pursuit (MP) procedure [Mallat 1993b]. In this case, the MP algorithm will use a circulating dictionary corresponding to the convolution with the PSF.

Challenges of the new generation of radio interferometer

The amount of data generated by radiotelescopes is extremely large. The last generation, especially SKA, sees all the orders of magnitude of data rate or image size increasing significantly. The consequences are difficulties in processing, storing, and transferring data. A lot of work is beeing done to optimize each part of the data processing pipeline to improve the software implementation on the best suitable hardware.

The case of SKA

A significant technical enhancement characterizes the new generation of radio telescopes. These improvements have resulted in a larger collective area, greater sensitivity, and better spatial coverage. Consequently, these new radio telescopes enable us for making images of better quality, better resolution, and more spectral information. However, more complex and expensive algorithms are mandatory to correct data from these antennas and reconstruct images of the sky.

Beyond the computational difficulty of processing the data, the main problem remains in handling the amount of data generated by this new generation of antennas. With LOFAR, a SKA pathfinder, interferometric imaging observations can typically produce 35 TB/h of raw visibilities [van Haarlem 2013]. However, this problem increases on a different scale with the radio interferometer SKA. In the first construction phase of SKA, which should be finished in 2025+, the number of antennas, sensitivity, and spectral resolution will be the largest ever contemplated. For SKA1-Low, the number of antennas is expected to reach 130,000 [START_REF] Acedo | [END_REF]], and the generated signal rate will be 1.7 Pb/s.

At the base of these antennas, this data will be reduced through beamforming to be sent to the Central Signal Processor (CSP). The CSP is in charge of processing the data at a rate of 3 Tb/s to correct them and generate the visibilities. In these two stages, the calculations are carried out with Field Programmable Gate Arrays (FPGAs) and Graphical Processing Units (GPUs) while keeping track of the acquisition times using specific frames. Thus, with SKA1-Low, the data rate of the visibility generated by the CSP will be about 0.5 Tb/s to be transferred to the Science Data Processor (SDP) [START_REF] Dewdney | [END_REF]]. The SDP is located in the Supercomputer Center in Perth and Capetown. SDP software is in charge of processing visibilities to reconstruct the hyperspectral cube of the observed sky. Fig. 2.10 summarizes the flows of this first phase. With SKA1-Mid, the input rates should reach 9 Tb/s for the CSP and 6.8 Tb/s for the SDP.

HPC challenges

The SDP is in charge of the imaging and calibration part. However, the incoming data flow prohibits any medium and long-term storage capacity allowing only short-term buffering. Therefore, the visibilities must be treated in quasi-real time during the whole duration of acquisition and observation. In order to have the computing capacity to process information at this rate, the servers will require thousands/millions of computing cores distributed over several computing nodes. Visibilities will be read and written in Measurement Set (MS) format, which consists of a table-based format built after the Casacore Table Data System (CTDS) [Diepen 2015]. The buffers on the shared file system must support the simultaneous IO write/read of several million tables. Since the sky observation never stops, even during the day, the HPC server requires permanent data processing. With this type of continuous pipeline, the risk of failure of the entire HPC system due to the failure of a single worker is very high and must be avoided. Thus, a "Workflow Execution System" that allows the isolation of local errors and thus avoids a global failure is mandatory. In the context of SKA, DaLiuGe [Wu 2017] is the most popular and has been tested in almost real conditions by simulating the SKA data and then doing a first processing step on a few thousand nodes of the Summit HPC server [Wang 2020].

In addition to IO and memory capacity, the design of the computational capabilities is essential because of the high computational costs of imaging and calibration algorithms. Indeed, the calibration problem is computationally expensive and mandatory to correct strong image artifacts. Indeed, [Wijnholds 2014] has shown the high computational cost of this step and that the complexity increases with the number of antennas and the size of the FoV. Studies in the literature aim to improve the calibration robustness to reduce artifacts [Ollier 2018] and to improve the algorithmic methods to reduce the computational costs [Tasse 2014]. The second challenge, which is the core of this thesis, is the imaging problem. Indeed, most reconstruction methods in radio astronomy require iterative algorithms. However, most of these methods use the FFT operator in addition to the gridding and degridding operators to interpolate the visibilities in a uniform grid. As indicated in Eq 2.18, these operators scale with the number of pixels and the number of visibilities. In the context of SKA, we can easily see a bottleneck in these operations because of the amount of data to process. For example, the first phase of SKA, SKA-1 Mid, operating around 20GHz, will consist of 197 antennas. The order of magnitude of the images to be constructed will be greater than 2 30 pixels with up to 2 16 frequency channels [Scaife 2020]. With a standard observation of several hours, [Wijnholds 2014] has estimated that the computational cost to compute a dirty image would be P im = 356 PFlops.

To compare, the LOFAR radio telescope, which is already very large, requires only P im = 2 PFlops. Currently, no HPC system can handle this computing power in real time. Moreover, data backup is not possible either because of the data volume.

State-of-the-art

The reconstruction of the sky in the SDP stage is the most critical step of the pipeline because it has the highest computational cost and is I/O-bound. Indeed, the constraint of computing in quasi-real time means that the data flow coming to the SDP is currently higher than the computing capacity. Therefore, in order to reduce the computation time to generate an image, we can divide the problem into three subparts :

• The imaging framework. Also called imager, it manages the sky reconstruction algorithms and high-level hardware resources.

• The amount of data to process.

• The optimization of intensive computations on an accelerator platform.

Chapter 2: Problem setup and state-of-the-art

Imagers

Imaging frameworks aim to process the data as quickly as possible while making the best reconstruction of the sky from the visibilities. There are many imagers whose specificities often depend on the targeted radio telescope. Thus, none of them can handle SKA's computational specificities.

The most famous imager is the Common Astronomy Software Applications (CASA) [Jaeger 2008]. It is considered the state-of-the-art imager and is responsible for democratizing the MS data format. It allows multi-core CPU operation as well as multi-node use. However, it is limited by its lack of direction dependant effects corrections and by its lack of diversity in deconvolution algorithms. Moreover, it does not incorporate the latest accelerator platform implementations or more efficient algorithms. It is, therefore, slower and less accurate than other imaging frameworks. One of the most famous imagers is WSClean [Offringa 2014]. It is an imaging framework for large FoV correcting visibilities using the w-projection [Cornwell 2008b] and w-stacking algorithms [START_REF] Young | [END_REF]]. Moreover, it allows CPU multi-core parallelization and multi-node by using MPI. It has recently upgraded for corrections for direction-dependent effects on accelerator platforms such as GPUs and FPGAs [START_REF] Veenboer | [END_REF]]. DDFacet [Tasse 2018] is an imager for wide FoV and wide frequency band sky images. Thanks to its calibration framework KillMS, DDFacet is the only imager to take into account the direction-dependent Jones matrices and corrected visibilities from the third generation of calibration. During the reconstruction of the sky, the corrections are applied by dividing the image into several facets, considering that the correction to apply is constant within the facet.

The framework allows a multi-core CPU parallelization based on facets. The distribution of the computation on several nodes is the subject of Chap 4. Many other imagers have been developed over the year with smaller development teams. Thus, they are less maintainable, and fewer reconstruction options are available. Among these, some have improved image reconstruction quality [Carrillo 2014,Bester 2021], and others have speeded up the computation [Bao-qiang 2019].

Data averaging

One way to reduce the computing time during imaging is to reduce the amount of visibility computed. During the gridding operation, [Muscat 2014] proposed to avoid the convolution for visibility whose coordinates are too close to the previously processed visibility to avoid performing the same computation twice. Another method to reduce the number of computations, the most used and straightforward one, is to average all visibilities in time. For example, a 10s averaging will average all visibilities in values and coordinates in intervals of 10s to reduce the total number of data to process. However, this technique is known to decorrelate the data, causing artifacts to appear on the image when averaging too much [Thompson 2001]. Baseline Dependent Averaging (BDA) has improved this technique [Wijnholds 2018a]. BDA averages the visibilities according to the baseline. It will apply a stronger averaging for small baselines and a weaker averaging for large baselines. This averaging can be done using specific convolution functions [Atemkeng 2016b].

Hardware accelerator

An active branch of research is focused on the acceleration of algorithmic blocks used by most sky reconstruction algorithms. The interpolation block is the hottest topic because it is often computationally expensive, C grid = O(P log 2 P + C 2 supp M ). Most of these works focus on GPU implementations. Also a hot topic in MRI, this interpolation is done as a non-uniform Fourier transform [Zhili [START_REF] Yang | [END_REF]. In radio astronomy, GPU implementations are mainly done on the gridding algorithm to use the fast Fourier transform algorithm.

Generic image-processing convolution on GPU has been studied intensively for years and is even studied in many tutorials to learn programming. However, because visibilities are not sampled uniformly in radio astronomy imaging, the gridding and degridding are rather a sampled convolution than an image convolution. As a result, memory access patterns are less predictable, and thus creating an image is much more computationally expensive. Implementations on accelerators started with the democratization of the GPU for scientific computing around the 2010s.

Edgar et al. [Edgar 2010] have developed a GPU gridder, focused on NVidia platforms with Common Unified Device Architecure (CUDA) language, for the Murchison Widefield Array (MWA) telescope. Their implementation linked a CUDA thread for each grid point (a pixel). Then, each thread searches across all the visibilities which one has an impact on the corresponding grid point by checking them one by one.

Amesfoort et al. [START_REF] Van Amesfoort | [END_REF] focused on maximizing the Bandwidth (BW) device memory obtained. Therefore, in their implementation, each thread block writes to a private grid stored on the device memory to avoid atomic writes. However, this technique is limited by the size of the grids that can be used and the GPU memory.

Romein [Romein 2012] has introduced a more efficient work distribution strategy for gridding. This strategy was designed to minimize device memory access. The grid is divided into subgrids, where the size of each subgrid is the size of the convolution kernel used. We then create a number of threads equal to the grid points in the subgrid. Each thread carries a large number of grid points in the grid but only one in the subgrid. Then we assume that the convolution kernel moves slowly on the grid, meaning that the position of the visibilities varies slowly in time. Each thread can accumulate results on registers until the grid point is no longer in the subgrid, accumulating the result in the device memory.

Merry's method [Merry 2016] is based on Romein's for the GPU implementation. Each CUDA thread handles several grid points on the same subgrid. This

Shift

Fourier space

Image space

Shift

Fourier space Fourier space Figure 2.11: Simplified Image Domain Gridding illustration. First, a group of visibilities is shifted to low frequencies. Then, a DFT is performed on the shifted visibilities to make a low-resolution image. A Fast Fourier Transform is applied to the low-resolution image to make gridded visibilities. Finally, gridded visibilities are shifted back and added to their original place. implementation is done in the form of loop unrolling. Moreover, because of the large size of the convolution kernels, each subgrid is divided into several tiles. Threads first deal with one tile, and when it is finished, they move on to the next one.

More recently, [van der Tol 2018] introduced a novel gridding method, the Image Domain Gridding (IDG), which avoids the costly step of computing oversampled convolution kernels. A simplified explanation of gridding with IDG is illustrated by Fig. 2.11. First, a group of visibilities is partitioned in a subgrid where the position of the central pixel of the block is given by (u 0p , v 0p ). The block of visibilities is then shifted from (u 0p , v 0p ) to the origin of the Fourier grid, such as

V (u, v) = P p=1 δ(u -u 0p , v -v 0p ) × V p (u, v),
(2.37) with V p (u, v) the partitioned visibility. Because there is a low number of samples and the subgrid is much smaller than the main Fourier grid, the low-resolution image I is built using a DFT, such as

I = F -1 { V }. (2.38)
The low-resolution images are then transformed back to the Fourier space using an FFT. The result is shifted back and added to the large uv grid. IDG method makes DDEs correction costless. Moreover, the highly parallel structure of the algorithms makes it ideally suited for massively parallel hardware [Veenboer 2017[START_REF] Veenboer | [END_REF]].

Conclusion

Imaging in radio astronomy is a complex field characterized as a large-scale problem. Indeed, due to the size of the images to be reconstructed or the number of visibilities to be processed, the imaging phase and all the previous steps have a very high computational cost. In addition, the latest generation of radio telescopes tends to increase the rate of visibilities and the complexity of the corrections during the imaging and calibration phases. There are numerous imaging algorithms based on different mathematical optimization methods. They all have in common being iterative and involving computationally expensive operators such as gridding and degridding. A lot of work has been done to reduce the reconstruction computation time, especially on computational optimizations of gridding and degridding operators. With the computational challenges of the next generation of radio telescopes, illustrated by SKA, it is crucial to persist in these optimizations to reduce the computation time and the amount of data to process. The next generation will also have a lot of corrections to perform on visibilities. Therefore, reducing the computation time for imagers, currently offering decisive wide-field and wide-band correction, is essential. The work of this thesis will therefore focus on these two axes for the acceleration of imaging algorithms in order to reduce the computation time. One part will be on the low-level acceleration of specific algorithmic blocks, gridding and degridding, on the accelerator platform. The second part will be on the high-level acceleration of the DDFacet imaging framework to be distributed as efficiently as possible on different computational nodes.

Chapter 3

Fast Grid to Grid interpolation In this work, we address the problem of the computational cost of interpolation on a Fourier grid in the context of radio interferometry imaging. We have seen in chapter 2 that the new generation of radio telescopes makes the imaging step a bottleneck. The main reason is the iterative aspect of these imaging algorithms and the use of interpolation methods, gridding and degridding, for a considerable number of samples. The computational cost of these operators is proportional to the number of data and the size of the convolution kernel used by the interpolation. In order to decrease the computational time of this part of the imaging algorithms, it is necessary to tackle all the fronts: the interpolation algorithm and its implementation.

This chapter is organized into three sections. The first is a brief review of CPU and GPU architectures in terms of parallelization paradigm and memory hierarchy. In the second, we have explicitly included a publication submitted to "Astronomy & Computing". It provides details on the decomposition of the gridding and degridding operators and on the algorithmic implementation of the proposed Grid to Grid method. Moreover, we compare the GPU implementation of the method with the state-of-the-art. The last section provides the results of the CPU implementation of the G2G algorithm.

36

Chapter 3: Fast Grid to Grid interpolation

Multi-core and many-core architectures

Multi-core and many-core processors are microprocessors on a single integrated circuit that read and execute program instructions. The main difference is the number of computing cores allowing to execute concurrent instructions. A multicore processor has up to several dozen cores. A many-core processor has more than a hundred or even a thousand cores.

Multi-core CPU architecture

Historically, software development was pushed to do Single Instruction Single Data (SISD) on single-core architectures. The CPU executed a single instruction on a data stream at a given clock cycle. The concurrency tasks were left to the task scheduler, which divided the processor into several processes [START_REF] Akhter | [END_REF]]. From the 80s to the 2000s, increasing clock frequency and cache size was enough to improve performance without significant changes in the architecture. However, as illustrated in Fig. 3.1, the frequency is directly related to energy consumption. Therefore, limitations in semiconductor technology have limited further increase in clock rate. In the 1990s, multi-core processors began to emerge. The multi-core made it possible to overcome the limitation of increasing frequency with the disadvantage of a much more complicated design and heat dissipation problems. The cores working in parallel allowed to make instruction-level parallelism. Therefore, multi-core processors are essentially Multiple Instruction Multiple Data (MIMD), meaning that the different cores run different threads operating on different parts of the memory. Initially, MMX instructions only performed operations on integers. Then the Streaming SIMD Extensions (SSE) and Advanced Vector eXtensions (AVX) allowed operations on float in vectors of 128, 256, and 512 bits registers.

Multi-core and many-core architectures

Many-core GPU architecture

Historically developed in the 90s and used for 3D rendering, the GPU was democratized in the 2000s. At this time, it was considered a General Purpose GPU (GPGPU) and a wide adoption by the scientific community came at the in the 2010s. The CPU is driven to process large sequential programs per CPU core with complex branching and wide memory access patterns. Following the SIMD paradigm, the GPU is driven to process simple operations on many elements. The GPU is then specialized for highly parallel computations such that more transistors are devoted to data processing rather than data caching and flow control [START_REF] Nvidia | Compute unified device architecture programming guide[END_REF]]. Fig. 3.2 illustrates the distribution of arithmetic resources and memories for a CPU and a GPU. The GPU dedicates much more transistors for processing, which is preferable for highly parallel computing. The memory accesses are supposed to be hidden by the calculations. Modern GPUs work very well in a context where most of the work is parallel. The CUDA nomenclature divides the work into a grid of several thread blocks. As illustrated in Fig. 3.3, the memory hierarchy follows the fineness of the elements. For a thread, the dedicated memory is in the form of a local register. Each block has access to a shared memory common to all threads of the block. Finally, global memory, or device memory, is common to all blocks, but the latency of memory accesses is much longer. In recent years, large projects involving astronomers, computer scientists, and engineers have developed new generations of radio telescopes to improve antenna sensitivity, resolution, and data quality. These improvements have led to a significant increase in data generation, involving a higher volume of data to process. Because the volume of data is directly related to the number of antennas, telescopes with many antennas will have more data to process.

Grid of blocks

Imaging is a critical phase of the data processing pipelines. It consists of adding the data generated by the telescope, called visibilities, on a Fourier-transformed grid to create an image of the sky. This step, called gridding, and its adjoint, the degridding, are critical steps in the imaging phase because they are computationally very expensive. Existing methods for sky reconstruction have all the common points of being iterative and require to grid and degrid the visibilities at each iteration. Even if algorithms like Cotton-Schwab CLEAN (Schwab, 1984) limit the number of these operations, the computational cost remains extremely high.

These gridding and degridding operators are good candidates to be parallelized on many-core accelerators like the GPU. Therefore, a lot of work has been done, mainly on the gridding operator, to decrease the main memory bandwidth and not rely on memory caches. Several works, such as (Merry, 2016), (Romein, 2012a), and (Muscat, 2014), used the advantage of the proximity of the coordinates between two data samples on the grid to reduce the memory accesses. However, most of the existing methods only focus on the gridding step.

This paper presents a method to reduce the algorithmic complexity of gridding and degridding operators by merging them to reduce the amount of data used. This method, based on FHD [START_REF] Sullivan | Fast holographic deconvolution: a new technique for precision radio interferometry[END_REF], also has the advantage of reducing the global memory footprint of these operations. Finally, this paper presents a GPU implementation of this new operator, taking into account the extension of this work with the w-correction, as well as the addition of a degridding operator implementation.

This paper is structured as follows. Section 2 presents the basics of imaging in radio astronomy. Section 3 explains in a theoretical way the Grid to Grid (G2G) method. Section 4 is dedicated to the parallelization of the method on GPU, compared to the state-of-the-art on the same subject. Section 5 shows the performances of this method on GPU. Finally, we conclude this work by discussing future works.

Radio Interferometric Imaging

Forward problem and synthesis imaging

A radio interferometer array is a network of antennas that generates measures of the radio emission of the observed sky x. Each antenna's pair is defined with a baseline b = (u, v, w) where u, v, and w are the coordinates in unit of wavelength λ. The measurement of one antenna pair of baseline b, or so-called visibility, is defined as

v(u, v, w) = x(l , m, n) n e -2iπ(ul +vm+w(n-1)) dl dm (1)
where x is the sky brightness distribution using a coordinate system (l , m, and n = √ 1 -l 2 -m 2 ) that indicates an angular position.

As the earth's rotation causes the baseline to change over time, the position of generated visibilities in the UVW space also changes during the time. Thus, the total number of visibilities generated during an observation is

M = N bl × N ch × N pol × T ∆t . ( 2
)
where N bl is the number of baselines, N ch the number of frequencies, N pol the number of polarization, T the duration of the observation, and ∆t the integration time for a sample. As current radio telescopes, such as SKA, have many antennas (197 for SKA-mid) and have to produce large images up to 30.000 × 30.000 pixels, the radio interferometric imaging problem is considered a large-scale problem.

If the array of antennas is coplanar (all visibilities fall on an arbitraly plane in the uvw space) and the FOV is small (n ≈ 1), each visibility v(u, v) is the 2D spatial Fourier transform of the sky distribution at frequency (u, v). This result is known as the van Cittert-Zernike theorem (Thompson et al., 2001), and Eq. (1) becomes v(u, v) =

x(l , m)e -2iπ(ul +vm) dl dm.

(3)

Since the array cannot cover the full (u, v) plane, radio interferometric imaging aims to recover the sky from incomplete visibility measurements, leading to an ill-posed linear inverse problem.

Imaging synthesis computes the so-called dirty image, which is the inverse Fourier Transform of the measured visibilities such as

y(l, m) = M i=1 v i e 2iπ(u i l+v i m) .
(4)

The most straightforward approach to do it, using a Direct Fourier Transform (DFT), is computationally too expensive. The computational cost of the DFT for M visibilities in N p × N p pixels image is

C DFT = O(N 2 p M) (5)
On the other hand, the FFT needs the visibilities to be uniformly sampled, which is not. The most common solution is to interpolate visibilities with non-uniform uv-coverage into a uniform one. This procedure is called gridding. Its computational cost for coplanar arrays and the FFT is

C grid = O(N 2 p log 2 N p + C 2 supp M) (6)
where C 2 supp is the size of the convolution function in 2D. After discretization, the dirty image can be represented by N p × N p = P pixels grid as y ∈ R P , and the visibility measurement as a vector v ∈ C M . Relation between y and v is described by

y = F † S † v (7)
where F † ∈ C P×P is the inverse Fourier transform matrix, S † ∈ C P×M maps M visibilities to P Fourier coefficients and is the gridding operator. The discrete forward problem is the adjoint of eq. 7 such as

v = SFx + n (8)
where F ∈ C P×P is the Fourier transform, S ∈ C M×P is the degridding operator, adjoint of S † , and n ∈ C M is an i.i.d. Gaussian noise.

Gridding operator

To use the FFT algorithm, the gridding (and similarly the degridding) operator aims to interpolate data with non-uniform coordinates on a uniform grid. This kind of problem is common in different fields of study. This method is, for example, known as the Non-Uniform Fast Fourier Transform (NUFFT) in the medical domain, such as MRI [START_REF] Fessler | Nonuniform fast fourier transforms using min-max interpolation[END_REF], [START_REF] Barnett | A parallel non-uniform fast Fourier transform library based on an "exponential of semicircle[END_REF], (Zhili Yang and Jacob, 2009).

In such interpolation, the true position b i = (u i , v i ) of the visibilities are approximated to the nearest neighbor on an oversampled grid of resolution (∆ u , ∆ v ). If K is the oversampling factor, the FFT steps of the grid g are (K∆ u , K∆ v ). The approximated position of b i becomes

b i = (p i ∆ u , q i ∆ v ). ( 9
)
In that case, using the gridding operator S † to compute the coordinate b k = (p k K∆ u , q k K∆ v ) of the grid g is equivalent to

g(b k ) = M i=1 C † b k -b i v i = M i=1 C † p k K∆ u -p i ∆ u , q l K∆ v -q i ∆ v v i ( 10
)
where C † is a discrete precomputed 2D kernel of size C supp × C supp in FFT step which contains K 2 ×C supp ×C supp values. The choice and the size of this interpolator to avoid aliasing effects and maximize accuracy is a well-known subject in the literature, see for instance [START_REF] Fessler | On NUFFT-based gridding for non-cartesian MRI[END_REF], [START_REF] Beatty | Rapid gridding reconstruction with a minimal oversampling ratio[END_REF], [START_REF] Thévenaz | Image interpolation and resampling[END_REF]. In our case, we used the Kaiser-Bessel function with a support C supp = 7.

Sky estimation

The sky x reconstruction for incomplete data v is an ill-posed inverse problem. Therefore, interferometric imaging is usually defined as

x = arg min x 1 2 ∥v -SFx∥ 2 + R(x) (11)
where ∥v -SFx∥ 2 is the data fidelity term and R a regularizer [START_REF] Giovannelli | Positive deconvolution for superimposed extended source and point sources[END_REF]Thiébaut and Young, 2017;[START_REF] Bester | A practical preconditioner for wide-field continuum imaging of radio interferometric data[END_REF]. Algorithms used to solve this optimization problem require computing many gradients of the criterion during their iterative processes. The computing cost of a gradient of the data fidelity term

∇J(x) S T D = F † S † (v -SFx) (12)
can be very high since it requires the computation of the gridding and degridding. The gradient computation ∇J(x) S T D corresponds to the residual image computation δy of the historical Cotton-Schwab CLEAN algorithm (Schwab, 1984). We proposed in this paper a new formulation to reduce the computing cost of the gradient evaluation without additional approximation. The gradient computation corresponds to the so-called "major loop". The deconvolution, called "minor loop," used in the overall reconstruction methods like (Schwab, 1984) Cotton-Schwab CLEAN is not studied here.

Grid to Grid interpolation

In this section, we propose a decomposition of the gridding and degridding operators in order to facilitate their fusion into the Grid to Grid operator. We also present a generalization and the particularities of the G2G method.

Gridding and degridding decomposition

A way to implement the gridding operator S † is to decompose it into three sub-operators : accumulation, convolution, and subsampling.

1. Accumulation A † .
Let's define g ′ a vectorized oversampled 2D grid of size KN p × KN p = K 2 P with thin pixel resolution (∆ u , ∆ v ). The oversampled pixel value associated with the coordinate b k on the thin grid is

g ′ (b k ) =        i v i if b i = b k , 0 otherwise.
Thus, g ′ is a grid of accumulated visibilities that are approximated at the same position on the thin grid. We therefore define an accumulation operator

A † ∈ {0, 1} K 2 P×M such as g ′ = A † v. A † is a sparse operator with one ′ 1 ′ per column, ′ 0 ′ otherwise. 2. Convolution C † .
Considering a 2D grid g of size KN p × KN p with thin step resolution (∆ u , ∆ v ) (same size as g ′ ). This step is a twodimensional discrete convolution between the thin grid g ′ and the convolution function C † , such as the value associated to the coordinate b k is

g(b k ) = i j C † (p k -p i )∆ u , (q k -q j )∆ v g ′ (b i j ).
Thus, we define the discrete convolution operator

C † ∈ C K 2 P×K 2 P such as g = C † g ′ . 3. Sub-sampling O † .
The FFT grid is coarser by a factor K regarding the oversampled grid. We define the subsampling operator O † ∈ {0, 1} P×K 2 P that reduces a grid of size K 2 P into a grid of size P. Thus, the FFT grid g is g = O † g. The subsampling operator has one ′ 1 ′ per line, ′ 0 ′ otherwise.

Combining these three sub-operators, we can define the FFT grid as g = O † C † A † v. Thus, the gridding operator can be written as

S † = O † C † A † . (13)
This decomposition of gridding into sub-operators is illustrated in Fig. 1. In practice, these operators are not instantiated as a matrix but are implicit, being applied by dedicated code.

The degridding operator S is the adjoint of the gridding S † . Thus, from Eq. ( 13), we can define the degridding decomposition as

S = ACO, (14)
where O ∈ {0, 1} K 2 P×P is an oversampling operator filling the oversampled grid with zeros, C ∈ C K 2 P×K 2 P is a 2D discrete convolution where the convolution function C is the flipped conjugate of C † , and A ∈ {0, 1} M×K 2 P is a sampling mapping operator that builds M visibilities at their approximate coordinates from K 2 P coefficients. If g is the Fourier grid to be degridded, the visibility vector v is v = ACO g. This decomposition is illustrated in Fig. 2.

Figure 1: Gridding decomposition. First, accumulation of visibilities on a 2D thin grid of size KN × K M, then a 2D discrete convolution with a kernel C † , and finally, a K subsampling.

Figure 2: Degridding decomposition. First, oversampling of a map g with a factor K, then a 2D discrete convolution with a kernel C, and finally, a mapping operator A maps the Fourier coefficients to M visiblities.

Gradient computation with G2G

The standard method to compute the gradient operator described by Eq. ( 12) can be developed such as it becomes

∇J(x) = F † S † v -F † S † SFx = y -F † S † SFx. ( 15
)
The Fast Holographic Deconvolution (FHD) [START_REF] Sullivan | Fast holographic deconvolution: a new technique for precision radio interferometry[END_REF] algorithm used this decomposition. We can note that FS † v is the dirty image y, which is computed only once. In this method, they build the holographic map H ∈ C P×P such as H = S † S. This method has the advantage of being fast because a simple matrix product is needed to compute the gradient. However, the computation time and the memory footprint required to compute H are prohibitive.

In our case, we use the same development of Eq. 15, as well as the gridding and degridding decomposition. Thus, the succession of degridding and gridding operators can be rewritten as

S † S = O † C † A † ACO. ( 16
)
from this equation, we built A * = A † A, and A * ∈ N K 2 P×K 2 P a diagonal matrix, such as a = diag( A * ) is a vector of size K 2 P.

The value of each element a i of a corresponds to the number of visibilities whose positions have been approximated on the corresponding oversampled pixel. Therefore, it is sufficient to do a term-wise product between g and a instead of going through the visibility vector v. This method, illustrated by Fig. 3, will be referred as the Grid to Grid method (G2G), and the gradient computation becomes

∇J(x) G2G = y -F † G2G(Fx), (17)
where

G2G(•) = O † C † A * CO• .
In a typical continuum imaging application, the STD method uses a different number of subbands for gridding and degridding. Typically, gridding requires fewer subbands than degridding. However, with the G2G method, the number of subbands must be identical and set as the greater of the two.

G2G compression ratio

We define M ′ as the number of non-zero elements in a. With M the number of raw data, practical cases show that M > M ′ or M ≫ M ′ depending on the oversampling factor K. Consequently, the G2G method avoids processing the vector of M raw visibilities but instead processes a smaller vector of size M ′ . Moreover, unlike FHD [START_REF] Sullivan | Fast holographic deconvolution: a new technique for precision radio interferometry[END_REF] who builds and stores S † S ∈ C P×P as a matrix, we only need to build and store the sparse vector a. Thus, from Eq. ( 10) and Eq. ( 16) the output grid g at coordinates b

k = (u k , v k ) becomes g(b k ) = M ′ i ′ =1 C † b k -b i ′ a i ′ i, j C b i ′ -b i j g(b i j ). ( 18
)
For the same oversampling factor, the gradient computed by the G2G method Eq. ( 15) and the STD method Eq. ( 12) is the same. Furthermore, the computing cost of the gradient using the G2G method is given by

C G2G = O(2M ′ C 2 supp + M ′ + 2P log P + P) (19)
while the computing cost of the standard approach is

C S T D = O(2MC 2 supp + 2P log P + M). (20)
Hence, for the same K, the presented method is more efficient without additional errors when

M ′ < 2MC 2 supp + M -P 2C 2 supp + 1 (21)
With a data set for a standard sky observation (8h), empirical measurements show that MC 2 supp >> P -M. Thus, Eq. 21 can be simplified such that

M ′ ≲ M. (22)
In cases where MC 2 supp >> P -M is false, we can no longer simplify, and Eq. 21 must be respected to ensure a decrease in the computational cost.

We define as the compression ratio, the ration between M ′ and the raw number of visibilities M, such as

c α = M ′ M × 100. ( 23
)
Figure 3: Degridding and gridding decomposition using the G2G method. An element-wise multiplication is done between the diagonal of the operator A * and the oversampled grid.

Note that the compression provided by the G2G method can be seen as a generalization of the compression provided by the Baseline Dependent Averaging method [START_REF] Wijnholds | Baseline-dependent averaging in radio interferometry[END_REF][START_REF] Atemkeng | Using baseline-dependent window functions for data compression and field-of-interest shaping in radio interferometry[END_REF]. The BDA method averages visibilities over time as a function of baseline, whereas G2G averages visibilities over their position in the Fourier plane. Moreover, BDA's compression leads to an approximation that can be chosen for higher or weaker compression. With G2G, the compression comes directly from the model.

Extension for w-term

Taking into account the effects related to the non-coplanarity of the antenna array (visiblities fall on different planes in the uvw space), the visibility equation from Eq. 1 is v(u, v, w) =

x(l , m, n) n e -2iπ(ul +vm) G(l, m, w) dl dm ( 24) with G(l, m, w) = e -2iπ w( √ 1-l 2 -m 2 -1) .

(25) [START_REF] Cornwell | The non-coplanar baselines effect in radio interferometry: The w-projection algorithm[END_REF] has shown that a reprojection is allowed to and from any position in the (u, v, w) space to and from w = 0 plane by convolution using a known kernel. This technique, called w-projection, uses the convolution theorem from the Eq. 24 such as

v(u, v, w) = v(u, v, w = 0) * G(u, v, w), (26)
where G(u, v, w) is the Fourier transform of G(l, m, w). This convolution can be done during the gridding, which makes the convolution function more complex. In order to avoid computing the convolution function on the fly for each w value, the closest w-layer is chosen. A vast number of convolution functions are precomputed and stored. Moreover, the support size of the kernel depends on the image size and the w value [START_REF] Tasse | Applying full polarization A-Projection to very wide field of view instruments: An imager for LOFAR[END_REF]. Thus, the computational cost of the convolution increases regarding this particular correction. The G2G interpolation method has no significant changes by including the w correction. Indeed, using the decomposition used in Fig. 3, only the discrete convolution operators C † and C change depending on the convolution function, which is related to w. As illustrated in Fig. 4, a matrix A * must be computed for each layer W, such that A * w i matches all w values approximated at w i .

As shown in Fig. 5, the values of w vary smoothly over time. The construction of the sparse diagonals of A * w i is, therefore, not a problem.

Compression approximation

The plots of the variation of baseline length in the (u, v, w) space are curves. The projections of these curves in the (u, v) space become arcs. The u, v, and w coordinates are given from the cartesian coordinates system (X, Y, Z), a terrestrial reference frame that never change from a local observer on earth, of the baselines as

         u v w           = 1 λ           sin H 0 cos H 0 0 -sin δ 0 cos H 0 sin δ 0 sin H 0 cos δ 0 cos δ 0 cos H 0 -cos δ 0 sin H 0 sin δ 0                     X Y Z           (27)
where (H 0 , δ 0 ) are the coordinates of the observation's phase center on the celestial sphere. It has been shown by (Thompson et al., 2001) that the track formed in time by the baseline projected in the (u, v) space is an ellipse. Thus, by transforming Eq. 27 and removing the time aspect, we obtain the following ellipse equation

u 2 +       v -Z λ cos δ 0 sin δ 0       2 = X λ 2 + Y λ 2 . ( 28
)
The portion of the ellipse traced during the observation depends on the azimuth, the elevation and the latitude of the baseline, the declination of the source, and the hour angle covered during the observation. Fig. 6 shows the superposition of the uv track and the ellipse built with the same parameters. We can estimate the ellipse's perimeter formed from Eq 28 by taking the semimajor axis a =

Samples

X 2 λ + Y 2 λ and the semiminor axis b = sin(δ 0 ) X 2 λ + Y 2 λ .
The different formulas of Ramanujan [START_REF] Villarino | Ramanujan's perimeter of an ellipse[END_REF] allow us to approximate the ellipse's perimeter P with very low errors depending on a and b. We can, therefore, estimate the theoretical compression factor c α as a function of the ellipse's perimeter for an observation time P ∆H and the resolution ∆ u that depends on the oversampling factor.

In the case of a single baseline i at a specific frequency, the idea is to divide the number of intervals in the portion of the ellipse's perimeter P ∆H,i ∆ u by N ri , the number of visibility acquired during the observation period for this specific baseline. As the number of intervals might be higher than the number of samples, meaning no compression, the compression number is limited by 1. Thus, the compression estimation is

c α,i = min P ∆H,i ∆ u N ri , 1 × 100. ( 29
)
To estimate the compression for all baselines, take the average of all c α,i such as c α = mean( c α,i ).

(30)

GPU Implementation

In this section, we will introduce in more detail the state of the art of different GPU implementations of gridding algorithms. We will then see the specificity of the implementation of the G2G method on GPU.

Related work

Generic image-processing convolution on GPU has been studied intensively for years and is even studied in many tutorials to learn programming [START_REF] Cheng | Professional CUDA c programming[END_REF]. However, in radio astronomy imaging, the gridding and degridding are rather a sampled convolution than an image convolution. As a result, memory access patterns are less predictable, and thus, creating an image is much more computationally expensive.

GPU gridding implementations started with the democratization of the GPU for scientific computing around the 2010s. (Edgar et al., 2010) have developed a GPU gridder, focused on NVidia platforms with CUDA language, for the Murchison Widefield Array (MWA) telescope. They recognized that accumulating visibilities convoluted on the grid was not thread safe. Moreover, adding the visibilities directly on the grid significantly increases the device-memory access cost. Their implementation linked a CUDA thread for each output grid point (a pixel). Then, each thread searches for all the visibilities which one has an impact on the corresponding grid point by checking them one by one. However, this approach was inefficient, and only 60 visibilities out of 130,000 were concerned. To overcome this problem, they sorted the visibilities according to the (u, v) coordinates and grouped them into bins so the threads could search for the right visibilities more efficiently. (Van Amesfoort et al., 2009) focused on maximizing the bandwidth device memory obtained. Therefore, in their implementation, each thread block writes to a private grid stored on the device memory to avoid atomic writes. However, this technique is limited by the size of the grids that can be used and the GPU memory. (Romein, 2012b) is considered as the reference for recent GPU implementations and introduced a more efficient work distribution strategy for gridding. This strategy was made to minimize device memory access. The grid is divided into subgrids, where the size of each subgrid is the size of the convolution function. For each CUDA block, the number of threads created is equal to the number of grid points in the subgrid. Each thread is focused on a large number of grid points in the grid but only one in the subgrid. Because visibilities coordinates move slowly in time for the same baseline, looping over visibilities in the CUDA kernel also makes the convolution function move slowly. Each thread accumulates the results of the convolved visibility on registers until the grid point is no longer in the subgrid. When the grid point is out of the subgrid, the register is added to the device memory using an atomic add operation. Thus, this strategy minimizes writing in the device memory by using local registers.

The main improvement of (Merry, 2016) is the introduction of thread coarsening for large convolution functions. Thread coarsening is similar to loop unrolling but applied to parallel work items instead of sequential ones. In its implementation, each thread handles several adjacent grid points within a subgrid. Thanks to thread coarsening, memory accesses within a thread improved because the data for each grid point is close to each other in the device memory. However, this method requires significantly increasing the number of registers per thread.

Finally, the IDG also made a GPU implementation (Veenboer et al., 2017). This is one of the few works proposing implementations for gridder and degridder. However, the IDG algorithm is fundamentally different from the classical gridder. The principle is to do many NUFFT on small grids with very little data each time (van der Tol et al., 2018). In this case, we are not doing a sample convolution on the image but rather a succession of sine, cosine, and shift. This method, therefore, allows DDEs (primary beams, ionospheric screens, etc.) to be corrected transparently and without additional computational cost.

Grid to Grid sequential details

First, we present a sequential implementation of G2G that serves as the basis for the CPU code. Alg. 1 is the sequential pseudo-code for working with data and grids for one polarization. The diagonal diag(A * ) is stored in an array nhit which contains the value of each non-zero element a i with its corresponding coordinates on the diagonal. The getNhit function get the information of the A * operator through the nhit array regarding the for loop iteration.

The algorithm processes sequentially, in the order of storage, the M ′ non-zero elements diag( A * ). The distance between two consecutive samples for the same baseline is very small during an observation. Therefore, the distance between diag( A * ) elements is also very small, which allows a low spatiality between two iterations.

Grid to Grid GPU implementation details

we present the work distribution strategy, which implements the G2G method with a grid in the Fourier plane as input and a new grid in the Fourier plane as output. Our STD and G2G GPU implementation are based on Romein's previous implementation (Romein, 2012a) to reduce device-memory access. Additionally, we implement a reduction method to decrease the synchronizations between threads by optimizing the SIMD synchronous instruction within a warp. The G2G method inputs a Fourier grid and outputs another Fourier grid. On the other hand, gridding takes a visibility vector as input and outputs a Fourier grid. The degridding takes a Fourier grid as input and produce visibilities.

The strategy is as follows. We decompose the input and output grids into subgrids, where each subgrid has the size of the convolution function. In the example of Fig. 7, the grids are of size 12 × 12, and the convolution function is 4 × 4. In reality, both are much larger. We create a number of threads equal Algorithm 1: G2G -sequential pseudo code.

Data: (...,K, Igrid, Ogrid, nhit)

for i in 1:M' do u, v, w, a i , ch = getNhit(i, nhit) ; c i = (0,0); for u tap ← [-hal f size] to [+hal f size] do for v tap ← [-hal f size] to [+hal f size] do convU, convV = getConv(u tap , v tap , K); wc = convFctConj[w][convU][convV] ; c i += Igrid[u + u tap ][v + v tap ][ch] × wc; end end c i = c i × a i ; for u tap ← [-hal f size] to [+hal f size] do for v tap ← [-hal f size] to [+hal f size] do convU, convV = getConv(u tap , v tap , K); w = convFct[w][convU][convV] ; Ogrid[u + u tap ][v + v tap ][ch] += c i × w; end end end
to the number of grid points covering the size of the convolution function, 16 threads in this case (the case of larger kernels will be discussed later). Each thread manages the grid points assigned to it. In our case, a thread manages the ⋆ points and another the • points, which are equivalent grid points on the input grid and the output grid, as illustrated by Fig. 7. As seen in Sec. 4.2, the coordinates of two consecutive a i elements of the nhit array are close to each other. Therefore, the convolution function moves slowly as well, and the grid points of the input grid can be stored in registers. In the same way, the output grid points are first accumulated in registers before being updated in global memory once the thread has to deal with another grid point. With a convolution function of size n × n, this method allows reading and writing with atomicAdd instruction only 2 × n grid points simultaneously. Indeed, in the case of a diagonal move with visibilities close to each other, at most one row and one column will come out of the new grid points, so there would be 2 × n reading and 2 × n writing.

In order to accumulate the value of the compressed visibility on registers, it is necessary to reduce the portion of visibility calculated by each thread. However, the reduction has a low parallelization potential and can become a bottleneck in our algorithm. To minimize the time spent on this reduction, we need a number of threads equal to a power of 2. CUDA has designed several types of methods to optimize the reduction operator [START_REF] Harris | Optimizing parallel reduction in cuda[END_REF]. We have chosen a strategy to optimize this reduction to minimize the mandatory synchronization steps.

The hardware limits the number of threads per block and is often 1024 threads. The W-projection algorithm requires convolution functions of size C 2 supp > 1024. In this case, like (Merry, 2016), each thread handles several grid points adja-+ Figure 7: Mapping of work item matching the input and output grids. Each CUDA thread manages the same grid point in the convolution kernel area for the input and output grid. For example, one CUDA thread is in charge of all the grid points tagged by •, and another CUDA thread by x.

cent within a subgrid. Moreover, the memory accesses are done in loop unrolling to benefit from better memory access for the adjacent grid points. However, the increase in the number of registers prevents us from doing full thread coarsening.

We used a separable approximation of the convolution function. Moreover, like Merry and Romein, we assumed that the convolution function was polarization-independent. The results can, therefore, not be directly applied to the A-projection algorithm. A tunable factor is the number of data of diag A used per thread block. Arbitrarily, we have used 1024 elements per block. As Merry pointed out, higher numbers make more use of spatial coherence between adjacent elements but reduce parallelism.

A pseudo-code of the G2G GPU implementation is represented by Alg. 2. This simplified code shows the case with a single polarization. The pseudo-kernel loops the 1024 diag( A * ) elements, such that Nhit is an array that contains the coordinates and the non-zero values of the diagonal elements.

It is important to note that the STD gridder, STD degridder, and G2G implementations all have the same level of kernel optimization. Thus, comparing the performance of these implementations in the result part can be done without bias.

Results

This section aims to show the performance of the G2G algorithm implemented on GPU and CPU. We will study the accuracy and the memory footprint related to the oversampling factor. Then, we will see the performances of the implementation on CPU and GPU before comparing it to other state-of-Algorithm 2: Kernel G2G -pseudo code.

Data: (..., Igrid, Ogrid, nhit) in i = (0,0) ; add i = (0,0); c i = (0,0);

for i in n do u t , v t ,
in i = Igrid[gridU][gridv] ; end w = convFct[w][th convU][th convV]; wc = convFctConj[w][th convU][th convV]; c xx = reduce(in i ×wc); add i = c i × w × a i ; end atomicAdd(Ogrid[prev gridU][prev gridV], add i );
the-art algorithms. Finally, we will study the roofline model of our implementation in order to show the different possible improvements.

Experiment setup

For the GPU part, the experiments have been done on a BullSequana XH2000 with an NVIDIA A100-SXM GPU using CUDA. The number of CUDA threads per block used depends on the size of the convolution kernel. Each thread block handles 1024 visibilities, so the total number of blocks is the number of visibilities, 4 000 000, divided by the number of visibilities per block 1024 (as the result is not an integer, the last block deals with less than 1024 visibiltites). For the CPU part, the experiments have been done on BullSequana XH2000 system based on Intel(R) Xeon(R) Platinum 8358 32 cores dual socket. Thus, the code was parallelized on 64 cores.

We will not study the pre-processing necessary for constructing diag( A * ), as it can be done while building the dirty image. Nevertheless, we found that pre-processing performance to make the diag( A * ) matrix is highly dependent on ordering the data by baseline, as in any other gridder implementation.

The data set used a simulation of a 8h observation with the VLA-D with a 1s integration time for 64 frequency channels. The observation has 4 000 000 quad-polarized raw visibilities per frequency channel. The size of the grid used is 1280 × 1280 pixels following the CASA tutorial1 . The G2G method is advantageous for compact configurations like VLA-D or Nenu-FAR. It would be less beneficial for instruments like VLBI.

Oversampling factor K accuracy

The performance of the G2G method is directly related to the oversampling factor K used for the resolution of the oversampled grid, the uv coverage, and the convolution function. For the same oversampling factor, the gradient computed with the STD method is the same as one with the G2G method, such that ∇J G2G (x K= j ) = ∇J S T D (x K= j ).

(31)

Thus, the computational cost is reduced without depreciating the quality of the computed gradient. First, we looked at the accuracy related to the oversampling factor. The oversampling factor K = 64 is the default oversampling factor used in the WSClean imager [START_REF] Offringa | WSClean: an implementation of a fast, generic wide-field imager for radio astronomy[END_REF]. We will therefore take the gradient from Eq. 15 computed for this K as a reference. Fig. 8 and9 show the maximum error

MaxAE(K) = max( ∇J G2G (x K=64 ) -∇J G2G (x K= j )|) (32)
and the mean absolute error between gradients

MAE(K) = 1 P P i=0 |∇J G2G (x i,K=64 ) -∇J G2G (x i,K= j )|. ( 33
)
There is a direct impact between the reconstruction quality of the gradient and the oversampling factor. The smaller K is, the stronger the approximation to the true coordinates of the visibilities will be. But we know that a significant approximation tends to decorrelate the data and introduce artifacts in the image (Thompson et al., 2001). It is a known and already observed result [START_REF] Offringa | Precision requirements for interferometric gridding in the analysis of a 21 cm power spectrum[END_REF]. Moreover, the error is also dependent on the chosen convolution function.

The problem of accuracy for recent radiotelescopes such as LOFAR or MeerKAT, and soon SKA, is general for both the G2G and STD methods. Indeed, these two methods achieve identical accuracy for the same oversampling factor as they use the same convolution function for the interpolation. New gridding techniques have recently been developed to increase accuracy with a very low oversampling factor [START_REF] Ye | High accuracy wide field imaging method in radio interferometry[END_REF][START_REF] Barnett | A parallel non-uniform fast Fourier transform library based on an "exponential of semicircle[END_REF].

Memory footprint

We compare the memory footprint between the G2G method Eq. 17 and the STD method Eq. 12. In the case of STD, the memory footprint is mainly composed of visibilities for all polarity at all frequencies and their respective uvw coordinates. The memory footprint is estimated by

Mem STD = (N r .N ch .N pol + 3.N r ).N bytes , ( 34
)
with N bytes the number of bytes used for the precision of the data (single or double float), N r the number of visibilities for a single frequency channel and a single polarization, N ch the number of frequency channels, and N pol the number of polarizations. The memory footprint of the G2G method is mainly occupied by the dirty image as well as the compressed diag( A * ) calculated in part 3.2. This compressed array has five columns: three for the approximate uvw coordinates, one for the slice number, and one for the value of the corresponding element. Therefore, the memory footprint is estimated by

Mem G2G = (N S lices .N x .N y .N pol .N bytes ) + (5.n coo .N int32 ) (35)
with n coo is the length of the compressed diag( A * ) matrix, N S lices is the number of slices in the hypercube, N bytes is the number of bytes for the floating numbers, and N int32 is the number of bytes for an int32. The memory footprint remains constant regardless of K for the STD method. For the G2G method, the memory footprint varies with K. The smaller the K, the smaller the memory footprint because c α will be smaller. Furthermore, we compare two cases. The first is with a dirty image with as many slices as frequency channels. For each K, the memory footprint is smaller than the STD method. The main weight comes from the dirty image. In the case of the dirty image with one slice for all frequencies, the memory footprint decreases drastically. For each K, the main weight of the memory footprint is α. With K = 8, the memory footprint used with this method represents 4.4% of the memory footprint used with the STD method.

The memory footprint results follow the same logic as the compression ratio results in Tab. 2. The lower the sampling factor, the higher the compression. For K = 64, the required data volume is 27.8% of the initial data set. For K = 8, c α is only 3.4%. The difference between c α and c α is very small, so the theoretical compression gives reliable information of the gain more quickly. The G2G method obtains compression gain under M ′ < M. In practice, this is the case for most radio telescopes. The more compact the array, like VLA or NenuFar, the greater the compression gain. For less compact arrays, such as LOFAR or SKA, a compression gain is still expected, but less by an order of magnitude. On the other hand, the memory cost for noncompact networks working on large images can be dominated by saving the dirty image on several frequency bands. As a result, despite lossless data compression and reduced computational cost, the memory cost with some networks can be very high.

K

GPU performances

We will use three metrics to measure the CPU and GPU implementation performance. The first metric is the execution time in seconds. The second metric, define as giga Grid Point per Second (GPS), is the number of pixels updated or read per second during the execution of the kernel such that

GPS G2G (K) = #FourierPoint G2G .C 2 supp t G2G (K) . × 10 -9 , ( 36
)
with t G2G (K) is the GPU execution time for the G2G kernel regarding K, and #FourierPoint G2G = M ′ .N pol is the number of Fourier-point to be gridded and degridded. The last metric define as the compressed Million Fourier-point Throughput (MFT) is the number of data processed in the G2G kernel, such as

MFT G2G (K) = #FourierPoint G2G t G2G (K) × 10 -6 . (37)
Fig. 10 and Fig. 11 show the execution time of the G2G method for several convolution function sizes as a function of the oversampling factor for the CPU and the GPU implementations. In both cases, we can see the influence of the oversampling factor in the execution time because the number of elements to be processed depends on K. Moreover, the increase in execution time as a function of the convolution function is globally well proportional to its size. The exception is the case with C supp = 64 in GPU, where each thread takes care of several pixels and tends to accelerate the CUDA kernel. We can notice an order of magnitude difference in favor of the GPU implementation compared to the 64-core parallelized CPU implementation. This result confirms the interest in GPU implementation, even when we can use high-performance CPUs with many cores. Fig. 12 shows the evolution of GPS as a function of K for several convolution function sizes. For each C supp , the overall trend is an increase in performance as the oversampling factor is higher. However, the evolution remains choppy, and some K performs better than others. Surprisingly, this is not with K being a power of 2. The performances are separated into two blocks. The first one is with C supp = 16 and C supp = 32, whose performances are slightly lower than the others. The second with C supp = 8 and C supp = 64, whose last one confirms the good performances on execution time.

Fig. 13 illustrates the data throughput of the G2G method as a function of K for several convolution functions. As expected, the larger the convolution function, the lower the data throughput. Indeed, more CUDA threads are allocated to process a sample. As a result, the number of blocks executed in parallel is lower, thus reducing the throughput. On the other hand, the throughput remains relatively stable as a function of K.

MFT G2G

Figure 13: GPU kernel -MFT G2G Fourier-point throughput of G2G method regarding the oversampling factor K.

Comparison with STD

We also compare the GPU state-of-the-art gridder and degridder implementations following the Romein implementation with the G2G method. In this comparison, the STD method is built from the GPU implementations of these gridder and degridder. Moreover, we will only take into account the CUDA kernel times without taking into account the memory transfer times. We take this assumption for the case where the kernel times would cover the memory transfer costs. The metrics used to compare the performance of the G2G method with the STD method are the same as the previous one. The GPS is calculated as

GPS S T D = #FourierPoint S T D .C 2 supp t grid + t degrid × 10 -9 , (38)
with #FourierPoint S T D = N r .N ch .N pol is the number of visibilities to be gridded and degridded. The Fourier-point throughput is calculated as

MFT S T D = #FourierPoint S T D t grid + t degrid × 10 -6 . (39)
For the G2G and STD methods, we take an oversampling factor K = 8 for the convolution functions.

Fig. 14 compared the performances of the two methods. For C supp = 64, the STD method is slightly superior by 2%. For all other cases, the G2G method is better. For the case of C supp = 8, the MFT is 38% better for the G2G method than for the STD method, and goes down to 14% better for C supp = 32. It should be noted that these experiments greatly benefit the STD method since the raw visibilities used are taking great advantage of the GPU implementation to reduce the memory bandwidth. Finally, Figure 14: GPU kernel -GPS comparison between G2G and STD methods for

K = 8.
we compare the Fourier-point throughput processed per second for the different methods. The results are presented in Fig. 15. The performances for STD and G2G follow the previous result, and G2G has better performances for small convolution functions.

MFT

Figure 15: GPU kernel -Fourier-point throughput comparison between G2G, STD.

Roofline model

The roofline model [START_REF] Williams | Roofline: an insightful visual performance model for multicore architectures[END_REF] is a tool to visualize the possible limitations of an algorithm with respect to the maximum theoretical performance of the target architecture. The model is characterized by the device's peak performance in FLOPs and by the attainable memory bandwidth. This model aims to indicate the possible bottleneck of the specific algorithm. The roofline model for the GPU implementation of the G2G, the gridding and degridding algorithms are represented by Fig. 16 The three points are in the compute-bound area without being limited by the theoretical power of the GPU. The gridding reaches 22% of the theoretical performance, the degridding 11%, the STD 14%, and G2G between with 17%. The STD point was built from the gridding and degridding points. G2G improves the arithmetic intensity compared to STD. The result also shows a large margin of performance improvement before reaching the theoretical peak since the bandwidth does not limit the methods. In order to improve the performance, it is necessary either to increase the number of operations performed in the kernel time or to decrease the kernel time without changing the number of operations. The memory access time to access the convolution function coefficients and the lack of parallelism of the reduction operation are known to be bottlenecks. We can therefore expect much better results by improving these two points. One possible improvement is integrating the full Merry's implementation, Merry (2016), with coarsening threads.

Conclusion

We presented the Grid to Grid method used in radio interferometric imaging that can be seen as a generalization of the BDA. The gridding and degridding operators have a high computational cost. The G2G method shows the computational cost reduction by merging these two operators in a theoretical and experimental way. The main cost reduction factor is the oversampling factor used for the convolution function and the fineness of the approximation of true uv coordinates of the visibilities. This method makes no degradation of the reconstruction quality compared to classical gridding and degridding algorithms for the same oversampling factor K. Moreover, it can be used in addition to other averaging methods.

We have shown the impact of the oversampling factor on the precision and the error generated by it. We have also shown its impact on the memory footprint and execution time. Consequently, a compromise to be made on the oversampling factor to balance the computation acceleration while limiting the error. The G2G method showed a compression gain for all types of radio telescopes, especially if the array was compact. For noncompact arrays, we still expect a compression gain and a lower computational cost but potentially a higher memory cost due to dirty image storage. GPU implementation of the G2G methods is well suited by increasing the performance with a factor 10 regarding the CPU implementation. Finally, we showed that the GPU implementation of G2G was more efficient than the GPU implementation of a gridding/degridding w-projection algorithm with the same optimization level. The Fourier-point throughput is up to a 37% increase. On the other hand, the limitations come from the absence of direction-dependent effects processing as with A-projection. This algorithm requires much more memory access to load the convolution functions, which are time, frequency, baseline dependent, and not separable.

In future work, we plan to improve the CUDA kernel at several levels. The first is using texture to store convolution functions to improve accuracy. The second is to fully implementing Merry's GPU optimization is also a promising direction. In addition, we plan to implement the w-stacking method (Offringa and al., 2014), which allows using a low oversampling factor and thus works on larger images. Moreover, a higher level improvement is the possibility of using a hybrid version of w-projection and w-stacking. Finally, port the CUDA code to HIP, [START_REF] Kondratyuk | Gpu-accelerated molecular dynamics: State-of-art software performance and porting from nvidia cuda to amd hip[END_REF], in order to use the code on Nvidia or AMD platforms in a transparent way.

CPU experiments

The purpose of the experiments is to show additional results on the performances of the implementation of the G2G algorithm on CPU. The previous section only showed CPU computation time in order to be compared with GPU results. This section will provide CPU implementation details, parallelization scaling, and performance comparison with state-of-the-art gridder and degridder.

The experiments were done on a BullSequana XH2000 system based on Intel(R) Xeon(R) Platinum 8358 -32 cores dual-socket. The test case used is the same as the one used for the GPU version in the previous section.

CPU implementation

The CPU work distribution is more straightforward than the GPU implementation. In order to optimize the code as much as possible and to be able to reuse it easily according to the parameters of the image to be reconstructed, the code is based on a C++ template. The implementation of G2G follows the simplified code represented by Alg. 5. This algorithm is based on the pixel-to-pixel equation from Sec.3.2

g(b k ) = M ′ i ′ =1 C † b k -b i ′ a i ′ i,j C b i ′ -b ij g(b ij ).
(3.1) Eq. 3.1 give an explicit solution of the value of the pixel at coordinates bmb k ) = (u k , v k ). The algorithm build iteratively the grid g by looping on the M ′ non-zero coefficient of diag(A * ).

The double for loops Alg. 5 that read the pixels of the input grid and write to the output grid is responsible for the quadratic cost of the w-projection algorithm. The computational cost C w-grid = O(C 2 supp M ) depends on the size of the convolution kernel. Moreover, in the case of quad-polarization, numerous and expensive multiplications between two complex numbers are required. We have therefore optimized this step by taking advantage of the vectorization of these calculations with the AVX instructions available with recent Intel processors. Indeed, the compiler does not systematically vectorize such operations because of the depth of the code. Also, considering recent hardware, we have made these optimizations only with AVX instructions rather than SSE instructions. Indeed, AVX performances are, at worst, equivalent to SSE instructions. Indeed, performances of AVX instructions are, at worst, similar to SSE instructions and can get a significant speed-up factor when data is reused [START_REF] Jeong | [END_REF]]. Thus, an example of quad-polarization complex vector (c 1 , c 2 , c 3 , c 4 ) multiply by the same convolution function complex value p 1 is shows in Algo 4. Values c 1r and c 1i are, respectively, the real and imaginary parts of c 1 . The result of complex multiplication is

(c 1r + c 1i )(p 1r + p 1i ) = (c 1r p 1r -c 1i p 1i ) + (c 1r p 1i + c 1i p 1r ), (3.2)
where the left bracket of the result is the real part and the right bracket is the imaginary part. We can check how the Alg.4 works by looking at the first two elements

CPU Results

Fig. 3.4 shows the scaling of the G2G implementation parallelized on several cores with a convolution kernel and an oversampling factor K. The scaling is not linear when the parallelization is over 32 NUMA cores. Fig. 3.5 and Fig. 3.6 show the execution time, and the processed visibility throughput as a function of the oversampling factor of the G2G algorithm for several convolution kernel sizes parallelized on 64 cores and an oversampling factor K = 8. We find similar results to those on GPU. The execution time increases as a function of the oversampling factor while the processed visibility rate remains stable. However, the results are an order of magnitude lower than the GPU implementation. with each time parallelization on 64 cores. In each case, the performance of the G2G method is superior, up to 45% for large convolution kernels. The performance difference is more significant for the CPU than the GPU implementation. However, as for the execution time, the performance is an order of magnitude lower than the GPU implementations.

Conclusion

This chapter presents the Grid to Grid algorithm, merging the two operators gridding and degridding, to reduce the computational cost. We have shown that the order of the operators when computing the gradient of data fidelity term of a reconstruction algorithm in radio interferometric imaging allows for reducing computational costs. Moreover, we can extend this algorithmic optimization to other reconstruction problems, including MRI imaging. We have shown that this method reduces the computational cost and the memory footprint required for the calculations as a function of the oversampling factor K, which directs the accuracy of the interpolation. We also showed that this method could be applied to radio interferometric imaging by considering the correction of the coplanarity of the antenna array.

In a second step, we presented the implementation of this method on NVidia GPU and Intel CPU architectures. We have highlighted the effects of the oversampling factor on the reconstruction error and the execution time. The GPU implementation showed that our method was globally more efficient, with performances up to 37%, compared to the state-of-the-art. The CPU implementation, parallelized on multiple cores, also performed better than standard implementations. Overall, the GPU version is an order of magnitude better than the CPU version.

Algorithm 5: CPU G2G -pseudo code.

Data: (...,K, Igrid, Ogrid, nhit) In this chapter, we will look at optimizations for the coarse-grained parallelization of an imaging algorithm. We will propose a multi-core implementation of a DDFacet imaging framework based on the CLEAN algorithm to reduce the timeto-solution to reconstruct the sky image. First, we will present DDFacet in more detail and all the specificities of this framework. Then, we will see its algorithmic decomposition and propose an implementation on an HPC server with distributed memory, as well as a new multi-core implementation based on the same paradigm. Finally, we will show the results of this implementation. The results will focus on scalability, limitations, and perspectives for improving performance. 1

1 c xx ,

The Direction-Dependent Facets Imager (DDFacet)

In this section, we will introduce in more detail the specificities of the DDFacet imager. First, we will see the interest of the facets used in DDF to make observations on large FoV while correcting the perturbations. Then, we will study the specificities of the DDFacet implementation.

Non-coplanar and coplanar facet imaging

Faceting aims to approximate a wider FoV with many smaller FoV. [Cornwell 1992] proposed a method based on polyhedron-like faceting. In the polyhedron approach, each narrow-field facet is tangent to the celestial sphere at its own phase center (l φ , m φ ). It is usually categorized as non-coplanar faceting as each facet is set to its own tangent plane, and the transformation is done in the uvw space, not in the image domain.

To synthetise an image with a new phase tracking center from the original phase tracking centre (l 0 , m 0 , n 0 ), the shift theorem can be used. Let (l -l φ , m -m φ ) the coordinates in the reference frame of facet φ, and (l ∆ , m ∆ , n ∆ ) = (l φ -l 0 , m φm 0 , n φ -n 0 ). The simplfied RIME, Eq. 2.11, becomes

V pq = lm G p X(l -l φ , m -m φ )G H q e -2iπ(u(l -lφ)+v(m-mφ)+w(n-nφ)) dl dm n = lm G p X(l, m)G H q e -2iπ(ul +vm+w(n-n 0 )) dl dm n e 2iπ(ul ∆ +vm ∆ +wn ∆ )
(4.1) Thus, we can multiply each visibility by a complex exponential. From Eq. 2.15, 1 The contents of this chapter have been partially published in [Monnier 2022]. v i e 2iπ(u i (l-lφ)+v i (l-lφ)) e 2iπ(ul ∆ +vm ∆ +wn ∆ ) .

(4.2)

However, the intensity is still measured relative to the uvw coordinates of the original phase center. [Taylor 1999] has shown that in this case, the facets created have different phase centers, but their tangent plane is the same as the original image plane phase center, as illustrated by Fig. 4.1. In this figure, the first facet is tangent to the observation's phase center, and ε 1 represents the distance between the image plane and the celestial sphere. The second facet is centered on θ i , the angle from the phase tracking center, parallel to the original image plane, and ε 2 is the distance between the image plane and the celestial sphere.. The new image plane remains parallel to the original image plane. Image adapted from Rick Perley [Taylor 1999] In order to make facets tangent to their phase center, new (u ′ , v ′ , w ′ ) coordinates must be computed from the (u, v, w) coordinates using rotation matrices. These new coordinates are calculated after the visibilities are phase shifted such that

Celestial sphere

   u ′ v ′ w ′    = R(α φ , δ φ )R T (α, δ)    u v w    (4.3)
where the angular coordinates α and δ are the right ascension and the declination of a point in the sky. In effect, ensuring the facets are tangent to the new phase centers ensures that the height difference between the sky and the projected facet remains comparable between corresponding pixels [Hugo 2016]. Polyhedron facets, or non-coplanar facets, are illustrated in Fig. 4.2.

One of the biggest issues of polyhedron faceted imaging is the difficulties regarding the deconvolution's minor cycles. Indeed, as the visibilities' coordinates of each facet are rotated, these rotations are preserved after the Fourier transform. The PSF of each facet is also rotated. Thus, each facet has its own PSF, and the subtraction phase in CLEAN minor cycles requires careful consideration.

One way to solve this problem is coplanar faceting by correcting the noncoplanarity during gridding within each facet. In order to take into account the Celestial sphere phase error δn φ , coupled with the w-term, Kogan [Kogan 2009] showed that the phase error could be approximated by a first-order Taylor expansion such that

δn φ = 1 1 -l 2 φ -m 2 φ (l φ (l -l φ ) + m φ (m -m φ )) (4.4)
This leads to new coordinates u ′ and v ′ such as

u ′ = u -w l φ 1 -l 2 φ -m 2 φ v ′ = v -w m φ 1 -l 2 φ -m 2 φ
Tasse [Tasse 2018] showed a more accurate approximation of δn φ by fitting a low-order 2D polynomial. These approximations allow pre-computing a convolution kernel per facet that will be applied in the same way as the w-projection.

DDFacet specificities and features

Astrophysics fluxes are not uniformly distributed over the observed sky. The DDE calibration must therefore cluster groups of sources to determine the Jones matrices in all these directions. DDFacet makes images from images irregularly tesselated. The tesselation to cluster the sources is done from a Voronoi tesselation-like algorithm. Thus, each facet is mapped to a polygon resulting from the tessellation, and then masked when the facets are combined to make the dirty image. Experience shows that we need to split the sky model into a few tens of directions to be able to describe the spatial variation of the direction-dependent Jones matrices [van Haarlem 2013]. Moreover, it is common to compute a PSF per facet such as the PSF for facet i is

y P SF φ i = τ φ i 0 1 , (4.5)
where τ φ i is the local convolution approximation for facet φ i , and 0 1 is a vector of the size of the facet fill with O but the center of the facet. There is also the possibility of reducing the computational cost in two ways. The first one is averaging visibilities with the BDA algorithm to reduce the number of data while minimizing the decorrelation [Atemkeng 2016a, Wijnholds 2018b]. The second method is to use sparsification, based on the Compress Sensing theory. The idea behind this is that in the first major cycles, in a CLEAN-type algorithm, the deconvolution tends to be highly sparse because there are only a few extremely bright sources, and they are usually far from each other. The sparsification allows taking randomly only a small percentage of visibilities for gridding during the first major cycles. The information is sufficient to create the model of the brightest sources. The percentage of visibilities selected increases with the next major cycles. DDFacet incorporates natively a few deconvolution algorithms and compensates for the DDE. For instance, one can use the standard Multi-Scale Multi-Frequency Synthesis(MSMFS) CLEAN algorithm [Rau 2011]. As it is well-known that the CLEAN-like algorithms are not robust in the deconvolution of extended emission, the SubSpace Deconvolution algorithm (SSD) is implemented. As described in [Tasse 2018], we can also find an extension of SSD using a genetic algorithm (SS-DGA). For the next of this chapter, every mention of the deconvolution algorithm will reference the MSMFS algorithm implementation described in [Tasse 2018].

McMn parallelization

This section describes the parallelization paradigms implemented in DDFacet. First, we will introduce the simplified sequential algorithm to introduce notations. Then, we will see the multi-core parallelization with a modified implementation compared to the native version to simplify the implementation of the multi-node parallelization. Finally, we will propose a multi-node parallelization to reduce the latency to generate a reconstructed sky image.

Sequential algorithm

The visibilities for a specific observation and a specific frequency band, consisting of several frequency channels, are stored in the "Measurement Set" (MS) format. For observations on large frequency bands, the data of the same observation are therefore separated by frequency band in several MS files. The structure of DDFacet's image reconstruction method is described by Alg. 6. This algorithm follows the CLEAN Cotton & Schwab method, presented by Alg. 3. The K th major cycles correspond to the update of the residual image, such as

δy (k) = F † S † (v -SF x (k) ), (4.6)
and the gridding and degridding steps incorporate the direction-dependent corrections. This sequential algorithm is also illustrated in Fig. 6. The color identifies the data domain. For instance, the residual image δy (k) and the model image x in blue are in the image domain. The red color represents the discrete Fourier domain, like the Fourier transform of the model g or the grid obtained after gridding g. Finally, the green color represents the data with continuous coordinates like the raw visibilities v, the model visibilities v, and the residual visibilities δv.

The main difference between Eq. 4.6 and what is done in DDFacet is in the implementation. The image is divided into I facets, so the FFT, gridding, and degridding operators are performed independently between all facets of the same image cube through for loops. Moreover, the residual image is built iteratively by processing the visibility vector v M S,j for each of the J th MS file separately and independently. If the residual image is a hyperspectral cube, then the data for each frequency is accumulated on different slices of the cube. Otherwise, each gridded data will accumulate on the same slice before going into the image domain. The deconvolution step is deliberately left generic because several algorithms can be integrated into this function. The iterations of this deconvolution function are called Minor Cycles. In our case, the chosen algorithm is the MSMFS algorithm.

Multi-core parallelization

In order to facilitate the implementation of distributed memory parallelization, a new implementation of multi-core parallelization has been made using Concurrent-Futures2 . This new implementation, following the same parallelization paradigm as the native version [Tasse 2018], is illustrated in Fig 4 .4.

Concurrent-Futures provides a high-level interface for asynchronously executing callable methods. The interest in using such a library is to ensure better maintainability in the code. Moreover, many libraries for parallel computing, like Dask [Rocklin 2015], are based on Concurrent-Futures. This implementation also paves the way to the future if it is necessary to change technology.

In the case of DDFacet, the image is divided into facets whose calculations are independent of each other. Thus each loop implying facets, as gridding, degridding, and the computation of the residual image δy of Alg. 6 is parallelized. The parallelization framework is based on an asynchronous behavior between the computational and I/O phases. The main process of the imager manages a bunch of workers, where each worker is a process, depending on the number of cores available or set by the user. The main process then dispatches the different jobs in a dedicated queue. A compute queue is dedicated to the computational jobs like gridding, degridding, and FFT for each facet φ, and an I/O queue for the I/O relative jobs.

Multi-node parallelization with data distribution

A node is an individual computer consisting of one or more processors which contains multiple cores sharing the same memory. Thus, the naive multi-node parallelization paradigm used to increase the image production rate is to distribute tional load. Each node computes its metadata regarding its sub-dataset, such as BDA computation, weight computation, and DDE Jones matrix. The metadata is then gathered for all nodes to know the reconstructed hyperspectral cube parameters. After that, the computational core of each node is identical, using the multi-core parallelization described in Sec 4.2.2. The FFT, degridding, and gridding operators will be applied depending on the dataset. After the gridding of all MS for each node, the subgrids g (k)

1..N are transferred and reduced to the master node such as

g (k) = N i=1 g (k) i (4.7)
Only the master node is in charge of moving into the spatial domain and applying

Implementation

The software implementation has been done in Python using the Concurrent-Futures interface of mpi4py3 based on Message Passing Interface (MPI). MPI is a standard message passing interface for applications and libraries running on concurrent computers with a logically distributed memory [Walker 1996]. MPI allows the introduction of a communicator so that a group of processes can make point-to-point or collective communication. Very often, MPI processes, in the same communication, execute the same job with the specificity of its rank number. We call rank 0 the process considered as the master and rank N the Nth MPI process in a given communication. Communication can be more complex than a simple point-to-point data transfer. The most useful collective communications are :

Reduce : takes an array of input elements on each process and returns an array of output elements to the master process. The reduction operation can be the sum of the elements, the minimum element, the maximum element, etc.

Gather : takes elements from many processes and gathers them into one or multiple processes.

Broadcast : one process sends the same data to all processes in a communicator.

In our case, the Concurrent-Futures interface of mpi4py allows spawning MPI processes from rank 0 in the same way as Concurrent-Futures for multi-core parallelization. Thus, this multi-node implementation makes two levels of DDFacet parallelism that are managed consistently between the two. The first native one is based on facets parallelization and communication is made using shared memory. The second one is based on Measurement Set parallelization, which allows for a distributed memory system. Moreover, MPI allows better use of the computational resources within a node. Indeed, it is possible to spawn several MPI processes on a single computing node. Thus, when the processor has a large number of cores or when we use few facets to reconstruct an image, such as

I f acets N cores < 1, (4.8)
it results in an under-utilization of the computational resources. Therefore, we can take advantage of this two-level parallelization by spawning several MPI processes per compute node. In this case, each MPI processes share the computation resources equally and is optimal when

N mpi/node ≥ N cores I f acets . ( 4

Gridding and degridding speed-up α

We are interested in the scalability factor of the main part distributed over several MPI processes: gridding and degridding. For this, we define the scalability factor α p (n) as is the time of the master for p MPI processes per node distributed on 1 node. It is expected that α p (n) should scale as n the number of nodes. The scalability factor for the case of 1 MPI process per node is illustrated in Fig. 4.9. We can see that the scalability factor follows the theoretical curve entirely. This result shows us the efficiency of the distribution of the calculations on several nodes and that the related overheads are negligible. Fig. 4.10 and Fig. 4.11 show the scalability curves for 2 MPI processes per node and 4 MPI processes per node. The scalability factor also follows the theoretical curve in the case of 4 MPI processes per node. In the case of 2 MPI per node, the scalability curve is slightly worse but still follows the theoretical curve very closely. The speedup of each case with respect to the reference time on 1 node with 1 MPI process, represented by Fig. 4.12, is defined as

α p (n) = t GDG p,1 t GDG p,n , ( 4
α ref p (n) = t GDG 1,1 t GDG p,n . (4.12)
On one hand, each case still shows a good linearity. On the other hand, the curves do not show a perfect acceleration from 2 MPI per Proc. It might be explained by the fact that the computational resources are shared between the different MPI processes.

Global speed-up b eta

Fig. 4.13 shows the total speedup using the reference of 1 MPI process per node. The speedup is calculated as

β p (n) = t T 1,1 t T p,n , (4.13)
where t T p,n is the total execution time of DDFacet, with p the number of MPI processes per node and n is the number of node. At any point, the speedup is significant, reaching β 1 (16) = 5.7. However, the speedup reaches an asymptote at 5.3 in all cases. As a result, the speedup slowly increases when parallelization is done on a high number of nodes. Besides, we see the interest in using several MPI processes per node for our test case. Indeed, with 4 MPI processes per node, we can see a decrease in the computational time. For example, we have β 4 (1) = 2.01, showing an improvement of the multi-core parallelization thanks to the implementation of the distributed version. The total execution time and the profiling of the different steps of the algorithm are illustrated by Fig 4 .14. For each measure, the profiling was done on rank 0, i.e., the master MPI process that manages the deconvolution. For the case of 1 MPI process on one node, the gridding and degridding steps represent 75%, which confirms the result of [Tasse 2018]. We can also see that the computation time is divided into four main blocks: Gridding and degridding, Initiatization, Deconvolution and Others. This last one is an accumulation of all the short unmeasured times, some latencies of memory access to the shared file system, and saving in the cache. This result highlights the need to distribute this core computation to reduce the total execution time. As expected, the total execution time tends to decrease with distribution on several nodes, demonstrating the interest of this implementation to reduce the time-to-solution to produce an image.

Following the parallelization diagram in Fig. 4.7, we can see that the distributed parallelization steps concern gridding, degridding, FFT, and the initialization step. The same corresponding timings significantly decrease when parallelized on several nodes and several MPI processes per node. This result confirm the scalability results show in the previous subsection. We see that on a large number of nodes, the ratio of the execution time of the gridding and degridding only represents a few percent of the total execution time. As a result, the execution time is no longer limited by the gridding but rather by the other incompressible times.

As soon as MPI is used, the communication and reduction times "MPI" appear, as well as the synchronization times with other processes "Wait". These two times remain very low but tend to increase with the number of MPI processes. This increase seems logical since MPI operations are done with a more significant number of processes and therefore take more time. Furthermore, the synchronization can vary greatly depending on the balance of MS among MPI processes. Finally, the deconvolution time remains relatively stable for each case. The only minor variations come from some IO latency. Otherwise, the number of cycles and the computational load is identical for each case. This incompressible time becomes one of the limits for reducing the computation time when parallelizing on a large number of nodes. (4.14) and the total speedup as

Analyses

γ G (n) = t GDG 1,n t GDG 4,n ,
γ T (n) = t T 1,n t T 4,n , ( 4.15)
The gridding speedup remains constant with a factor γ around 2, which demonstrates an improvement of the multi-core parallelization. However, the total execution time speedup decreases. We can therefore conclude that the main limiting factor for the reduction of the execution time is the fixed cost times and the overheads. These times, as in Others, contain pre-process time or data movement which are more difficult to reduce without a deeper analyze which is not the aim of this work.

The limitations of the total speedup of Fig. 4.13 can be analyzed using Amdhal's law [Amdahl 1967], whose formula for the theoretical acceleration of a spot is given by where S latency is the theoretical speedup, p is the proportion of executing time that the parallelized part occupied initially, and s is the number of nodes used.

S latency = 1 (1 -p) + p s , ( 4
In our case, only the gridding and initialization steps are relevant to studying the theoretical asymptote. The results give us a theoretical curve limited to a value of S latency = 9, 09, and if we take a hypothetical parallelization on 24 nodes, the theoretical result is S latency = 6, 79. Compared with our experimental results, our implementation has added overhead to the execution, not allowing us to reach the theoretical speedup. The main known overheads are the communication time and the MPI waiting time.

The parallelization efficiency η, illustrated by Fig. 4.16 is defined as

η p (n) = β p (n) n .
(4.17)

For all three cases, the efficiency decreases until it reaches an efficiency of 21% for η 1 (24). On a smaller number of nodes, the efficiency for 4 MPI processes per node and 2 MPI processes per node is higher than 1. We can therefore look at the balance between the speedup and the hardware resources. The maximum speedup is obtained for β 1 (16) = 5.7, which has an efficiency of η 1 (16) = 35%. However, η 4 (4) has an efficiency of 1 for a speedup of β 4 (4) = 4.1. The speedup is thus slightly lower using 4 times fewer resources. The efficiency curve, therefore, provides a good metric for choosing a balance between speedup and hardware resources to ensure optimal parallelization.

Memory footprint

Finally, we can also observe the memory consumption during gridding and degridding, as indicated by the Tab 4.1. This maximum memory consumption on the master node shows an increase in the memory footprint when increasing the number of MPI processes on the same node. This increase is justified by reading several MS per node when several MPI processes coexist on the same node.

Conclusion

This chapter presents McMn DDFacet, a wide-band wide-field spectral deconvolution framework on a distributed memory HPC system. We have presented a two-level parallelization paradigm with facet-based parallelization for the multi-core shared memory aspect and MS-based parallelization for the multi-node distributedmemory aspect. The multi-core parallelization is implemented with the concurrentfuture library to ensure asynchronous parallelization. The communication is done with dictionaries and arrays stored in shared memory. The multi-node parallelization is performed on MPI calls of the mpi4py Concurrent Future library. This implementation allowed us to distribute the most expensive computational blocks, such as gridding and degridding, ensuring scalability. Using a LOFAR dataset, the total time to generate a reconstructed image was reduced by a factor up to 5.7. Moreover, we observed excellent gridding and degridding parallelization scalability on multiple nodes. We have also shown that our implementation can improve the performance of multi-core parallelization by using multiple MPI processes on a node in the case of a large number of cores compared to the number of facets. This implementation also shows the interest of having compute nodes with large processors to be able to optimize the parallelization on several cores in a more efficient way. Finally, the overall performance presented in the experiment is limited by the size of the test case used. A much larger test case would have given better performances in terms of speed-up. The perspective for the future is to continue working on ways to reduce the execution time. As this work was focused on gridding and degridding parallelization, the next track is to apply our multi-node framework to distribute the deconvolution computation on several nodes. Moreover, to optimize the computational load between several nodes because of the compression applied by the BDA algorithm. Finally, to implement optimization for gridding and degridding on acceleration platforms for DDFacet. Chapter 5

General conclusion

Summary of the contributions

This manuscript has detailed the research effort conducted with the objective of reducing the computation time of imaging algorithms in radio interferometry based on both algorithmic methods and hardware implementation. This work was motivated by the approach of the construction of a new generation of radio telescopes. These new instruments will bring extreme computational constraints that must be understood to design the computation servers. It is, therefore, crucial to tackle the problem as a whole. Consequently, we propose solutions to reduce the computation time at all levels of the problem. For the algorithmic part, we propose a method to optimize expensive computational blocks to reduce the computational cost of imaging algorithms. We also propose several parallelization strategies for multiscale hardware implementations. First, a fine-grained parallelization to optimize the implementation of a computationally intensive kernel on dedicated acceleration platforms. Then, a coarse-grained parallelization to optimize image reconstruction algorithms computation on high-performance computing servers.

In Chapter 3, we have presented a decomposition of the interpolation operators, gridding and degridding, into sub-operators allowing a visual and mathematical representation of their functioning. In the first step, from these decompositions, we proposed a method named Grid to Grid (G2G) allowing to merge these two operators to decrease the memory footprint and the computational cost. In our case, we have applied it in the context of radio astronomy to correct the effects associated with the antennas' non-coplanarity. With a practical case of the VLA telescope, we have shown that we could predict the reduction of the computational cost and the memory footprint required for the computations using the G2G method. We also showed that these gains must be balanced between gains and oversampling factor error. The higher the gains, the higher the error. Finally, we have confirmed that the interpolation operators were more efficient on GPU than on CPU with an order of magnitude difference in performance and execution time, whatever the size of the convolution kernel. Moreover, we have shown that the GPU implementation of 82 Chapter 5: General conclusion the G2G method matches state-of-the-art gridders with a gain of almost 40%, as well as more recent methods like IDG.

When Chapter 3 showed results based on an algorithmic method with finegrained parallelization on GPU, Chapter 4 focused on coarse-grained parallelization of an imaging framework on an HPC server. We used the DDFacet imager, which allows us to correct Direction-Dependent Effects (DDE) by splitting the sky image into several facets. In order to have a homogeneous computation distribution management, we proposed a new implementation of multi-core parallelization, based on the same parallelization paradigm as the native version, using the concurrent-future technology to launch concurrent instructions. This implementation was used as a basis to distribute computations on several processes. We have chosen the MPI technology to manage the computation parallelization on coarser-grain in order to take advantage of the communications that can be done inter-node or intra-node in a completely transparent way. The multi-core parallelization is based on the independence of the calculations between the facets of an image. MPI parallelization is based on the independence of the calculations between the different wavelengths of the visibilities. Our proposal allowed us to distribute the most time-consuming computations, gridding and degridding, from 80% of the total execution time to only 13% when distributed over several AMD EPYC 7742 64-core processors. The implementation showed very good scalability for distributed computational blocks and improved multi-core parallelization thanks to the MPI implementation.

Perspectives

Several areas of improvement emerge from the work explored in this manuscript. Many improvement and optimization axes are possible for the GPU implementation of the Grid to Grid method presented in chapter 3. These concern two main points: better memory management and better low-level optimization. Indeed, one slowing down factor is the reading of the convolution kernels, which requires many memory accesses. Using more efficiently all the available memory areas, such as the texture or the constant memory, would reduce latency and the error related to the approximation. A low-level code analysis also optimizes the number of registers used to use the hardware better and decreases the time spent in the kernel. The w-stacking or hybrid methods (w-stacking and w-projection) are also interesting tracks that would lend themselves well to the G2G method. This method also opens the door to all image reconstruction algorithms that are more sophisticated but slower than the CLEAN algorithm. Reducing the computation time to compute the gradient of the data fidelity term, mandatory in many of these algorithms, could make these methods more achievable on a large scale. Moreover, using a preconditioner in gradient descent algorithms, such as conjugate gradient, could further accelerate the computation by decreasing the number of iterations required for reconstruction. In addition, this proposition has the advantage of being applicable to all Fourier synthesis problems, such as MRI imaging problems. Thus, implementing G2G can be considered for all reconstruction methods requiring this type of interpolation.

Chapter 4 showed the limits of our implementation with the test case we used. It would be interesting to do additional tests with different configurations to go further. The most interesting configuration would be with an even bigger test case, an authentic simulation of the data acquired with SKA, for example. In this way, we could show the excellent scalability of the whole framework in this type of test case. Moreover, the possibility of distributing the computations has opened the door to a lot of improvement in parallelization. Considering the large images reconstructed with SKA, the distribution of the computations to calculate a single grid is possible.

On the other hand, few deconvolution algorithms allow their computation to be distributed on several nodes. However, some, like SSD, can. In order to additionally reduce the computation time, the opportunity of optimizing this type of algorithm with a distribution of calculations is an essential step. Moreover, the developed parallelization framework can be reused for other algorithms. The deconvolution step is a mandatory step and also very expensive in terms of computation time. Fine and coarse grain optimizations on implementing such algorithms likely lead to significant gains in computation time.

Résumé en Français

L'objectif d'un radio interféromètre est de simuler un très grand télescope à partir d'un réseau d'antennes espacé spécialement. L'imagerie est une technique permettant de reconstruire une image du ciel observé à partir des signaux reçus par les différentes antennes du réseau. Les données utilisées pour la reconstruction se font à partir d'une observation indirecte. Le ciel observé est dans le domaine spatial et les données utilisés, les visibilités, sont situé dans le domaine de Fourier. Le problème de reconstruction est un problème "mal-posé" car les mesures ne couvrent par l'entièreté du plan de Fourier, impliquant donc une infinité de solutions, et qu'elles sont également corrompue par des effets liés à la propagation des signaux dans l'atmosphère terrestre. À cause de la grande dimension du problème de reconstruction de l'image du ciel, des algorithmes itératif, introduisant de l'information a priori sur le ciel, sont nécessaires pour la reconstruction. De nombreuses méthodes différentes existent comme les méthodes baysienne, le compress-sensing, ou la méthode historique CLEAN. Chacune d'entre elles a ses spécificités quant à la qualité ou la vitesse de reconstruction, mais ont presque toutes un point commun : la nécessité d'interpoler les visibilités sur une grille uniforme afin d'utiliser des algorithmes de transformée de Fourier rapide.

Dans le modèle « backward », l'interpolation, nommée gridding, consiste à étaler la visibilité, dont la vraie coordonnée est continue, sur une grille uniforme en utilisant un noyau de convolution. Dans le modèle « forward », l'interpolation, nommée degridding, est l'opération adjointe et regroupe l'information sur une large zone centrée sur la position de la visibilité. Afin de minimiser l'erreur lors du gridding et degridding, ces interpolations convolutives sont faites en utilisant un noyau de convolution spécifiquement construit pour l'opération.

Les radiotélescopes ont pour caractéristiques communes de générer un débit de données extrêmement important rendant très lourd le traitement et la conservation des visibilités. Avec la nouvelle génération de ces interféromètres, le débit de visibilités à traiter augmente d'un ou plusieurs ordres de magnitude. Le cas du radiotélescope SKA est caractéristique, le débit est de plusieurs TB/s de visibilités à traiter pour reconstruire des images de plus de 30.0000 x 30.000 pixels. En plus de la correction des données extrêmement coûteuse en temps de calcul, la recon-88 Chapter 5: General conclusion struction des images est un véritable challenge. Le coût calculatoire des opérateurs d'interpolation, le gridding et degridding, des algorithmes de reconstruction est très élevé et devient un goulot d'étranglement à cause de l'aspect itératif des algorithmes. De ce fait, la reconstruction d'une image se compte en heures ou en jours de calculs intensifs sur des supercalculateurs.

Les motivations de cette thèse portent sur les différentes manières de diminuer le temps de calcul des méthodes d'imagerie en se concentrant sur deux aspects. Le premier est l'aspect algorithmique afin de diminuer le coût calculatoire nécessaire par les opérateurs impliqués dans les méthodes de reconstruction du ciel. Le second aspect est l'implémentation hardware venant avec une parallélisation à grain fin sur des accélérateurs multi-coeurs et many-coeurs et une parallélisation à gros-grain sur serveur HPC.

Nous présentons une méthode de réduction du coût calculatoire des opérateurs de gridding et degridding en les fusionnant en un opérateur unique. Cette méthode, nommée Grid to Grid (G2G), s'appuie sur la succession des deux opérateurs, ainsi que de l'approximation des coordonnées des visibilités sur la grille de Fourier. Cette méthode peut être adaptée pour corriger les effets liéés à la coplanarité du réseau d'antennes. En complément, nous proposons une implémentation CPU et GPU de cette méthode. Avec nos expériences, nous montrons que G2G permettait de réduire grandement le coût calculatoire ainsi que l'empreinte mémoire sans pénaliser la qualité de la reconstruction. Le facteur de surréchantillonage utilisé pour réduire l'erreur d'approximation des coordonnées sert de balance entre la réduction du coût calculatoire et la précision de l'interpolation pour diminuer l'erreur. De plus, nous montrons que les implémentations CPU et GPU de la méthode G2G sont plus efficaces que l'état de l'art en termes de débit de données traiter et en temps de calcul.

Nous présentons également une distribution multi-coeur multi-node sur serveur HPC du framework d'imagerie DDFacet. Cet imager permet de reconstruire des images du ciel en prenant en compte les effets de corruptions des données liés à la direction d'observation. Afin de les corriger, l'image est divisée en plusieurs images plus petites, nommées facettes, pour simuler l'observation d'un grand champ de vue. La distribution des calculs est divisée en plusieurs niveaux. Le premier niveau est une parallélisation multi-coeur pour systèmes à mémoire partagée basé sur l'indépendance des calculs entre les facettes. Le second est une parallélisation multi-noeud pour système à mémoire distribué basé sur l'indépendance des calculs du gridding et degridding entre les différentes fréquences d'observation. Ce double niveau de parallélisation a pour effet de réduire grandement le temps d'exécution. Nous montrons que sur un grand nombre de noeuds, le temps d'exécution est limité par les parties non parallélisables du framework. De plus, l'accélération n'étant pas linéaire, nous montrons qu'il est possible de choisir un optimum entre l'accélération et les ressources de calculs utilisés.

En résumé, l'ensemble de nos travaux au cours de cette thèse nous a permis en premier lieu de poser un état des lieux faisant le pont entre les communautés de l'astrophysique, du traitement du signal et optimisation hardware. Ensuite,

  ]. . . . . . . . . . . . . . . . . . . . . . . . 2.6 (a) : Relative position of the 27 antennas of the VLA radiotelescope. (b) uv track for a snapshot observation. (c) uv track for a 4-hour observation. (d) uv track for a 10-hour observation. . . . . . . . . . 2.7 Single point source image quality for a VLA observation of (a) 10 min, (b) 1 hour, (c) 2 hours, and (d) 8 hours. . . . . . . . . . . . . .2.8 Several calibration regimes, from [Wijnholds 2010]. (Left) narrow FoV, same ionospheric and tropospheric phases. (Right) large FoV, several calibrator sources. . . . . . . . . . . . . . . . . . . . . . . . . 2.9 The non-diagonality of S † S. . . . . . . . . . . . . . . . . . . . . . . . 2.10 Data rate of the different pipeline stages for the first phase of SKA. Image Domain Gridding illustration. First, a group of visibilities is shifted to low frequencies. Then, a DFT is performed on the shifted visibilities to make a low-resolution image. A Fast Fourier Transform is applied to the low-resolution image to make gridded visibilities. Finally, gridded visibilities are shifted back and added to their original place. . . . . . . . . . . . . . . . . . . . . . . 32 3.1 Evolution of CPU clock rate and power consumption for eight generations of Intel x86 microprocessors. From [Patterson 2016]. . . . . 36 3.2 Distribution of chip resources for a CPU vs a GPU. Figure from [Nvidia 2007]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3 GPU memory hierarchy. Each thread has private local memory, each block of thread has shared memory available for all the threads in the block, All threads have acces to the same global memory. From [Nvidia 2007]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4 CPU -G2G scaling. . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.5 CPU -G2G computation time regarding the oversampling factor K. 54 3.6 CPU -G2G visibilities throughput regarding the oversampling factor K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.7 CPU -Giga grid point read-write comparaison between G2G and STD methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.1 If the visibilities are phase-shifted to the facet center without baseline rotation, the image plane is phase shifted to the location of the new phase center. The new image plane remains parallel to the original image plane. Image adapted from Rick Perley [Taylor 1999] . . . . . 61 4.2 Polyhedron imaging aims to approximate the celestial sphere by a polyhedron. Image adapted from Rick Perley [Taylor 1999]. . . . . . 62 4.3 DDFacet sequantial imaging pipeline. . . . . . . . . . . . . . . . . . . 66 4.4 Multi-core parallelization . . . . . . . . . . . . . . . . . . . . . . . . 67 4.5 Concurrent-futures. As the length of the process pool is two, two jobs at the same time can be processed. The other jobs are waiting in the dedicated queue. If the job return a result, it is stored in the result queue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.6 Existing DDFacet multi-node parallelization. Each node generates an independent image using different dataset of a specific observation. 68 4.7 Multi-core and multi-node parallelization based on independent MS files using distributed memory system. . . . . . . . . . . . . . . . . . 68 4.8 Left : HPC system based on dual-socket node and common Shared file System. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.9 Scalability α 1 (n) of the gridding/degridding using one MPI process per node. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.10 Scalability α 2 (n) of the gridding/degridding using 2 MPI process per node. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 α 4 (n) of the gridding/degridding using 4 MPI processes per node. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 speedup α ref (n) of the gridding/degridding regarding the number of nodes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Total execution time SpeedUp regarding the number of nodes used. . 4.14 Profiling of the master MPI process for a complete execution regarding the number of MPI processes per node (1, 2, and 4), and the number of nodes used for the parallel execution. . . . . . . . . . . . 4.15 Speedup of gridding/degridding and total execution time between 4 MPI processes and 1 MPI process regarding the number of nodes. . 4.16 Parallelization efficiency regarding the number of nodes. . . . . . . . 4.17 (Left) Dirty Image. (Right) Restored image after 20.000 iterations of MSMFS Clean. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2

 2 Figure 2.1: (left) Ruber's dish radio telescope, https://en.wikipedia.org/wiki/ Radio_telescope. (right) Arecibo dish radio telescope, https://en.wikipedia. org/wiki/Arecibo_Telescope.

Figure 2

 2 Figure 2.2: Stations of the LOFAR radiotelescope.

Figure 2

 2 Figure 2.3: Image quality of SKA-mid (left) and the best current facility operating at the same frequency range of the JVLA (right). From [ska 2021].

Figure 2

 2 Figure 2.6: (a) : Relative position of the 27 antennas of the VLA radiotelescope. (b) uv track for a snapshot observation. (c) uv track for a 4-hour observation. (d) uv track for a 10-hour observation.

Figure 2 . 7 :FFigure 2

 272 Figure 2.7: Single point source image quality for a VLA observation of (a) 10 min, (b) 1 hour, (c) 2 hours, and (d) 8 hours.

Figure 2

 2 Figure 2.9: The non-diagonality of S † S.

  the residual image δy to the dirty image y. */ 3 δy ← y; /* return value of the brightest pixel i max and its pixel coordinates (k, l). */ 4 i max , k, l ← f indP eak(δy); /* Loop until the number of iteration n reach a limit set by the user or the brightest pixel i max reach a threshold (usually noise level value). */ 5 while i max > threshold and n < minor max do /* Add ratio γ of i max to the sky model at the position of the brightest pixel (k, l). l) ← x(k, l) + γ × i max ; /* Update δy by subtracting a ratio of the PSF y P SF centered at the brightest pixel location (k, l). */ 7 δy ← δy -γ × y P SF (k,l) ; 8 i max , k, l ← f indP eak(δy); 9 n ← n + 1; end /* Final Sky model is the convolution of x with clean PSF added to δy.

  the residual image δy to the dirty image y. */ δy ← y; /* Init the PSF y P SF as the center of the PSF until the first sidelobe. */ y P SF ← small(y P SF ); /* return value of the brightest pixel i max and its pixel coordinates (k, l). */ i max , k, l ← f indP eak(δy); /* Loop until the number of major iteration m reach a limit set by the user. */ while m < major max ; do while i max > threshold and n < minor max do /* Add ratio γ of i max to the sky model at the position of the brightest pixel (k, l). */ x(k, l) ← x(k, l) + γ × i max ; /* Update δy by subtracting a ratio of the PSF y P SF centered at the brightest pixel location (k, l). */ 10 δy ← δy -γ × y P SF (k,l) ; 11 i max , k, l ← f indP eak(δy); 12 n ← n + 1; end /* Evaluate a new δy by subtracting the convolution of x and y P SF to y. */ 13 δy ← yy P SF * x; 14 m ← m + 1; end /* Final Sky model is the convolution of x with clean PSF added to δy. the residual image δy to the dirty image y. the PSF y P SF as the center of the PSF until the first sidelobe. */ y P SF ← small(y P SF ); /* return value of the brightest pixel i max and its pixel coordinates (k, l). */ i max , x, y ← f indP eak(δy); /* Loop until the number of major iteration m reach a limit set by the user. */ while k < major max do while i max > threshold and n < minor max do /* Add ratio γ of i max to the sky model at the position of the brightest pixel (k, l). */ x(x, y) ← x(x, y) + γ × i max ; /* Update δy by subtracting a ratio of the PSF y P SF centered at the brightest pixel location (k, l). */ δy ← δy -γ × y P SF (x,y) ; i max , x, y ← f indP eak(δy); n ← n + 1; end /* Generate model visibilities v by Fourier transform and degrid δy by grid and inverse Fourier transform residual visibilities δv. */ δy ← f ourierGrid(δv); end /* Final Sky model is the convolution of x with clean PSF added to δy. */ x ← x * clean P SF + δy;

  Figure 2.10: Data rate of the different pipeline stages for the first phase of SKA.

  -core and many-core architectures . . . . . . . . . . . . 36 3.1.1 Multi-core CPU architecture . . . . . . . . . . . . . . . . . . 36 3.1.2 Many-core GPU architecture . . . . . . . . . . . . . . . . . . 37 3.2 Fast Grid to Grid (G2G) interpolation . . . . . . . . . . . . 38 3.3 CPU experiments . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.1 CPU implementation . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.2 CPU Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

Figure 3 . 1 :

 31 Figure 3.1: Evolution of CPU clock rate and power consumption for eight generations of Intel x86 microprocessors. From [Patterson 2016].

Figure 3

 3 Figure 3.2: Distribution of chip resources for a CPU vs a GPU. Figure from [Nvidia 2007].

Figure 3 . 3 :3

 33 Figure 3.3: GPU memory hierarchy. Each thread has private local memory, each block of thread has shared memory available for all the threads in the block, All threads have acces to the same global memory. From [Nvidia 2007].

Figure 4 :

 4 Figure 4: Extension of the grid to grid method including the w correction.

Figure 5 :

 5 Figure 5: Variation of the w coordinate value in the uv coverage. Different shades of grey are used for a different set of w values. w varies slowly in time during the observation. A simulation of a VLA-D 8 hours observation pointed at RA=19 : 25 : 59, Dec= 21.06.26.



Figure 6 :

 6 Figure 6: uv-coverage for a single baseline and the corresponding ellipse (dashline).

Figure 8 :

 8 Figure 8: Mean absolute error of the gradient computation regarding the oversampling factor K.

Figure 9 :

 9 Figure 9: Maximum absolute error of the gradient computation regarding the oversampling factor K.

Figure 10 :Figure 11 :

 1011 Figure 10: CPU -G2G computation time regarding the oversampling factor K.

Figure 12 :

 12 Figure12: GPU kernel -GPS G2G regarding the oversampling factor K.

Figure 16 :

 16 Figure 16: Roofline model.

  Figure 3.4: CPU -G2G scaling.

Figure 3 Figure 3

 33 Figure 3.5: CPU -G2G computation time regarding the oversampling factor K.

Figure 3

 3 Figure 3.7: CPU -Giga grid point read-write comparaison between G2G and STD methods.

  a new phase center gives the following dirty image y equationy(l -l φ , m -m φ ) = M i=1

Figure 4 . 1 :

 41 Figure 4.1: If the visibilities are phase-shifted to the facet center without baseline rotation, the image plane is phase shifted to the location of the new phase center.The new image plane remains parallel to the original image plane. Image adapted from Rick Perley[Taylor 1999]

Figure 4 . 2 :

 42 Figure 4.2: Polyhedron imaging aims to approximate the celestial sphere by a polyhedron. Image adapted from Rick Perley [Taylor 1999].

  Fig 4.5 shows an example of a compute queue connected to a process pool of 2 workers. If the job returns a result, the worker pushes it to the result queue.

Algorithm 6 :Figure 4

 64 Figure 4.3: DDFacet sequantial imaging pipeline.

Figure 4 Figure 4

 44 Figure 4.6: Existing DDFacet multi-node parallelization. Each node generates an independent image using different dataset of a specific observation.

Figure 4 Figure 4

 44 Figure 4.9: Scalability α 1 (n) of the gridding/degridding using one MPI process per node.

Figure 4 Figure 4

 44 Figure 4.11: Scalability α 4 (n) of the gridding/degridding using 4 MPI processes per node.

Figure 4

 4 Figure 4.13: Total execution time SpeedUp regarding the number of nodes used.

Figure 4

 4 Figure 4.14: Profiling of the master MPI process for a complete execution regarding the number of MPI processes per node (1, 2, and 4), and the number of nodes used for the parallel execution.

Fig 4 .

 4 Fig 4.15 shows the gridding speedup and the total execution time for each node between 4 MPI processes and 1 MPI process, such as the gridding speedup is defined as

  Figure 4.15: Speedup of gridding/degridding and total execution time between 4 MPI processes and 1 MPI process regarding the number of nodes.

  Figure 4.16: Parallelization efficiency regarding the number of nodes.

Figure 4 .

 4 Figure 4.17: (Left) Dirty Image. (Right) Restored image after 20.000 iterations of MSMFS Clean.

ii Contents Fast Grid to Grid interpolation 35 3

  .1 Multi-core and many-core architectures . . . . . . . . . . . . . . . . 36 3.1.1 Multi-core CPU architecture . . . . . . . . . . . . . . . . . . 36 3.1.2 Many-core GPU architecture . . . . . . . . . . . . . . . . . . 37 3.2 Fast Grid to Grid (G2G) interpolation . . . . . . . . . . . . . . . . . 38 3.3 CPU experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.1 CPU implementation . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.2 CPU Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 The Direction-Dependent Facets Imager (DDFacet) . . . . . . . . . . 60 4.1.1 Non-coplanar and coplanar facet imaging . . . . . . . . . . . 60 4.1.2 DDFacet specificities and features . . . . . . . . . . . . . . . 62 4.2 McMn parallelization . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.2.1 Sequential algorithm . . . . . . . . . . . . . . . . . . . . . . . 63 4.2.2 Multi-core parallelization . . . . . . . . . . . . . . . . . . . . 64 4.2.3 Multi-node parallelization with data distribution . . . . . . . 64 4.3 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

	Multi-core multi-node parallelization	59
	4.1 General conclusion	81

4.3.1 Setup and test case . . . . . . . . . . . . . . . . . . . . . . . . 71 4.3.2 Gridding and degridding speed-up α . . . . . . . . . . . . . . 72 4.3.3 Global speed-up b eta . . . . . . . . . . . . . . . . . . . . . . . 74 4.3.4 Analyses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.3.5 Memory footprint . . . . . . . . . . . . . . . . . . . . . . . . 77 4.4 5.1 Summary of the contributions . . . . . . . . . . . . . . . . . . . . . . 81 5.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

List of Figures

2.1 (left) Ruber's dish radio telescope, https://en.wikipedia.org/wiki/ Radio_telescope. (right) Arecibo dish radio telescope, https://en. wikipedia.org/wiki/Arecibo_Telescope. . . . . . . . . . . . . . . 2.2 Stations of the LOFAR radiotelescope. . . . . . . . . . . . . . . . . . 2.3 Image quality of SKA-mid (left) and the best current facility operating at the same frequency range of the JVLA (right). From [ska 2021].

3 Challenges of the new generation of radio interferometer . 27 2

Contents 2.1 Context of study . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.1 Radio Interferometer . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.2 Radio Interferometric Measurement Equation . . . . . . . . . 9 2.1.3 Synthesis imaging . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.4 Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Radio Interferometric Imaging . . . . . . . . . . . . . . . . . 16 2.2.1 Inverse Problem . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2.2 CLEAN -the historical method . . . . . . . . . . . . . . . . . 20 2.2.3 Other methods . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2..3.1 The case of SKA . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.3.2 HPC challenges . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3.3 State-of-the-art . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.

4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

  w, a i , ch = getNhit(i, nhit) ; gridU, gridV = getGrid(u t , v t ); convU, convV = getConv(u t , v t ); if gridU! = prev gridU or prev gridV! = gridV then atomicAdd(Ogrid[prev gridU][prev gridV], add i );

Table 1 :

 1 The experimental results are presented in table 1 for several oversampling factors. Memory footprint needed to compute the gradient in GB.

		K	8	16	24	32	40	64.
	STD	single double 7.7 3.8	3.8 7.7	3.8 7.7	3.8 7.7	3.8 7.7	3.8 7.7
	G2G 64 slices double 3.27 3.42 3.56 3.70 3.83 4.22 single 1.7 1.85 1.99 2.13 2.27 2.66
	G2G 1 slice	single 0.17 0.32 0.46 0.60 0.73 1.12 double 0.19 0.34 0.48 0.62 0.76 1.15

Table 2 :

 2 Theoretical c α and experimental c α compression ratio.

  c xy , c yx , c yy = (0,0); , v t , w, a i , ch = getNhit(i, nhit) ; /* Get oversampled coordinates u t , v t , the w value, the frequency channel, and a i for the i th element of the G2G sparse vector diag(A * ).Ogrid[u + u tap ][v + v tap ][ch][yx] += c yx × w; Ogrid[u + u tap ][v + v tap ][ch][yy] += c yy × w; end end end the DDEs related to the new generation of antenna thanks to the 3rd generation calibration. We propose a new Multi-core Multi-node (McMn) implementation to run the DDFacet imager on single or multi-node architectures with shared and distributed memory systems.

	2 for i (in parallel) M' do
			*/
	4	u, v = scale(u t , v t ); /* Scale u t and v t regarding the pixel
		resolution.	*/
	5 6 7	c xx , c xy , c yx , c yy = (0,0); for u tap ← [-half _size] to [+half _size] do for v tap ← [-half _size] to [+half _size]; /* Loop over the
		convolution kernel size. do	*/
	8	convU, convV = getConv(u tap , v tap , u , v, K); /* Get the
		convolution kernel indexes for the current iteration.
		*/
	9	wc = convFctConj[w][convU][convV] ; /* Get the complex
		convolution kernel value.	*/
		c xx += Igrid[u + u tap ][v + v tap ][ch][xx] × wc;
		c xy += Igrid[u + u tap ][v + v tap ][ch][xy] × wc;
		c yx += Igrid[u + u tap ][v + v tap ][ch][yx] × wc;
		c yy += Igrid[u + u tap ][v + v tap ][ch][yy] × wc; end end
		c xx = c xx × a i ; /* Element-wise multiplication between a i and
		its corresponding oversampled pixel.	*/
		c xy = c xy × a i ;
		c yx = c yx × a i ;
		c yy = c yy × a i ; for u tap ← [-half _size] to [+half _size] do for v tap ← [-half _size] to [+half _size] do
		convU, convV = getConv(u tap , v tap , u , v, K); /* Get the
		convolution kernel indexes for the current iteration.
		*/
		w = convFct[w][convU][convV] ; /* Get the complex
		convolution kernel value.	*/
		Ogrid[u + u tap ][v + v tap ][ch][xx] += c xx × w; /* Write in the
		output pixel grid.	*/
		Ogrid[u + u tap ][v + v tap ][ch][xy] += c xy × w;

3 u t

Table 4 .

 4 1: Maximum Memory consumption during gridding/degridding.

x ← x * clean P SF + δy;In order to save memory and accelerate the algorithm, the Clark-CLEAN[Clark 1980] variant was developed and is simplified in Alg. 2. The main idea is to use an approximation of the PSF by keeping only the central part of the PSF. This approximation remains relatively accurate because the sidelobes fall very quickly, and there is not much information left beyond them, especially if the uv-coverage is good. However, this ultimately leads to artifacts in the model. One uses masks to isolate the area of expected emissions to prevent these artifacts. Moreover, the residual is re-computed at each iteration to remove the potential remaining artifacts by

x ← x * clean P SF + δy;

https://casaguides.nrao.edu/index.php?title=VLA CASA Imaging-CASA6.2.0

https://docs.python.org/3/library/concurrent.futures.html

https://mpi4py.readthedocs.io/en/stable/mpi4py.futures.html

nous avons travaillé sur tous les aspects algorithmiques et implémentation pour diminuer le temps d'exécution des méthodes de reconstruction du ciel. Cela est parti des coûteux opérateurs de gridding et degridding jusqu'à des optimisations d'implémentation sur plateforme accélératrice et serveurs HPC.

Acknowledgments

Acknowledgment

This work was supported by grants from Région Ile-de-France and ANR DARK-ERA (ANR-20-CE46-0001-01).

of vec1 and vec2, respectively c 1 and p 1 complex values. After the initialization of the 256-bits vectors l.1, l2, and l.3, we have 4. vec3 = (c 1r p 1r , c 1i p 1i ) 5. vec2 = (p 1i , p 1r ) 6. vec2 = (p 1i , -p 1r ) 7. vec4 = (c 1r p 1i , -c 1i p 1r ) 8. output = (c 1r p 1r -c 1i p 1i , c 1r p 1i + c 1i p 1r ) // Here l.8 and l9 are fuse for simplification.

The output matches the result of the Eq. 3.2.

Algorithm 4: Quad polarization single float complex multiplication AVX instructions.

Data: (c 1 , c 2 , c 3 , c 4 , p 1 ) //with c the quad polarized complex value, p the complex weight.

Chapter 4

Multi-core multi-node parallelization Sky image reconstruction from the data is done using an imaging pipeline called an imager. An imager uses image reconstruction algorithms whose algorithmic blocks are highly differentiable from one another. Applying direction-dependent corrections during the gridding and degridding operators is computationally very expensive. [Tasse 2018] have shown that these operators can reach 80% of the total execution time. We have presented the state-of-the-art and a method to speed up the computation time for the gridding and degridding operators in chapters 2 and 3. However, facing tomorrow's computational challenges, distributing the computations over several computational resources becomes mandatory. Several imagers already allow parallelizing the computations on several CPUs, like CASA and WS-Clean. In our case, we are interested in DDFacet, an imager that allows to correct the deconvolution algorithm to the residual image δy, as the MSMFS algorithm is hardly parallelized on multi-node. The resulting hypercube x (k) is then broadcast to all nodes. Each node receives only the slices whose frequency band is associated

Chapter 4: Multi-core multi-node parallelization

The pseudo-code of the MPI implementation is described by Alg. 7.

Experiments

This section presents the experimental results of McMn parallelization of DDFacet on a distributed memory HPC server. We will first look at the hardware and the test case used. Then, we will present the results.

Setup and test case

The test case is an 8-hours observation performed with the LOFAR radio telescope covering the 75 -120 MHz frequency band. In order to increase the amount of data used, we artificially inflated the number of MS files to be processed. With an initial number of 24 MS files, for a total of 500 GB, we increased to 48 MS files, for a total of 1.1 TB of data to process. The image to be reconstructed was 10.000 × 10.000 pixels with a single slice for all frequencies. The stopping parameters of the imaging framework were 20.000 cumulative iterations of minor cycles made with the MSMF CLEAN algorithm. After experimentation, these 20.000 iterations were reached after three major cycles. The image was divided into 11 × 11 = 121 facets for this test case. The image corresponding to this test case is illustrated by Fig. 4.17, where the image on the left is the dirty image, and the image on the right is the reconstructed image.

The experiments were done on a BullSequana XH2000 system based on AMD EPYC 7742 64-core processors. Each node is a dual-socket node which, as illustrated by Fig. 4.8, is two processors connected by a PCI-express link. Both CPUs can read and write to the RAM of the other chip. Moreover, they share the same local memory inside the node. In practice, it can be seen as a big processor with double the cores. Thus, in our case, we have nodes with 128 physical cores. Hyperthreading on Intel CPUs or Simultaneous Multi-Threading (SMT) on AMD CPUs are named for the process of separating the physical core into virtual cores, allowing a core to do multi-threading more efficiently. This optimization allows gains of up to 50% on the efficiency of core usage. When Hyperthreading is active, each core can manage two threads. In our case, AMD processors allow to use of SMT, and the Operating System (OS) sees 256 logical cores.

Our test case aims at reconstructing an image of 121 facets from nodes with 256 logical-cores processor. As the multi-core parallelization is facet based, using 1 MPI process on the node, we will not be able to use all the available resources of the node. Indeed, with a single MPI process, we only use 121 cores out of the 256 available. However, using only 121 cores out of the 256 available logical cores means that the potential gains from using SMT are not achieved. With 2 MPI processes per node, we still arrive at an underutilization with 242 logical cores used out of the 256. We can fully exploit all available resources with 4 MPI processes per node. Thus, we obtained results for the experiments using 1 MPI process per node, 2 MPI processes per node, and 4 MPI processes per node. For the cases with 2 and 4 MPI processes, we made sure to bind the cores used by each process on a specific processor of the dual-socket node in order to minimize the memory transfers to the RAM of the other socket.

In order to balance the computational load among the nodes, the number of MS files to be processed by each MPI process in a single run must be the same. As the test case is 48 MS files, using 1 MPI process per node, parallelization could be

Publications

Softwares

• DDFacet, A facet-based radio imaging package. Github : https: //github.com/saopicc/DDFacet, branch : mpi_mergeMaster.

• Grid to Grid, MOnnier Radio Interferometric Software (MORIS), Github : On-demand.

Titre:

ExaSKA : Parallélisation sur serveur de de calcul intensif pour le radiotelescope SKA

Mots clés:

Calcul Parallèle, Radio Astronomie, Problème Inverse, Interpolation, Exascale, GPU

Résumé:

La

Title:

ExaSKA : Parallelization on a High Performance Computing server for the exascale radiotelescope SKA

Keywords:

Parallel Computing, Radio Astronomy, Inverse Problem, Interpolation, Exascale, GPU

Abstract:

Radio astronomy allows us to make observations of the sky at wavelengths invisible to the naked eye and thus obtain a lot of information about the universe that surrounds us. The antenna arrays that make up these observation instruments generate a continuous data flow at rates that make their processing and storage almost impossible. In addition to the very large amount of data to process, the reconstruction of sky images requires iterative algorithms that are costly in terms of time and computational power, making the largest supercomputers insufficient to do it in real-time. This thesis explores different ways to reduce the computation time for image reconstruction by focusing on the algorithmic aspect and the implementation aspect. We have proposed an algorithmic method to reduce the computational cost of expensive operators. We have also proposed an implementation on an accelerator platform and a large-scale implementation on several computational nodes of a supercomputer.

Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France