Oleg Mitrofanov
email: o.mitrofanov@ucl.ac.uk

Filip Dominec

Petr Kužel

John L Reno

Igal Brener

U-Chan Chung

Cathy Elissalde

Mario Maglione

Patrick Mounaix

Near-field probing of Mie resonances in single TiO 2 microspheres at terahertz frequencies

Keywords:

We show experimentally that poly-crystalline TiO 2 spheres, 20-30 μm in diameter, exhibit a magnetic dipole Mie resonance in the terahertz (THz) frequency band (1.0-1.6 THz) with a narrow line-width (<40 GHz). We detect and investigate the magnetic dipole and electric dipole resonances in single high-permittivity TiO 2 microspheres, using a near-field probe with a sub-wavelength (~λ/50) size aperture and THz time-domain spectroscopy technique. The Mie resonance signatures are observed in the electric field amplitude and phase spectra, as well as in the electric field distribution near the microspheres. The narrow line-width and the subwavelength size (λ/10) make the TiO 2 microspheres excellent candidates for realizing low-loss THz metamaterials.

Introduction

The concept of metamaterials has enabled design of electromagnetic functionalities, which do not occur in natural materials. This has led to the realization of intriguing electromagnetic phenomena, such as cloaking [START_REF] Schurig | Metamaterial electromagnetic cloak at microwave frequencies[END_REF] and negative refraction [START_REF] Smith | Metamaterials and negative refractive index[END_REF]. One of the key elements for metamaterials is a sub-wavelength size resonator with an effective magnetic response. Metallic split-ring resonators (SRR) are often used for that purpose [START_REF] Tao | Recent progress in electromagnetic metamaterial devices for terahertz applications[END_REF]. However, the SSRs show a less than ideal behavior due to their anisotropy, polarization dependent response and spectral broadening caused by metallic Ohmic losses. For these reasons, arrays of dielectric resonators have attracted much interest recently. A magnetic response can be produced in a high-permittivity, sub-wavelength size dielectric particle as a result of the lowest order resonance, known as the 1st Mie mode or the magnetic dipole (MD) mode [START_REF] Merlin | Metamaterials and the Landau-Lifshitz permeability argument: Large permittivity begets highfrequency magnetism[END_REF][START_REF] Kuznetsov | Magnetic light[END_REF][START_REF] Geffrin | Magnetic and electric coherence in forward-and back-scattered electromagnetic waves by a single dielectric subwavelength sphere[END_REF]. An array of such resonators exhibits a macroscopic magnetic response with negative values of permeability μ in the wavelength range determined by the resonator size [START_REF] Lewin | The electrical constants of a material loaded with spherical particles[END_REF][START_REF] O'brien | Photonic band-gap effects and magnetic activity in dielectric composites[END_REF][START_REF] Popa | Compact dielectric particles as a building block for low-loss magnetic metamaterials[END_REF][START_REF] Zhao | Experimental demonstration of isotropic negative permeability in a three-dimensional dielectric composite[END_REF][START_REF] Zhao | Mie resonance-based dielectric metamaterials[END_REF][START_REF] Vynck | All-dielectric rod-type metamaterials at optical frequencies[END_REF][START_REF] Ginn | Realizing optical magnetism from dielectric metamaterials[END_REF][START_REF] Němec | Tunable terahertz metamaterials with negative permeability[END_REF][START_REF] Shibuya | Terahertz metamaterials composed of TiO 2 cube arrays[END_REF]. Similarly, a negative dielectric response ε<0 can be realized using the second Mie mode, known as the electric dipole (ED) mode [START_REF] Kuznetsov | Magnetic light[END_REF][START_REF] Geffrin | Magnetic and electric coherence in forward-and back-scattered electromagnetic waves by a single dielectric subwavelength sphere[END_REF][START_REF] Lewin | The electrical constants of a material loaded with spherical particles[END_REF][START_REF] O'brien | Photonic band-gap effects and magnetic activity in dielectric composites[END_REF][START_REF] Popa | Compact dielectric particles as a building block for low-loss magnetic metamaterials[END_REF][START_REF] Zhao | Experimental demonstration of isotropic negative permeability in a three-dimensional dielectric composite[END_REF][START_REF] Zhao | Mie resonance-based dielectric metamaterials[END_REF][START_REF] Vynck | All-dielectric rod-type metamaterials at optical frequencies[END_REF][START_REF] Ginn | Realizing optical magnetism from dielectric metamaterials[END_REF][START_REF] Shibuya | Terahertz metamaterials composed of TiO 2 cube arrays[END_REF].

High-permittivity materials (ε>20) are rare in the optical part of the electromagnetic (EM) spectrum. At terahertz (THz) frequencies however, some materials, for example rutile compounds, exhibit ε>100 and moderate losses [START_REF] Shibuya | Terahertz metamaterials composed of TiO 2 cube arrays[END_REF]. The MD resonance in the region of 1 THz is expected for a microsphere, 30 μm in diameter, made of TiO 2 (ε≈100). An ensemble of such TiO 2 microspheres was recently shown to exhibit a resonant spectral signature due to the MD Mie mode at 0.7-0.9 THz [START_REF] Němec | Resonant magnetic response of TiO 2 microspheres at terahertz frequencies[END_REF]. The inhomogeneous line-width broadening of 0.2 THz however obscured the intrinsic resonator properties. The resonance frequency of each microsphere depends sensitively on the resonator size and shape, as well as the crystallinity and porosity of TiO 2 , all of which varied within the ensemble [START_REF] Němec | Resonant magnetic response of TiO 2 microspheres at terahertz frequencies[END_REF].

In order to determine the intrinsic line-widths of the MD and ED modes, the TiO 2 microspheres should be probed individually. We estimate that a single TiO 2 sphere of radius r = λ/(2n) = λ/20 (assuming ε = n 2 = 100) placed in the focus of a typical THz beam produces small attenuation (~1%) of the beam intensity at the frequency of MD resonance. This small value of attenuation is due to the reduction of the scattering cross-section at the MD resonance for high-permittivity resonators. In the vicinity of a high-permittivity resonator, on the other hand, the electromagnetic energy of the incident wave is concentrated. Near-field microscopy methods therefore may provide an alternative method for probing the resonator modes. Concentration of THz surface (Zenneck) waves near high-permittivity subwavelength size particles, including a TiO 2 microsphere, was recently observed with THz near-field microscopy [START_REF] Navarro-Cía | Terahertz imaging of sub-wavelength particles with Zenneck surface[END_REF]. Spectroscopic signatures of the Mie modes however were not detected [START_REF] Navarro-Cía | Terahertz imaging of sub-wavelength particles with Zenneck surface[END_REF].

Here, we investigate the MD and ED modes in single TiO 2 resonators using THz nearfield microscopy in the transmission configuration [START_REF] Mitrofanov | Collection mode near-field imaging with 0.5 THz pulses[END_REF][START_REF] Knab | Terahertz near-field spectroscopy of filled subwavelength sized apertures in thin metal films[END_REF][START_REF] Mitrofanov | Probing terahertz surface plasmon waves in graphene structures[END_REF]. We exploit the effect of the resonator on transmission of a plane wave through a sub-wavelength size aperture [START_REF] Aydin | Split-ring-resonator-coupled enhanced transmission through a single subwavelength aperture[END_REF][START_REF] Sadaune | Enhanced transmission via a sub-wavelength hole aperture coupled to a high-permittivity resonator[END_REF]. We observe clear signatures for the MD and ED modes by applying THz time-domain spectroscopy analysis. The MD mode displays a narrow resonance, several times narrower than the line-width observed in the far-field ensemble study [START_REF] Němec | Resonant magnetic response of TiO 2 microspheres at terahertz frequencies[END_REF].

Results and discussion

Experimental

For our studies, poly-crystalline TiO 2 microspheres are prepared by assembling and annealing clusters of TiO 2 nano-crystals [START_REF] Shibuya | Terahertz metamaterials composed of TiO 2 cube arrays[END_REF]. The nano-crystals were mixed with ethanol to obtain a liquid suspension, which was then dried by spraying it through a flame. In this process the nano-crystals assemble into nearly spherical clusters, which solidify after annealing at 1200° C for 2 hours. The annealing time and temperature affect the packing density of nano-crystals in the clusters. Scanning electron micrographs show that the nano-crystals within the clusters are densely packed [Fig. 1(a)]. There are small voids between randomly oriented nanocrystals. Due to the sub-wavelength size of the nano-crystals and voids, and their random orientation, the dielectric function can be assumed to be isotropic. After comparing the size of similar microspheres to their resonance wavelength, we estimate the relative dielectric constant to be in the region of 90-100 [START_REF] Němec | Resonant magnetic response of TiO 2 microspheres at terahertz frequencies[END_REF].

Microsphere samples ranging in diameter between 20 and 30 μm were picked from the batch randomly and mounted on a 12.5 μm thick low-dielectric constant polymer film (polyvinyl chloride) with a thin layer of adhesive (cyanoacrylate). The thin-film substrate is chosen to reduce its effect on the resonator properties. The microsphere position is controlled by an automated XYZ translation stage.

The sample is illuminated by a ~0.5 mm diameter beam of single-cycle THz pulses (0.5-2.5 THz) generated by the process of optical rectification in a ZnTe crystal. A sub-wavelength size aperture in a 300 nm Au film on GaAs is positioned below the sphere as illustrated schematically in Fig. 1(b) [START_REF] Macfaden | 3 μm aperture probes for near-field terahertz transmission microscopy[END_REF]. The aperture size of 5 μm or 10 μm is chosen to be several times smaller than the microsphere diameter and substantially smaller than the wavelength (120-600 μm). A photo-conductive THz antenna is integrated into the system in the near-field zone of the aperture [Fig. 1(b)], ~1 μm below the aperture plane [START_REF] Macfaden | 3 μm aperture probes for near-field terahertz transmission microscopy[END_REF]. The antenna is oriented along the polarization direction of the incident THz beam (x-axis).

THz near-field spectroscopy and imaging of TiO 2 spheres

The effect of the sphere on the THz beam transmission through the aperture can be understood as the localization of the resonant field near the sphere due to (1) the reduced effective wavelength and (2) energy trapping inside the resonator. The localized field couples through the sub-wavelength aperture more efficiently than the incident THz beam, resulting in enhanced transmission at the frequency of the Mie resonance. Recent numerical studies showed that the power transmission through a λ/20 aperture can be enhanced by up to 50 times when the dielectric resonator is positioned in front of the aperture [START_REF] Kang | Numerical analysis of enhanced transmission through a single subwavelength aperture based on Mie resonance single particle[END_REF]. This effect was utilized to probe a dielectric cube resonator at microwave frequencies [START_REF] Sadaune | Enhanced transmission via a sub-wavelength hole aperture coupled to a high-permittivity resonator[END_REF].

The transmission enhancement effect is observed in the electric field map E x (x,y) measured by scanning a spherical TiO 2 particle ~27 μm in diameter (Sample T1) with respect to the aperture in the xy-plane at a distance z≈1 μm [Fig. 1(c)]. The image displays an enhanced electric field (negative value) in the region of the sphere at a moment immediately after the excitation by the THz pulse (t = 1.33 ps). The region of enhanced field has an oblong shape with the long axis aligned to the incident field polarization direction (x-axis). The oblong shape reflects the electric field distribution near the resonator rather than the resonator's slightly elongated shape. It will be shown later using numerical simulations that an ideal sphere produces a similar oblong field pattern. We note that the aperture size is several times smaller than the size of the enhanced field region and therefore the aperture shape (5 μm x 5 μm square) is not reflected in the observed images. The THz image shows weaker lobes above and below the sphere, corresponding to the x-component of the field switching its sign outside the sphere.

Analysis of the THz pulse time-domain waveforms reveals that the enhanced transmission through the aperture occurs only at selected frequencies. The waveform measured by the 5 μm aperture probe when the microsphere is directly above the aperture [Fig. 1(d), red trace] displays coherent regular oscillations lasting over several cycles. The oscillations are absent in the waveform of the incident reference pulse [Fig. 1(d), black trace], measured ~50 μm away from the sphere, at a location marked by the black cross, where the effect of the sphere is negligible. The oscillations indicate that a resonance is formed inside the sphere.

The ratio of the Fourier transform spectra E ms /E 0 shows that the transmission amplitude decreases close to 0 at 1.1 THz and then rises to the level of ~12 (~150 in intensity) at a slightly higher frequency (1.175 THz). The observed transmission amplitude signature, similar to recent microwave scattering parameter measurements on a high-permittivity resonator [START_REF] Sadaune | Enhanced transmission via a sub-wavelength hole aperture coupled to a high-permittivity resonator[END_REF], follows a Fano line shape [START_REF] Luk'yanchuk | The Fano resonance in plasmonic nanostructures and metamaterials[END_REF]:

( ) ( ) ( ) 2 0 2 2 0 0 2 , / 2 ms q E E ω ω ω ω Γ   + -     = Γ + - (1)
where the Mie resonance frequency ω 0 /2π = 1.175 THz and the resonance line-width parameter Γ/(2π) = 0.02 ± 0.01 THz [Fig. 1(e)]. q is the Fano parameter. It quantifies the aperture transmission enhancement factor at the resonance frequency: E ms /E 0 = q. The resonance line-width for the observed asymmetric line shape can be estimated by fitting the Fano line shape to the data. The Fano shape (Eq. ( 1)) changes to the Lorentzian shape for q = 0. The MD resonance full width at half maximum, which would be measured in a far-field transmission measurement, therefore is expected to be equal to Γ/(2π) = 20 ± 10 GHz. We note however that the fitted line-width is similar to the measurement resolution line-width of 40 GHz, which we estimate from the duration of the recorded waveform (25 ps). Therefore we can conclude only that the observed line-width for Sample T1 is less than 40 GHz.

The spectrum also shows a smaller peak at ~1.7 THz. The set of two resonances matches the lowest order Mie modes, the MD and ED modes [Figs. 2(a) and 2(b)], simulated for a spherical particle with ε≈100. The simulations presented in Fig. 2 are performed with the finite-difference time-domain method (FDTD) using a freely available software package [START_REF] Oskooi | MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method[END_REF]. Similar spectra with the first resonance frequency ranging from 1.1 to 1.6 THz and the second frequency ranging from 1.6 to 2.0 THz are measured for five nearly spherical TiO 2 particles, ranging in size from 20 μm to 30 μm. The ratio of the sphere radius to the MD resonance wavelength r/λ≈0.05 is consistent with the assumption that the average dielectric function ε≈100.

The observed resonances confirm that the enhanced coupling effect is due to the excitation of normal modes in the microsphere (TiO 2 has no material-related modes in this spectral range). The resonance frequencies decrease as the resonator size increases, whereas their ratio ω MD /ω ED ≈0.7 varies slightly from sample to sample. We attribute this variation to the non-ideal shape of the samples. Figures 2(c) and 2(d) illustrate the dependence of the resonance frequencies for an ellipsoid with semi-axes a, b and c.

The EM simulations in Figs. 2(a) and 2(b) also show that the modes are confined within the sphere, with only a small fraction of the mode energy leaking outside. To detect this confined field, the near-field probe must be positioned within a short distance from the sphere. The mode confinement is illustrated in the map of the transmission spectra (Sample T2) measured for z ranging from 0 to 40 μm [Fig. 3(a)]. As the microsphere is moved away from the aperture, the detected signatures of the MD resonance at 1.4 THz and the ED resonance at 2.0 THz decrease in amplitude. By fitting the Fano function to the experimental data [Fig. 3(b)], we find that the enhancement factor for the MD mode, q, decreases rapidly within short distance from the sphere, 0<z<15 μm [Fig. 3(a)]. Fitting a Fano line shape to the data in Fig. 3 allows us to quantify the effect of the probe on the resonator characteristics. The line-width parameter for the MD mode remains constant within the experimental error: Γ/(2π) = 0.05 ± 0.01 THz. The resonance frequency however shifts slightly from f 0 = 1.43 THz for the case of z = 0 μm to 1.37 THz for z>15 μm. Accuracy of the extracted parameters for the MD resonance is however limited due to the measurement resolution and deviation of the experimental line shape from the ideal Fano line shape.

To verify the observed trends, we simulate the electromagnetic field distribution for an ideal sphere (ε = 100, r = 12.5 μm) placed in close proximity to a metallic plane with a 5 μm aperture. The simulations are performed with the FDTD method using a commercially available software package [27]. The electric field distribution in the xz-plane intersecting the sphere at its center is shown in Fig. 3(c) for two separations between the plane and the sphere: z = 1 μm (left) and z = 4 μm (right). The displayed field amplitude |E| is normalized to its maximum value inside the sphere. For the separation of 1 μm, the metallic plane noticeably perturbs the field distribution outside the sphere. Inside the sphere however, the effect of the metallic plane is weak: only small asymmetry appears in the azimuthally symmetric field distribution of the MD mode. For the separation of 4 μm, the field inside the resonator exhibits practically no asymmetry due to the metallic plane. The near-field probe therefore preserves the MD resonance in TiO 2 microspheres. The simulated transmission spectra for the aperture with and without the sphere allow us to compare the numerical modeling with the experimentally recorded spectra. The normalized amplitude spectra E ms /E 0 display remarkable similarities with the experiment. In Fig. 3(b), we illustrate the agreement between the experiment and simulations for a sphere of radius r = 11 μm and ε = 94.2 + 1.16i. The line shapes and the resonance frequencies for the MD and ED modes are in close agreement with the measured spectrum. The transmission enhancement factor at the MD resonance is slightly higher for the simulated spectrum. The enhancement factor however depends sensitively on the distance between the sphere and the aperture, as shown in the inset of Fig. 3(a).

The effect of the metallic screen on the MD resonance characteristics (frequency and linewidth) is small for TiO 2 resonators. The resonance frequency shift due to the metallic probe is only 0.3% in the simulations. It is attributed to the strong mode localization by the highpermittivity material. We note that the shift increases for smaller permittivity resonators. It was recently reported that the MD resonance frequency shifts by ~0.6% for a ε = 20 sphere resting on a metallic plane in comparison with the sphere in free space [START_REF] Markovich | Magnetic dipole radiation tailored by substrates: numerical investigation[END_REF]. The dielectric substrate can affect the resonance properties as well [START_REF] Markovich | Magnetic dipole radiation tailored by substrates: numerical investigation[END_REF], however the low-permittivity and the sub-wavelength thickness of the substrate used in our experiments produce no noticeable effect. This is verified by comparing spectra measured for one sphere placed directly on the aperture screen with and without the polymer film substrate [Fig. 1(b)].

Effects of resonator non-sphericity

The MD mode in Sample T2 has a broader line-width than that in Sample T1 by a factor of 2 and it displays a 'shoulder' on the high-frequency side of the resonance. We attribute this line-width broadening to the resonator shape. All samples in this study deviate slightly from a perfect sphere [Figs. 1(e) and 3(b)]. While the MD resonance in an ideal sphere is 3-fold degenerate, the degeneracy is lifted in an ellipsoid. Three orthogonal MD modes can be excited if the incident wave polarization is not aligned with any of the ellipsoid principal axes. The spectral response of the resonator in this case is a superposition of projections of these modes on the detector antenna axis. The observed MD line-width therefore can be broadened by the resonator non-sphericity.

To verify the presence of MD modes with different orientations we use spatial mapping in the xy-plane at various moments after the THz pulse excitation. Panel (a) in Fig. 4 shows a THz image of a TiO 2 resonator (Sample T3) at the moment immediately after the THz pulse arrival. The field profile has an oblong shape, as the profile in Fig. 1(c). A similar shape is predicted by the FDTD simulations for an ideal sphere, as shown in the electric field map in Fig. 4(b), simulated for a moment immediately after the short pulse excitation.

Approximately 5 ps after the incident pulse, the experimental E x map displays a very different pattern, anti-symmetric with respect to the y-axis [Fig. 4(c)]. A mode of this symmetry can be excited by the incident THz pulse only in the case of a non-spherical resonator. Such a pattern does not appear in the simulations of an ideal spherical resonator. The origin of the pattern in Fig. 4(c) is determined by analyzing the detected THz spectra at two locations marked by the red and blue arrows in the images: at the resonator center (x 0 ,y 0 ); and at the anti-symmetric mode maximum (x 0 , y 0 -10 μm). In both cases, the spectra show resonances in the frequency region of the MD mode (~1.4 THz), however the line-shape for the anti-symmetric mode is reversed [Fig. 4(d)]. The change of the line shape correlates with the phase spectra shown in [Fig. 4(e)]. The spectral signatures confirm the presence of two modes.

The resonance frequency of the anti-symmetric mode (f 2 = 1.36 THz) is slightly different from that of the symmetric mode (f 1 = 1.42 THz), indicating that the MD mode degeneracy is lifted. While the strong resonance at f 1 corresponds to the MD mode with a large projection of its magnetic moment on the y-axis, the weaker resonance at f 2 correspond to a MD mode with a different orientation. The dipole moment for f 2 must have a non-zero projection on the zaxis, as the observed anti-symmetric field profile in Fig. 4(c) matches the corresponding E x distribution [inset of Fig. 4(e)].

The small difference in their resonance frequencies allows us to detect the presence of the two MD resonances with different orientations of their magnetic dipole moments in the temporally resolved field distribution. After the THz pulse excitation, the relative phase between the coherent oscillations of the two dipole moments increases over time. For f 1 -f 2 = 60 GHz, a π/2 phase shift is accumulated ~4-5ps after the excitation. Around that moment in time, the oscillating field of one of the modes is at its crest, while the field of the other mode is near its node. The required time for gaining the π/2 phase shift (for f 1 = 1.42 THz f 2 = 1.36 THz) matches the time delay t = 5 ps after the excitation for the field distribution in Fig. 4(c). This analysis supports the hypothesis that the field distribution captured in Fig. 4(c) corresponds to the field profile of a MD mode with a large projection of its magnetic moment on the z-axis.

Conclusion

The poly-crystalline TiO 2 microspheres investigated in this study exhibit the MD and ED Mie resonances in the range of 1.0-1.6 THz and 1.6-2.1 THz respectively. The MD resonance linewidth for a single resonator is less than 40 GHz. The narrow line-width and the subwavelength size (2r≈λ/10) make the TiO 2 microspheres excellent candidates for realizing the effective magnetic response for low-loss THz metamaterials. In comparison with SRR, the TiO 2 resonators can provide a narrower line-width at THz frequencies and eliminate the anisotropy and polarization dependence of the effective magnetic permeability μ [START_REF] Tao | Recent progress in electromagnetic metamaterial devices for terahertz applications[END_REF][START_REF] Singh | Sharp Fano resonances in THz metamaterials[END_REF].

The near-field scanning probe microscopy technique employed here enables us to probe individual high-permittivity TiO 2 resonators. The Mie resonance signatures, extracted using the THz time-domain spectroscopy method, appear as Fano resonances in the amplitude and phase spectra, whereas the electric field distribution near the microspheres helps identifying the resonance nature. The experimental observations allow us to correlate the variation of the MD resonance line-width with the resonator non-sphericity, which results in the splitting of the degenerate MD modes.

Fig. 1 .

 1 Fig. 1. (a) Scanning electron microscope image of a typical poly-crystalline TiO 2 sphere; (b) Schematic diagram of the experimental setup: a TiO 2 sphere attached to a thin polymer film (PF) is positioned at distance z above the sub-wavelength aperture (SWA) probe containing a photoconductive THz antenna [23]. (c) A map of the electric field (E x ) near a TiO 2 sphere immediately after the excitation (Sample T1). The dashed circle shows the sphere size and position. The black cross marks a location where the incident THz pulse waveform is measured for reference. (d) THz pulse waveforms detected in the near-field zone of the sphere (red trace) and at the cross location (black trace). (e) The amplitude transmission spectrum (red line) obtained from the waveforms in (d) and the Fano line shape (black line) fitted to the data; inset: an optical microscope image (60 μm × 60 μm) of Sample T1.

Fig. 2 .

 2 Fig. 2. (a) FDTD simulated MD (1.21 THz) and (b) ED (1.66 THz) modes in a TiO 2 sphere (ε = 100, radius r = 12.5 μm); the x-polarized incident EM wave propagates along the z-axis. The color represents the electric field (E x ) normalized to the maximum field within the sphere. (c) and (d) Resonance frequencies of the MD and ED modes in an elliptical resonator.

Fig. 3 .

 3 Fig. 3. (a) Dependence of the transmission amplitude spectrum E ms /E 0 on the distance z between the aperture and the microsphere (Sample T2); inset: the resonance frequency f 0 and the enhancement factor q as functions of z. (b) The amplitude spectrum for the closest distance z<1 μm measured by a 5 μm aperture probe (red line), the fitted Fano resonance line shape (black line) and the simulated spectrum of a TiO 2 sphere; inset: the optical microscope image (60 μm × 60 μm) of Sample T2. (c) The simulated field distributions for a sphere (r = 12.5 μm) placed in front of the near-field probe at the MD resonance; the metallic screen with a 5 μm aperture is 1 μm (left) and 4 μm (right) away from the sphere surface.

Fig. 4 .

 4 Fig. 4. Measured (a) and simulated (b) electric field maps E x (x,y) at z = ~1 μm from the TiO 2 sphere immediately after the excitation of the sphere by the THz pulse, t = 1.1 ps (image area: 100 μm × 100 μm). (c) E x (x,y) map of the same area at t = 5 ps. The spectral amplitude ratio (d) and phase (e) measured at the sphere center (x 0 , y 0 ) and at (x 0 , y 0 -10 μm); these location are indicated by the red and blue arrows in (a) and (c). The diagram in the inset of (e) shows the electric field distributions |E| and E x for the MD mode with the magnetic moment perpendicular to the xy-plane.

Acknowledgments This work was supported by the Royal Society [Grant No. UF080745], Czech Science Foundation [project No. 14-25639S] and by European Union funding under the 7th Framework Programme [project NOTEDEV]. Fabrication of THz near-field probes was performed, in part, at the Center for Integrated Nanotechnologies, an Office of Science User Facility operated for the U.S. Department of Energy (DOE) Office of Science. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under Contract No. DE-AC04-94AL85000.


Emmanuel Agullo
email: emmanuel.agullo@inria.fr

Olivier Beaumont
email: olivier.beaumont@inria.fr

Lionel Eyraud-Dubois
email: lionel.eyraud-dubois@labri.fr

Julien Herrmann
email: julien.herrmann@ens-lyon.fr

Suraj Kumar
email: suraj.kumar@inria.fr

Loris Marchal
email: loris.marchal@ens-lyon.fr

Samuel Thibault
email: samuel.thibault@labri.fr

Bridging the Gap between Performance and Bounds of Cholesky Factorization on Heterogeneous Platforms

Keywords:

We consider the problem of allocating and scheduling dense linear application on fully heterogeneous platforms made of CPUs and GPUs. More specifically, we focus on the Cholesky factorization since it exhibits the main features of such problems. Indeed, the relative performance of CPU and GPU highly depends on the sub-routine: GPUs are for instance much more efficient to process regular kernels such as matrix-matrix multiplications rather than more irregular kernels such as matrix factorization. In this context, one solution consists in relying on dynamic scheduling and resource allocation mechanisms such as the ones provided by PaRSEC or StarPU. In this paper we analyze the performance of dynamic schedulers based on both actual executions and simulations, and we investigate how adding static rules based on an offline analysis of the problem to their decision process can indeed improve their performance, up to reaching some improved theoretical performance bounds which we introduce.

I. INTRODUCTION

Linear algebra operations are the basis of many scientific operations. Our objective is to optimize the performance of one of them (namely the Cholesky decomposition) on a hybrid computing platform. The use of GPUs and other accelerators such as Xeon Phi are common ways to increase the computation power of computers at a limited cost. The large computing power available on such accelerators for regular computation makes them unavoidable for linear algebra operations. However, optimizing the performance of a complex computation on such a hybrid platform is very com-plex, and a manual optimization seems out of reach given the wide variety of hybrid configurations. Thus, several runtime systems have been proposed to dynamically schedule a computation on hybrid platforms, by mapping parts of the computation to each processing elements, either cores or accelerators. Among other successful projects, we may cite StarPU [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF] from INRIA Bordeaux (France), QUARK [START_REF] Yarkhan | QUARK Users' Guide: QUeueing And Runtime for Kernels[END_REF] and PaRSEC [START_REF] Bosilca | PaRSEC: A programming paradigm exploiting heterogeneity for enhancing scalability[END_REF] from ICL, Univ. of Tennessee Knoxville (USA), Supermatrix [START_REF] Chan | SuperMatrix: a multithreaded runtime scheduling system for algorithms-byblocks[END_REF] from University of Texas (USA), StarSs [START_REF] Planas | Hierarchical task-based programming with StarSs[END_REF] from Barcelona Supercomputing Center (Spain) or KAAPI [START_REF] Gautier | Kaapi: A thread scheduling runtime system for data flow computations on cluster of multi-processors[END_REF] from INRIA Grenoble (France). Usually, the structure of the computation has to be described as a task graph, where vertices represent tasks and edges represent dependencies between them. Most of these tools enable, up to a certain extent, to schedule an application described as a task graph onto a parallel platform, by mapping individual tasks onto computing resources and by performing data movements between memories when needed.

There is an abundant literature on the problem of scheduling task graphs on parallel processors. This problem is known to be NP-complete [START_REF] Garey | Computers and Intractability, a Guide to the Theory of NP-Completeness[END_REF]. Lower-bounds based either on the length of the critical path (the longest path from an entry vertex to an output vertex) or on the overall workload (assuming ideal parallelism) have been proposed, and simple list-scheduling algorithms are known to provide 2 -1/m-approximation on homogeneous platforms, at least when communication times are negligible [START_REF] Graham | Bounds for certain multiprocessing anomalies[END_REF]. Several scheduling heuristics have also been proposed, and among them the best-known certainly is heterogeneous early finish time (HEFT) [START_REF] Topcuoglu | Performanceeffective and low-complexity task scheduling for heterogeneous computing[END_REF], which inspired some dynamic scheduling strategies used in the above-mentioned runtimes. However, it remains a large gap between the theoretical lower-bounds and the actual performance of dynamic HEFT-like heuristics. Another way to assess the quality of a scheduling strategy is to compare the actual performance to the machine peak performance of the computing platform computed as the sum of the performance of its individual computational units. Rather than this machine peak performance which is known to be unreachable, one usually considers the GEMM peak obtained by running matrix multiplication kernels (GEMMs). For large matrices, the task-graph of a Cholesky factorization exhibits a sufficient amount of parallelism, and a sufficient number of GEMM calls for this bound to be reasonable. However, on small and medium size matrices, there is still a large gap between GEMM peak performance and the best-achievable Cholesky performance.

In this paper, we optimize the dynamic scheduling of the Cholesky decomposition of a dense, symmetric, and positive-definite double-precision matrix A, into the product LL T , where L is lower triangular and has positive diagonal elements, using one runtime system, StarPU, and provide better makespan bounds to prove the quality of our schedules. The contributions of the paper are:

• Better lower bounds on the makespan of a Cholesky factorization on a parallel hybrid platform; • Better dynamic schedules, based not only on HEFT but also on an hybridization of static and dynamic task assignments; • A very efficient schedule for a simple hybrid platform model, achieved by constraint programming.

• Numerous experiments to assess the performance of our schedules using the StarPU runtime. Note that what is done here using StarPU could have been done with other runtimes, provided that we are able to control their mapping and scheduling policies. Similarly, we could have chosen another dense linear algebra factorization such as the QR or LU decompositions.

II. CONTEXT A. Cholesky factorization

The Cholesky factorization (or Cholesky decomposition) is mainly used to solve a system of linear equations Ax = b, where A is a N ×N symmetric positive-definite matrix, b is a vector, and x is the unknown solution vector to be computed. Such systems often arise in physics applications, especially when looking for numerical solutions of partial differential equations, where A is positive-definite due to the nature of the modeled physical phenomenon. One way to solve such a linear system is first to compute the Cholesky factorization A = LL T , where L (referred to as the Cholesky factor) is a N × N real lower triangular matrix with positive diagonal elements. The solution vector x can then be computed by solving the two following triangular systems: Ly = b and L T x = y.

Algorithm 1: Pseudocode of the tiled Cholesky factorization

for k = 0 to n -1 do A[k][k] ← POTRF(A[k][k]); for i = k + 1 to n -1 do A[i][k] ← TRSM(A[k][k], A[i][k]); for j = k + 1 to n -1 do A[j][j] ← SYRK(A[j][k], A[j][j]); for i = j + 1 to n -1 do A[i][j] ← GEMM(A[i][k], A[j][k], A[i][j]);
To take advantage of modern highly parallel architectures, state-of-the-art numerical algebra libraries implement tiled Cholesky factorizations. The matrix A = (A ij ) 0≤i,j≤n is divided into n × n tiles (or blocks) of n b × n b elements, and the tiled Cholesky algorithm can then be seen as a sequence of tasks that operate on small portions of the matrix. This approach greatly improves the parallelism of the algorithm and mostly involves BLAS3 kernels whose library implementations are really fast on modern architectures. The benefits of such an approach on parallel multicore systems have already been discussed in the past [START_REF] Gustavson | High-performance linear algebra algorithms using new generalized data structures for matrices[END_REF], [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF], [START_REF] Quintana-Ortí | Programming matrix algorithms-byblocks for thread-level parallelism[END_REF]. Following the BLAS and LAPACK terminology, the tiled algorithm for Cholesky factorization is based on the following set of kernel subroutines: Note that no extra memory area is needed to store the L ij tiles since they can overwrite the corresponding A ij tiles from the original matrix. The tiled Cholesky algorithm can be written as in Algorithm 1.

In this sequential pseudocode, we can notice that some kernel subroutines depend on each other, while others can be processed in parallel. Such an algorithm is commonly represented by its task graph (or DAG) that depicts its actual dependencies. In this well established model, each vertex of the graph represents a call to one of the four kernel subroutines presented above. The edges between two task vertices represent a direct data dependency between tasks. Figure 1 depicts the task graph for the Cholesky decomposition of a 5 × 5 tiled matrix.

B. Multiprocessor scheduling 1) Static task allocation:

It is well known that the allocation of the tasks to the computing cores affect the performance and scalability, because of data locality and task heterogeneity. This problem has been adressed in the distributed memory context. For example, the ScaLAPACK library [START_REF] Dhillon | Lapack working note 95 scalapack: A portable linear algebra library for distributed memory computersdesign issues and performance[END_REF] first distributes the matrix tiles to the processors, using a standard 2D block-cyclic distribution of tiles along a virtual p-by-q homogeneous grid. In this layout the p-by-q top-left tiles of the matrix are topologically mapped onto the processor grid and the rest of the tiles are distributed onto the processors in a round-robin manner. It then implements an owner-compute strategy for task allocation: a task overwriting a tile is executed on the processor hosting this tile. This layout is also incorporated in the High Performance Fortran standard [START_REF] Koelbel | The High performance Fortran handbook, ser. Scientific and engineering computation[END_REF]. It ensures a good load and memory usage balancing for homogeneous computing resources [START_REF] Dhillon | Lapack working note 95 scalapack: A portable linear algebra library for distributed memory computersdesign issues and performance[END_REF]. However, for heterogeneous resources, this layout is no longer an option, and dynamic scheduling is a widespread practice.

These ideas also make sense in a shared-memory environment in order to take advantage of data locality. For instance the PLASMA [START_REF] Buttari | Lapack working note 191: A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] library provides an option for relying on such static schedules on multicore chips.

2) Dynamic task graph scheduling: Dynamic strategies have been developed in order to design methods flexible enough to cope with unpredictable performance of resources, especially in the context of real time systems, where on-line and adaptive scheduling strategies are required [START_REF] Chetto | Dynamic scheduling of real-time tasks under precedence constraints[END_REF], [START_REF] Manimaran | An efficient dynamic scheduling algorithm for multiprocessor real-time systems[END_REF]. More recently, the design of dynamic schedulers received a lot of attention, since on modern heterogeneous and possibly shared systems, the actual prediction of either execution and communication times is very hard, thus justifying the design of ad-hoc tools that will be described in Section IV.

As presented earlier, many scheduling heuristics have been proposed for DAGs since this problem is NP-complete. Most of these heuristics are list-scheduling heuristics: they sort tasks according to some criterion and then schedule them greedily. This makes them good candidates to be turned into dynamic scheduling heuristics. The best-known listscheduling heuristic for DAGs on heterogeneous platforms is certainly HEFT [START_REF] Topcuoglu | Performanceeffective and low-complexity task scheduling for heterogeneous computing[END_REF]. It consists in sorting tasks by decreasing bottom-level, which is the weight of the longest path from a task to an exit task (a task without successors). In a heterogeneous environment, the weight of a task (or communication) is computed as the average computation (or communication) time over the whole platform. Then, each task is considered and scheduled on the resource on which it will finish the earliest. HEFT turns out to be an efficient heuristic for heterogeneous processors. Other approaches have been proposed to avoid data movement when taking communications into account, such as clustering tasks into larger granularity tasks before scheduling them [START_REF] Sarkar | Partitioning and scheduling parallel programs for multiprocessing[END_REF].

III. MAKESPAN LOWER BOUNDS

Performance results for linear algebra computation are often accompanied with an upper bound in terms of FLOP/s (FLoating-point Operations Per Second), in order to assess the achieved efficiency. Since the theoretical peak performance is usually unreachable, particularly with GPUs, the common bound being used is the performance of a simple matrix multiplication (GEMM) since this is the most efficient dense linear algebra operation, and thus providing a good hint of some achievable performance. This bound takes into account the heterogeneity of the platform by summing up the obtained GFLOP/s (GigaFLOP/s) on the various processing elements. It however does not take into account the heterogeneity of the application, which is particularly important for small and medium matrices, for which a fair amount of the tasks are not GEMMs but much less efficient tasks such as POTRFs, especially on accelerators.

We here propose much more accurate bounds that take into account both heterogeneity of the computation resources and of the application kernels, by taking as input the execution time of any kernel on any type of resource. They also to a certain extent take into account the task graph itself, in terms of task dependencies.

A. Linear Programming formulation

The makespan lower bound computation is based on a relaxation of the scheduling problem, in which almost all precedence constraints are ignored. This formulation focuses on the number of tasks n rt of each type t (GEMM, SYRK, TRSM, POTRF) which are executed on each resource type r (CPU, GPU, ...). From the Cholesky task graph, we know the number N t of tasks of each type t that need to be performed on the whole platform, and from the platform we know the number M r of processing elements of each type r available to schedule the tasks. For each task type t and resource type r, the calibration mechanisms inside StarPU (described in Section IV-A) provide the execution time T rt of these tasks on this resource type. The basic area bound is obtained by solving the following linear problem: minimize the makespan l such that ∀t, all N t tasks of type t get executed over the various processing element types r:

r n rt = N t
∀r, the M r resources of type r complete all their tasks of various types t within the makespan l:

t n rt T rt ≤ l × M r ∀r, t n rt ∈ N +
It is clear that the optimal value l * of this linear program is a lower bound on the total execution time of the task graph, since any execution needs to execute all tasks. Ignoring the task graph precedences in this bound allows one to handle tasks of the same type with a couple of variables (one per resource type), instead of having one variable for each task in the graph, thus limiting the number of variables and reducing symmetries in the solution space. While being very naive, this formulation allows StarPU, without any input from the application beyond the normal task submission, to automatically generate it and solve it on the fly very quickly, right after the application execution, which thus allows one to print this theoretical bound along the measured performance in the application output.

Due to the actual timings of the different task types, this linear program always decides that all POTRF tasks should be executed on CPUs, since all other task types make much more efficient use of the GPU resources. However, in practice all POTRF tasks are on the critical path of the Cholesky graph, and hence this implies that the resulting lower bound is too optimistic for small matrix sizes, since it does not take dependencies into account. This interesting feature of the Cholesky task graph to contain a path with all n POTRF tasks can be used to strengthen the bound, without adding other variables in the linear program. In addition to the n POTRF tasks, this path contains n -1 of the n×(n-1)

2

TRSM tasks, and n -1 of the n×(n-1) 2 SYRK tasks. We can thus add the following constraint, which states that the execution time is necessarily larger than the time to execute all these tasks in sequence:

r n rP T rP + (n -1) × T * T + (n -1) × T * S ≤ l
In this constraint, T rP denotes the execution time of POTRF tasks on resource type r, and T * T and T * S denote the fastest execution time of TRSM and SYRK tasks: we do not model exactly on which resources these TRSM and SYRK tasks are executed, and thus underestimate their completion times, ignoring which resource they actually run on 1 . The resulting lower bound is called the mixed bound in the rest of the paper. This linear program has a very small number of variables and constraints (in particular, they are independent of the matrix size), and it can thus be solved very quickly.

B. Constraint Programming formulation

In addition to this lower bound computation, we have used a Constraint Programming formulation of the scheduling problem, in order to obtain good feasible solutions. These solutions provide both a comparison point for StarPU schedules and a limit for possible improvements of the lower bound. The formulation contains one boolean variable b ir for each task i and each resource type r (only one can be true for a given task), and one integer variable s i for each task i which represents the starting time of the task. The constraints are the following: minimize l such that ∀i, only one type of resource executes task i:

OnlyOne(b i1 , . . . , b iR ) ∀i, task i completes:

s i + r b ir T ir ≤ l
∀r, ∀t, at time θ the M r resources of type r are executing at most M r tasks:

|{i st s i ≤ θ < s i + r b ir T ir }| ≤ M r ∀i → j, dependency i → j is respected: s i + r b ir T ir ≤ s j
We have implemented this constraint programming formulation using CP Optimizer v12.4. The first constraint is expressed using the alternative constraint, and the third constraint uses the concept of cumulative functions to express the number of tasks which use resources of type r at time t. The other constraints are simple linear constraints and are easily expressed. The solver explores the solution space with an exhaustive search and backtracking, using constraint propagation to reduce the search space as much as possible.

Furthermore, providing the result of a HEFT heuristic as an initial solution allows the solver to explore good solutions more rapidly. We let the solver optimize for 23 hours and keep the best solution found in this duration. The obtained solutions are quite good compared to what is obtained with other heuristics, but the solver is unable to prove optimality.

Because it would otherwise be extremely costly to solve2 , this formulation does not take into account data transfers. With the usual platforms and the dense linear algebra operation being studied (the Cholesky factorization), data transfers are indeed not a concern: computation is dense enough for transfers to be largely overlapped with kernel computation.

C. Upper bounds on performance

Lower bounds on execution time also give upper bounds on the performance. Therefore, we have plotted different theoretical performance upper bounds of the Cholesky factorization in Figure 2, based on real execution timings of different tasks on the Mirage machine (described in section V-B). The critical path bound is calculated based on the critical path of the Cholesky task graph. While calculating the critical path, we have taken into consideration the fastest execution time of each task among the different resources. The area bound and mixed bound calculations are based on the description given in mixed bound subsection III-A. Since GEMM is the fastest kernel of the Cholesky factorization algorithm, we have also plotted the GEMM Peak. This plot shows that the mixed bound is the tightest upper bound among all upper bounds, and we will therefore compare the performance of our experiments only with the mixed bound in the experiment section.

The performance of the constraint programming solution (best solution found in 23 hours, but not a bound because CP is unable to prove its optimality in 23 hours for matrices larger than 5 × 5 tiles) described in section III-B will be discussed in section V-C3.

IV. TOOLS AND LIBRARIES

For this study, we used the Chameleon [START_REF]Chameleon, a dense linear algebra software for heterogeneous architectures[END_REF] implementation of the Cholesky factorization, running on top of the StarPU runtime system. We performed real executions on the target platform, and we additionally used the Simgrid [START_REF] Casanova | SimGrid: a Generic Framework for Large-Scale Distributed Experiments[END_REF], [START_REF] Stanisic | Modeling and Simulation of a Dynamic Task-Based Runtime System for Heterogeneous Multi-Core Architectures[END_REF] simulator, in order to reduce the experimentation time, improve reproducibility of the experiments, and also be able to modify the execution platform.

A. StarPU runtime system

StarPU [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF] is a runtime system aiming to allow programmers to exploit the computing power of the available CPUs and GPUs, while relieving them from the need to specifically adapt their programs to the target machine and processing units. The StarPU runtime supports a taskbased programming model. Applications submit computational tasks, forming a task graph, with CPU and/or GPU implementations, and StarPU schedules these tasks and associated data transfers on available CPUs and GPUs. The data that a task manipulates is automatically transferred between the local memory of the accelerators and the main memory, so that application programmers are freed from the scheduling issues and technical details associated with these transfers. In particular, StarPU takes care of scheduling tasks efficiently, using well-known generic dynamic and task graphs scheduling policies from the literature (see Section II-B), and optimizing data transfers using prefetching and overlapping, in particular. In addition, it allows scheduling experts, such as compiler or computational library developers, to implement custom scheduling policies in a portable fashion.

In this study, we specialize the StarPU scheduling algorithms to include a mixture of static and dynamic task assignments, based on the knowledge of the Cholesky task graph, to improve performance on small and medium size matrices. In the following, we call "small" a matrix with less than 10 × 10 tiles, "medium" a matrix with tile size between 10 and 20, and "large" a matrix with more than 20 × 20 tiles.

B. Chameleon dense linear algebra library

To cope with the increased degree of parallelism, a new class of dense linear algebra algorithms has been proposed, often referred as tile algorithms in the literature [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF], [START_REF] Quintana-Ortí | Programming matrix algorithms-byblocks for thread-level parallelism[END_REF].

These algorithms led to the design of new libraries in the past five years such as PLASMA [START_REF] Buttari | Lapack working note 191: A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF], FLAME [START_REF] Igual | The FLAME approach: From dense linear algebra algorithms to highperformance multi-accelerator implementations[END_REF] and DPLASMA [START_REF] Bosilca | Distributed dense numerical linear algebra algorithms on massively parallel architectures: DPLASMA[END_REF]. Although both static and dynamic versions of the algorithms have been initially implemented, the dynamic codes are now predominant since they proved to provide more flexibility. These dynamic codes rely on runtime systems (QUARK [START_REF] Yarkhan | QUARK Users' Guide: QUeueing And Runtime for Kernels[END_REF], Supermatrix [START_REF] Chan | SuperMatrix: a multithreaded runtime scheduling system for algorithms-byblocks[END_REF], PaRSEC [START_REF] Bosilca | PaRSEC: A programming paradigm exploiting heterogeneity for enhancing scalability[END_REF]) that have been specifically designed for the purpose of the numerical software (in the case of PLASMA, FLAME and DPLASMA, respectively).

The advantage of relying on specialized runtime systems is that they can be optimized for both the numerical algorithm and the target architecture. On the other hand, designing and maintaining a runtime system is a highly time consuming task, which makes it difficult to design a fullyfeatured specialized runtime system. The Chameleon library is based on the PLASMA tiled algorithms and code but relies on the StarPU generic runtime system instead of the specialized QUARK runtime system. One advantage is that it allows for handling heterogeneous architectures (whereas PLASMA and QUARK were initially designed for multicore chips). Another advantage when aiming at focusing on the impact of scheduling strategies is that it allows for running in simulation mode with the field-proven combination [START_REF] Stanisic | Modeling and Simulation of a Dynamic Task-Based Runtime System for Heterogeneous Multi-Core Architectures[END_REF] of StarPU and Simgrid.

C. Simgrid simulation engine

Simgrid [START_REF] Casanova | SimGrid: a Generic Framework for Large-Scale Distributed Experiments[END_REF] is a versatile simulation toolkit initially designed to study the behavior of large-scale distributed systems like grids, clouds, or peer-to-peer systems. It builds on fluid network models that have been proven as a reasonable alternative to both simple analytic models and expensive, difficult-to-instantiate packet-level simulations.

The Simgrid version of StarPU [START_REF] Stanisic | Modeling and Simulation of a Dynamic Task-Based Runtime System for Heterogeneous Multi-Core Architectures[END_REF] uses Simgrid to simulate the execution of an application within a single machine. The idea is to run the application normally, except that data transfers and computation kernel calls are replaced by a simple procedure accounting for the time they are expected to take, and gathered coherently by Simgrid. StarPU models each execution unit (CPUs and GPUs) by defining the time taken by each execution unit on each possible task/kernel [START_REF] Augonnet | Automatic Calibration of Performance Models on Heterogeneous Multicore Architectures[END_REF]. It also models the PCI buses between them, using offline bus bandwidth measurements, and relies on Simgrid to compute the interferences on PCI buses between the different transfers.

The resulting simulated times are very close to actual measurements on the real platforms [START_REF] Stanisic | Modeling and Simulation of a Dynamic Task-Based Runtime System for Heterogeneous Multi-Core Architectures[END_REF], and properly reproduce the various behaviors that can be observed for the various schedulers. This allows one to confidently run experiments with the Simgrid version of StarPU, which provides several advantages:

• The time to simulate execution is reduced, since no actual computation or data transfer is done. The Simgrid simulator itself is not parallel, so the whole execution gets serialized, but several simulations can be run in parallel for e.g. various matrix sizes or schedulers, and one then gets all the results in parallel. • The experiments do not depend on the availability of the platform, both in terms of quotas, and in terms of versions of the installed software, thus allowing reproducible experiments. This proved useful while performing the experimentation for this very article, since the platform became unavailable for a couple of weeks due to Air Conditioning issues. • The platform can be modified, for instance to change the available PCI bandwidth, the execution times of the kernels, etc. In Section V-C2, we use this feature in order to build a virtual "related" heterogeneous platform. In "actual execution mode", we perform the real execution on Mirage machine (described in section V-B) with StarPU runtime system. While in "simulation mode", we perform simulation on any machine with Simgrid version of StarPU runtime system by using the configuration files of target platform and expected execution time of kernels on each resource of the target platform.

V. EXPERIMENTS AND RESULTS

A. Schedulers

We have experimented with a few schedulers of StarPU, which are representative of state-of-the-art dynamic heuristics.

The random scheduler assigns tasks randomly over all the computation resources. It uses an estimation of the relative performance of the resources as coefficients to balance the randomness, so that GPUs will be assigned more tasks, according to their average acceleration ratio. This is thus representative of classical partitioning heuristics, which take into account the heterogeneity of the platform, but do not take into account the heterogeneity of the tasks.

The dmda (deque model data aware) and dmdas (deque model data aware sorted) schedulers use the minimum completion time heuristic to assign tasks to computational resources: each task is assigned to the processing resource which is estimated to complete it first, taking into account both the estimated computation time on the estimated target resource, and the possible required data transfer time. The difference between dmda and dmdas is that dmdas schedules tasks in order of their priorities, thus making it representative of the state-of-the-art HEFT heuristic [START_REF] Topcuoglu | Performanceeffective and low-complexity task scheduling for heterogeneous computing[END_REF], [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF].

We are calculating the priorities of different tasks in dmdas by estimating the longest path (in terms of execution time) from a task to an exit task (a task without successors) in Cholesky task graph. For longest path calculation, we have taken the fastest execution time of each task among the different resources into consideration.

The Cholesky factorization is a structured application, so we can estimate some extra information in advance by analyzing the task graph with the help of different tools. This information could be an exact schedule, priorities for some specific tasks, scheduling of some tasks on a particular worker/resource type, etc. In the following section, we inject more or less of these extra information as static knowledge, to influence the scheduling decisions and get better performance.

B. Experimental Setup

We have used a machine called Mirage to run and simulate our experiments. It has 2 Hexa-core Westmere Intel R Xeon R X5650 processors and 3 Nvidia Tesla M2070 GPUs. In the actual execution, we used only 9 CPU cores of the mirage machine so that the remaining 3 CPU cores can be used to fully exploit the critical resource (GPUs) of the system. To make the performance comparable we stick to 9 CPU cores in all of our experiments.

We have used Chameleon v1.0, StarPU v1.2.0 and Simgrid v3.10 for our experiments. We used Intel R Math Kernel Library 11.1, MAGMA 1.4.1 and CUBLAS 6.0 to perform actual executions.

C. Results

We have divided our experiments into two categories based on the types of configuration used. The first one is Homogeneous category where we have run and simulated the performance behavior with 9 homogeneous CPU cores and the second one is Heterogeneous category, where we used 9 CPU cores and 3 GPUs to run the tasks.

From previous work we are getting maximum performance in heterogeneous case with tile size equal to 960 [START_REF] Agullo | Faster, Cheaper, Better -a Hybridization Methodology to Develop Linear Algebra Software for GPUs[END_REF], [START_REF] Agullo | QR factorization on a multicore node enhanced with multiple GPU accelerators[END_REF], that is why we also kept the same tile size value throughout all our experiments.

For actual executions, we provide the average and standard deviation of 10 runs in the plots. In simulation mode, results are deterministic for all schedulers except for the random scheduler which relies on random allocation choices. The simulated plots therefore provide average and standard deviation values of 10 simulations with various seeds for the random scheduler.

1) Homogeneous case: For the homogeneous case, we provide the results of real execution runs of Cholesky factorization with the three different StarPU schedulers: random, dmda and dmdas.

From Figure 3, it is clear that the random scheduler does not perform well. This happens because it does not take into account the already assigned workload of the workers, and just selects a worker among all workers with equal probability. This shows that the scheduler needs to take scheduling decisions in some smart way. The other two schedulers which are based on data aware and early finish time strategies perform much better than the random scheduler. Figure 3 also shows that dmdas slightly underperforms compared to dmda for smaller number of tiles. This is due to the fact that dmdas is biased towards the longest path (path with more work) and chooses some tasks in the beginning which do not generate enough level of parallelism. But as time progresses, dmdas starts choosing tasks which releases a higher number of tasks, because these tasks would be the critical ones, which improves the overall performance of the execution.

We are also interested to know what the upper bound of the performance is, in order to determine how far these results are from that bound with different types of schedulers. Since actual executions add some runtime overhead and affect the performance, to mitigate this overhead we have compared the bound with simulated performance. Figure 4 shows that the behavior is very similar to the original execution, with a slight increase in performance, since we have removed the runtime overhead from the simulation. It also shows that the gap between mixed bound and achieved performance is significant for small matrices.

2) Heterogeneous Case: In this subsection, we consider all the processing units of the Mirage machine. 9 CPUs and 3 GPUs are used for the execution of tasks while the remaining 3 CPUs are used as drivers for the 3 GPUs. Table I shows the GPUs performance for each kernel with respect to CPUs performance, e.g.: GEMM is 29 times faster on GPU compared to CPU.

We divide our work into two parts. In the first part, we consider the impact of heterogeneity of resources by considering a heterogeneous platform with related performance. More specifically, we designed a fictitious hardware configuration, where execution time of each kernel on GPU is made to be exactly K times faster than the CPU execution time, and we call this case the heterogeneous related. The common accelerator factor K is an average over the actual measured acceleration factors, computed as follows :

K = N P * a P + N T * a T + N S * a S + N G * a G Total Number of Tasks
where, N P : total number of POTRF tasks a P : acceleration factor of POTRF on GPU N T : total number of TRSM tasks a T : acceleration factor of TRSM on GPU N S : total number of SYRK tasks a S : acceleration factor of SYRK on GPU N G : total number of GEMM tasks a G : acceleration factor of GEMM on GPU Here, the acceleration factor depends on the number of tasks and the number of tasks depends on the number of tiles. In the second part of our work, we show the achieved performance with the actual hardware with the help of both actual and simulated executions, and we call this case the heterogeneous unrelated case.

We are using the mixed bound (as explained in Section III-A) to compare the performance. The bounds do not take into account the communication constraints. Therefore, to be fair in the comparison we have used the simulated performance, where communication costs have been removed by modifying the platform file of our machine (one of the interesting features of the Simgrid version of the StarPU runtime system). Heterogeneous related case: Figure 5 shows the simulated performance with different schedulers on the fictitious heterogeneous platform. Here, we can see that the random scheduler performs very poorly because it assigns tasks randomly to the worker without knowing the already assigned workload of workers, which limits the number of ready tasks in the system, and introduces significant idle time on our critical resource (GPUs). We have also computed the mixed bound for this fictitious platform. The difference between simulated performance and mixed bound is once again significant for small and medium size matrices.

Heterogeneous unrelated case: First we compare the performance of different schedulers in actual execution and then between simulated performance and mixed bound. As shown in Figure 6, in actual executions, the random scheduler does not perform well because it is not taking data movement into account while making scheduling decisions: it assigns worker randomly for each task, which may select different resource types for data dependent tasks and result in lots of data movement from CPU memory to GPU memory and vice-versa. In addition, it is also not taking the affinity of tasks to resource (e.g.: GEMM/SYRK is more suitable to be executed on GPU) into account, which degrades the overall performance of the system. The other two schedulers perform comparatively better than the random scheduler because they take into account data transfers when assessing completion time in the HEFT-like scheduling strategies.

Here we can also see that dmda outperforms dmdas performance for some matrices, for the same reason as for the homogeneous case (choosing the critical task versus tasks which generate high level of parallelism). We are now again interested in determining how far we are from the peak performance of the application. Thus, we performed the simulation with different numbers of tiles. Figure 7 illustrates the comparison between bounds and achieved performance in simulation. Here we can also see that the performance difference between the best scheduler and the mixed bound is significant for small and medium size matrices.

Comparison between Heterogeneous related and unrelated case: In order to determine the impact of heterogeneity of speed-up of tasks on performance, we present a comparison between related and unrelated heterogeneous simulations. To this end, we scaled the mixed bound of the related case such that it perfectly matches with the mixed bound of the unrelated case, and also scaled all the performance values of the related case with the same factor. The obtained results are given in Figure 8, which can now be compared with the unrelated case of Figure 7.

Here we can see that unrelated speed-ups make the problem harder. That is why the gap between state-of-the-art schedulers performance and mixed bound is large in Figure 7 compared to Figure 8. Here, it is also clear that there is room for improvement in the case of small and medium size matrices in the heterogeneous case.

3) Scheduling with static knowledge: The significant gap between the performance of StarPU schedulers and the theoretical bound (mixed bound) for small and medium size matrices in Figure 7 highlights the following things:

• Either the dynamic schedulers of StarPU return a scheduling that can be improved for small and medium size matrices; • Or the theoretical bound is not tight enough; • Or both. Indeed, the dmda and dmdas schedulers take only dynamic decisions to map the ready tasks onto the processors depending on the state of resources and estimation of execution and communication times (also priorities among ready tasks in dmdas), without taking into account the overall task graph. These local choices may lead to bad decisions when the parallelism in the task graph is limited. We conducted some experiments to improve the overall performance with static information in the heterogeneous unrelated case.

Since GEMM and SYRK kernels are well suited to execute on GPUs, we enforced these kernels to be executed on GPUs as static information to the StarPU runtime system. This strategy improves the performance slightly for some matrices in simulation but the performance improvement was not significant and the reason for this is that the StarPU schedulers (dmda and dmdas) already choose GPUs to execute most of the GEMM and SYRK kernels.

We also analyzed the solution of the mixed bound and noticed that a significant portion of the TRSM kernels were mapped onto CPUs. Analyzing traces generated by dmda and dmdas schedulers reveals that both policies allocate very few TRSMs on CPUs. Since the mixed bound does not take all dependencies into account, it is not clear which TRSM kernels should be executed on CPUs in order to improve the performance. On the Mirage machine, with real timings of tasks, we found that the critical path of the Cholesky factorization passes through the diagonal and second diagonal tiles (sequence of POTRF → TRSM → SYRK → POTRF ..... → SYRK → POTRF). Therefore, we have evaluated the performance in simulation with dmdas scheduler where all the TRSM kernels which are at least k (1 ≤ k < Number of Tiles) tiles away from the diagonal are forced to execute on the CPUs (see Figure 9) and plotted the best obtained performance in Figure 10. We obtained best performance when all the TRSM kernels which are more than 6-8 tiles away from the diagonal are forced on CPUs. Figure 10 shows that providing information about the TRSMs triangular structure statically allows one to achieve better performance than present state-of-the-art schedulers for small and medium size matrices. We eventually used the constraint programming (CP) described in Section III-B to find an optimal solution and ran for 23 hours, but unfortunately we did not manage to get an optimal solution, particularly for large matrix sizes, which produce a very large constraint program (and thus this is not a performance bound). Nevertheless, for reasonable matrix sizes, it provides good and feasible solutions in that span of time. Theoretical performance value with CP solution (CP solution(23 hrs) in Figure 10) was better than the values what we are getting with state-of-the-art schedulers in simulation for small and medium size matrices. We thus injected the exact schedule obtained from CP solution in the simulation and obtained almost equal (difference is less than 1%) performance (CP solution in simulation in Figure 10) compared to theoretical performance value, which also shows the robustness of the Simgrid version of StarPU with simulation.

Performance improvement obtained in simulation by injecting static information to scheduler motivated us to conduct some actual execution with static information. Therefore, we conducted some actual execution by injecting the triangular information, force all TRSMs on CPUs which are at least k (1 ≤ k < Number of Tiles) tiles away from diagonal as static knowledge. Figure 11 shows the best obtained performance in actual execution among performance obtained with all possible values of k. We also conducted some experiments by injecting the CP schedule in actual execution for small matrices, however we did not achieve good performance improvement compared to what we are achieving in simulation. The CP formulation indeed does not account data transfers, since as described in Section III-B, solving a CP with data transfers has shown intractable for the purpose at stake. Actual execution with CP schedule thus adds lots of idle time on resources during data transfer, and consequently does not reproduce the same performance in actual execution. The simulated execution has however allowed us to show, at least in the case without data transfers, that some heuristics get relatively close to an achievable CP solution. We are currently extending the CP formulation to partially take data transfers into account, so that it can be used for real executions, but this is beyond the scope of this paper.

VI. DISCUSSION

A. dmda vs dmdas scheduler

We were expecting that dmdas would always perform better than dmda scheduler because it is also taking the HEFT priorities into account while making scheduling decision. Nevertheless, we found a few cases where dmda outperforms dmdas. We investigated the generated trace files with dmda and dmdas schedulers in order to determine the reasons of this behavior and we found that dmdas puts emphasis on critical path rather than parallelism, since it selects some tasks in the beginning which are critical but are not generating enough level of parallelism. That introduces some idle time on the critical resource (GPUs) and degrades the overall performance of the system, which is a known defect of the HEFT scheduler in general. Figure 12 shows traces with dmda and dmdas schedulers. We conducted some experiments in simulation by injecting only the mapping information (i.e. only the CPU/GPU information, not the exact task order) of the feasible solution statically obtained by constraint programming, and let the scheduler decide the precise ordering and worker dynamically. This extra information about resource allocation did not improve the performance of the system compared to the performance obtained by dmda and dmdas schedulers, which indicates that the feasible solution is highly dependent of the precise ordering chosen by constraint programming. This shows that heuristics required to achieve this performance are probably very complex, probably even beyond only backfilling.

VII. CONCLUSION

In this work, we have bridged the gap between theoretical performance bounds and actually achieved performance on the dense Cholesky factorization. On the former side, we have proposed improved bounds which take into account both resource and task heterogeneity, as well as critical paths. On the latter side, we have introduced some static information into the dynamic task scheduler of StarPU, which brought the performance closer to the theoretical bounds, and very close to what a statically-optimized schedule can achieve. We have also shown that the performance achieved by such statically-optimized schedule depends on precise non-intuitive task ordering, which thus can not be reached by simple list-scheduling heuristics, even with backfilling.

We will verify the results on other hardware platforms, and apply the same methodology to other dense linear algebra algorithms, but we also plan to try other classes of applications, notably more irregular applications such as sparse linear algebra or FMM (Fast Multipole Method).

More generally, this work opens a bridge to closer interaction between applications and tasks schedulers. We have shown that while generic heuristics such as HEFT achieve very good performance, application-specific scheduling hints can noticeably improve performance. We aim at generalizing and formalizing this type of information, so that scheduling experts can easily analyze achieved performance, optimize the schedule statically, and try to inject more or less application-specific scheduling hints into the scheduler, such as "this proportion of TRSM tasks should be run on CPUs", or "these TRSM tasks should be run on CPUs", etc. The code produced for the purpose of the study will be reversed in the StarPU runtime system and Chameleon dense linear algebra library.

•

  POTRF: This LAPACK subroutine is used to perform the Cholesky factorization of a symmetric positive definite tile A kk of size n b × n b producing a lower triangular tile L kk of size n b × n b .• TRSM: This BLAS subroutine is used to apply the transformation computed by POTRF to a A ik tile by means of a triangular system solving. • GEMM: This BLAS subroutine computes a matrixmatrix multiplication of two tiles A ik , A jk and subtract the result to the tile A ij . The old value of A ij is overwritten by the new one.• SYRK: This BLAS subroutine executes a symmetric rank-k update on a diagonal tile A kk .

Figure 1 .

 1 Figure 1. Task graph for the Cholesky decomposition of a 5 × 5 tiled matrix.

Figure 2 .

 2 Figure 2. Heterogeneous theoretical performance upper bounds

Figure 3 .

 3 Figure 3. Homogeneous actual performance

Figure 4 .

 4 Figure 4. Homogeneous simulated performance

Figure 5 .

 5 Figure 5. Heterogeneous related simulated performance

Figure 6 .

 6 Figure 6. Heterogeneous unrelated actual performance

Figure 7 .

 7 Figure 7. Heterogeneous unrelated simulated performance

Figure 8 .

 8 Figure 8. Heterogeneous related simulated scaled performance

Figure 9 .

 9 Figure 9. TRSMs forced on CPUs

Figure 10 .

 10 Figure 10. Heterogeneous unrelated simulated performance with static knowledge

Figure 11 .

 11 Figure 11. Heterogeneous actual performance with static knowledge

Figure 12 .

 12 Figure 12. GPU Traces for 8 × 8 tiles

  Therefore, we get different acceleration factors with different number of tiles. Acceleration factors for 4, 8, 12, 16, 20, 24, 28 and 32 tiles matrices are 17.30, 22.30, 24.30, 25.38, 26.06, 26.52, 26.86 and 27.11 respectively.

It is possible to include additional variables to the linear program to have more precise values, but this does not provide a better bound unless we take more dependencies into account, which requires adding too many variables and constraints and makes the linear program intractable.

We also have written a version of the constraint programming formulation which takes data transfer times into account but we could not obtain results at the scale of interest for this paper.

ACKNOWLEDGMENT Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, being developed under the INRIA PlaFRIM development action with support from LaBRI and IMB and other entities: Conseil Régional d'Aquitaine, Université de Bordeaux and CNRS. This work has been partially supported by the ANR SOLHAR project, funded by the French Research Agency.