Computational advances in gravitational microlensing: a comparison of CPU, GPU, and parallel, large data codes

To assess how future progress in gravitational microlensing computation at high optical depth will rely on both hardware and software solutions, we compare a direct inverse ray-shooting code implemented on a graphics processing unit (GPU) with both a widely-used hierarchical tree code on a single-core CPU, and a recent implementation of a parallel tree code suitable for a CPU-based cluster supercomputer. We examine the accuracy of the tree codes through comparison with a direct code over a much wider range of parameter space than has been feasible before. We demonstrate that all three codes present comparable accuracy, and choice of approach depends on considerations relating to the scale and nature of the microlensing problem under investigation. On current hardware, there is little difference in the processing speed of the single-core CPU tree code and the GPU direct code, however the recent plateau in single-core CPU speeds means the existing tree code is no longer able to take advantage of Moore's law-like increases in processing speed. Instead, we anticipate a rapid increase in GPU capabilities in the next few years, which is advantageous to the direct code. We suggest that progress in other areas of astrophysical computation may benefit from a transition to GPUs through the use of"brute force"algorithms, rather than attempting to port the current best solution directly to a GPU language -- for certain classes of problems, the simple implementation on GPUs may already be no worse than an optimised single-core CPU version.


Introduction
Gravitational microlensing is the deflection of light by matter in the regime where multiple-imaging and high magnification events occur, but the angular separation of images is too small to resolve. The observable effect is a change in brightness of a source over a period of time, with the details of the light curve depending on the complexity of the lens distribution and the surface brightness profile of the source. Gravitational microlensing has been used to study high magnification events on local scales, viz. the searches for dark matter compact objects in the Galactic bulge and halo (see for example references in and citations to Alcock et al. 1993;Aubourg et al. 1993;Udalski et al. 1993;Calchi Novati et al. 2002;Wyrzykowski et al. 2009) and extrasolar planets (e.g. Gould & Loeb 1992;Bond et al. 2004;Sumi et al. 2010), and quasar microlensing on cosmological scales (e.g. Vanderriest et al. 1989;Irwin et al. 1989;Kochanek 2004;Bate et al. 2008;Dai et al. 2010). For recent reviews of microlensing on all scales, see Wambsganss (2006), Kochanek et al. (2007), Gould (2009), Mao (2008) and Gaudi (2010). 6 Research undertaken as part of the Commonwealth Cosmology Initiative (CCI: www.thecci.org), an international collaboration supported by the Australian Research Council Email addresses: nbate@astro.swin.edu.au (N.F. Bate), cfluke@astro.swin.edu.au (C.J. Fluke) Microlensing on local scales occurs in the low optical depth regime, where the source is affected by at most a few microlenses along the line of sight. Large numbers of stars must therefore be monitored in order to observe a single microlensing event. Conversely, the microlensing optical depth at the location of lensed quasar images is likely to be ∼ 1. This microlensing signal is therefore more complex, as it is the result of a large number of microlenses acting on the source essentially all the time (Wambsganss, 2006). In the remainder of this work, we focus on high optical depth microlensing.
A key element in many microlensing calculations is the determination of the magnification map; a typical example is shown in Figure 1. This map provides an estimate of the pointsource brightening over a finite region of a two-dimensional source plane. The magnification of an extended source is found by convolving a source profile (intensity and geometry) with the magnification map. This approach has been used to probe the structure of quasar emission regions, including accretion discs (e.g. Kochanek 2004;Eigenbrod et al. 2008;Floyd et al. 2009) and broad emission line regions (e.g. Lewis & Ibata 2004;Abajas et al. 2007).
The limitation in many microlensing calculations is the computational overhead in producing magnification maps with sufficient resolution and high (statistical) accuracy. Fitting observed light-curves obtained over long periods of time requires large maps with high pixel resolution to resolve small sources (such as quasar X-ray emission regions). Similarly, fitting Figure 1: A sample microlensing magnification map for a microlensing model with the following parameters: total convergence σ = 0.4, shear γ = 0.0, N pix = 1024 2 pixels in the source plane and N rays = 2000 light rays per source pixel (on average). These parameters correspond to N * = 558 lenses. The source plane has a side length of 20 Einstein radii. The greyscale covers the range from low magnification (black) to high magnification (white). multi-wavelength observations requires both large maps and high pixel resolution in order to simultaneously model emission from small optical emission regions and larger infrared emission regions. The situation is even worse when attempts are made to simultaneously model microlensing of accretion discs and the much larger broad emission line regions, a difference in scale of 2-3 orders of magnitude. Finally, many maps must be generated in order to ensure statistical independence of simulated measurements. It is only reasonable to sample any one magnification map as many times as you can place a source on it without overlapping.
In this paper we compare three techniques for computing magnification maps: a hierarchical tree code for single-core CPUs; a parallel, large data tree code for use on a CPU-based cluster supercomputer; and a direct, inverse ray-shooting code for use on graphics processing units (GPUs). We consider the accuracy of the codes, to determine whether the approximations inherent in the tree code are valid over a wider range of parameter space than has been feasible to test before, and compare run-times between the single-core tree code and the GPU code over a range of astrophysically-motivated parameters. While the "brute force" approach to ray shooting has not been considered a feasible option on CPU architectures, we show that the GPU implementation is highly competitive, and is indeed faster than the tree code in many circumstances. The expected increase in GPU speeds in the years ahead bodes well for advances in gravitational microlensing simulations. Our work hints at the potential benefit to other fields of astrophysical com-putation from the implementation of simple, "brute force" algorithms as a means of accelerating the adoption of GPUs in astronomy for time-consuming computations.
The remainder of this paper is set out as follows. In Section 2, we describe the process of generating magnification maps, and summarise key features of three gravitational microlensing codes which can be used for this task. In Section 3, we present a series of code comparisons (accuracy and speed). In Section 4, we discuss the relative merits of hierarchical and direct rayshooting approaches, and suggest the classes of problems to which each is best-suited. We present our conclusions in Section 5.

Generating magnification maps
In general, magnification maps are calculated using the (twodimensional) gravitational lens equation which relates the locations of the source, y, and an image, x, via the deflection angle, α(x). The magnification is: where A = ∂y/∂x is the Jacobian matrix of equation (1). Equation (1) gives a one-to-one relationship between an image location and the corresponding source location (x → y), but a one-to-many relationship for the alternative case of mapping a source location to its multiple images (y → x). While analytic solutions exist for a small number of simple lens models (e.g. Schneider et al. 1992), in general a numerical approach is used to calculate µ over a finite region of the source plane. This can be achieved most easily using inverse ray-shooting (Kayser et al. 1986;Schneider & Weiss 1986;Schneider & Weiss 1987). Here, light rays are propagated backwards from the observer, through the lens plane, where they are deflected using equation (1), and then mapped to a grid of pixels in the source plane. This approach avoids the need to determine the location of every image. By counting the number of light rays reaching each source plane pixel, N i j , compared to the average per pixel if there was no lensing, N rays , an estimate of the per-pixel magnification, µ i j , is obtained: For cosmological microlensing problems, the appropriate model for the deflection angle describes N * compact lens masses in a smooth background mass distribution, acted on by an external (macro-model) shear γ. This model, an extension of the one proposed by Chang & Refsdal (1979), is suitable for considering the impact of lensing due to individual stars where the external shear is due to the large-scale galactic mass distribution. The lens equation in this case is: where the convergence σ = σ c + σ * combines contributions from smooth matter σ c , and compact objects σ * . This is the model for the deflection of light by a screen of mass that we will use throughout this paper. The relationship between σ * and the number of lenses N * is: If the deflecting mass distribution were smoothed out, a circle in the image plane would map onto an ellipse with major-to-minor axis ratio |1 − σ + γ|/|1 − σ − γ| in the lens plane. The area of this 'minimum' ellipse is set by the size of the image plane. In practice, the graininess of the matter distribution can cause rays to be scattered outside the receiving area. Similarly, rays from a long way outside the shooting region can be scattered into the receiving square in the source plane. For that reason, the lenses in the codes considered here are distributed in a circle of area A with a diameter larger than the long axis of the minimum ellipse.
The number of calculations, Φ, required to obtain a mangification map with N pix pixels scales as: where N FLOP is the number of floating point operations for each deflection calculation, and N rays is chosen to achieve a desired level of accuracy. It was quickly realised that the direct rayshooting approach was not practical on single CPUs: a scenario, Φ fiducial ∼ 3×10 16 , with N pix = 2500 2 pixels, N rays = 500, N * = 10 6 lenses, 1 and N FLOP = 10 would take months to years to calculate (Wambsganss, 1999) on hardware that was available at the time. Wambsganss (1990;1999) introduced a method based on a hierarchical tree code (Barnes & Hut, 1986). Here, the contribution of each lens depends on its distance from the light ray -lenses that are further away are grouped together into pseudo-lenses with a higher mass. This approach reduces N * , thus increasing the speed of the overall calculation, but at the cost of a slight error in the magnification map due to the approximation. This approach has been used to constrain both the size and temperature profile of quasar accretion discs (e.g. Wambsganss, Paczynski, & Schneider 1990;Rauch & Blandford 1991;Morgan et al. 2006;Pooley et al. 2007;Keeton et al. 2006;Anguita et al. 2008;Eigenbrod et al. 2008;Bate et al. 2008;Floyd et al. 2009), dark matter fractions in the lensing galaxies (e.g. Pooley et al. 2009;Dai et al. 2010), and possibly even the orientation of the quasar with respect to the line of sight (Poindexter & Kochanek, 2010).
Despite the much-improved speed compared to direct rayshooting, the hierarchical method is complex to implement, and requires a certain amount of expert knowledge. These facts make slight modifications or additions to the code quite challenging. Parameters such as the required accuracy (which determines whether to use a particular lens cell, or its four subcells, based on the distance between the ray and the cell's centre of mass) must be determined experimentally by repeatedly 1 Relevant for the case σ ∼ 1 running the code until acceptable results are obtained. Additionally, the generation of the tree data structure imposes a memory overhead which exceeds typical single CPU memory sizes for large N * . A parallel hierarchical tree code, suitable for running billion-lens calculations, has now been developed by Garsden & Lewis (2010). This approach overcomes the single-CPU memory limitation by distributing the computation over a number of nodes in a CPU-based cluster supercomputer.
Alternatives to the hierarchical method include the use of polygonal cells (Mediavilla et al., 2006), optimised to map areas of the image plane to the source plane thus reducing N rays required for a given accuracy, and the contouring method for point sources (Lewis et al. 1993;Witt 1993), which maps an infinite line plus N * closed loops from the image to source plane. The polygonal cell approach results in its own computational overhead in determining the correct polygonal lattice, and suffers from the need for a substantial number of polygons in the vicinity of caustics, which results in lengthy computation times. The contouring method has limited applicability to extended sources. Wyithe & Webster (1999) presented an extension of the contouring method to cover extended sources, however the additional computations required slowed the code to the point where it could no longer compete with the simplicity of inverse ray-shooting.
A recent alternative is the appearance of the GPU as a mathematical co-processor.
Originally developed to enhance the rate of generating graphics for the computer game industry, GPUs are revolutionising scientific computing (Fournier & Fussell 1988;Tomov et al. 2003;Venkatasubramanian 2003;Owens et al. 2005). Their highly parallel processing architecture, accessed through flexibile programming libraries such as NVIDIA's CUDA 2 or the crossplatform OpenCL 3 standard created by the Khronos Group, is enabling speed-ups of O(10) to O(100) times over a broad range of algorithms. Indeed, the notion of general purpose computing on graphics processing units (GPGPU) has motivated a reinvestigation of time-consuming processing tasks.
From a conceptual and implementational point of view, direct ray-shooting is a very simple algorithm with a high degree of parallelism: the deflection of each light ray is independent of every other light ray, and the deflection of a light ray due to a single lens is independent of all other lenses. The existence of such fine-grained parallelism means that ray-shooting can be considered an "embarassingly parallel" algorithm -ideal for GPGPU. Thompson et al. (2010) reported on a CUDA implementation of direct ray-shooting, and demonstrated that their code could calculate a magnification map for the Φ fiducial scenario in ∼ 7 hours using a four-GPU NVIDIA S1070 Tesla unit. Moreover, they were able to calculate a magnification map with one billion lenses (N pix = 512 2 and N rays = 42) in 24 hours, with a peak sustained processing rate of 1.28 Tflop/sec. While the CUDA implementation of direct ray-shooting gave an O(100) processing speed-up compared to using the same method on a single CPU, a critical comparison omitted from Thompson et al. (2010) was between the tree code and the GPU code.
In this paper, we compare: • The "industry standard" single-core tree code by Wambsganss (1990;1999), which we refer to as CPU-T; • The parallel, large-data tree code by Garsden & Lewis (2010), which we refer to as CPU-P; and • The CUDA direct ray-shooting code by Thompson et al. (2010), which we refer to as GPU-D.
We now describe salient points of these three approaches in more detail, including implementation issues that are relevant for the code comparison in Section 3.

CPU-T
The most computationally intensive part of the ray-shooting algorithm is the calculation of the deflection of individual rays by each lens. The hierarchical tree code attempts to minimise the number of calculations necessary by grouping together lenses that are sufficiently distant from a given light ray. A group of lenses is treated as a single pseudo-lens with a total mass equal to the sum of the individual lens masses, located at their centre of mass.
The first step in the calculation is to divide up the lenses into a cell structure, starting with the root cell, which is the size of the lens plane. Each cell that contains more than one lens is subdivided into four cells with half the side length of the parent cell. This process is continued iteratively until only cells containing zero or one lens remain. This calculation need only be carried out once per magnification map, and typically takes only a small fraction of the code's total run-time. A determination must be made for each light ray regarding which lenses need to be included individually and which can be grouped together into pseudo-lenses. The lens equation then becomes: where N L is the number of lenses to be included directly, and N PL is the number of pseudo-lenses with total mass m PL j and centre of mass position x PL j . We note that the actual implementation of the tree code uses higher order multipoles of the mass distribution in order to increase accuracy. This calculation scales on average as O(logN * ).
Additional techniques are used in the standard hierarchical tree code in order to improve the accuracy and speed of the calculation. For example, "test rays" are used to determine the cell structure, and then the same cells are used for a number of actual rays surrounding the test ray. Also, the deflection due to pseudo-lenses is only calculated directly for the test rays; for the actual rays, pseudo-lens deflections are interpolated between those calculated for the test rays. For full details, see Wambsganss (1999).

CPU-P
A tree code requires the following data structures to be stored in memory: 2-D magnification map array; lens array (mass and 2-D location for each lens); and the cell tree (stored as an array of indices), which provides a significant memory overhead as N * increases. In order to increase the scale of microlensing simulations beyond those that were achievable with CPU-T, where typical single-core memory permits a maximum of 20 million lenses, Garsden & Lewis (2010) implemented a distributed, parallel version of the tree code (CPU-P). They take advantage of the inherent parallelism of the ray deflections by distributing rays between multiprocessors of a CPU-based cluster supercomputer, using a simple but effective approach to static load balancing between computing nodes. To handle the large memory requirements for billion-lens simulations (∼ 200 GB), they use a combination of local memory, distribution of data between nodes (thus sharing the memory requirements), and disk-based storage (although this limits processing speeds somewhat by the requirement of costly file accesss).
The CPU-P code implements the same cell-tree and microlensing algorithms as CPU-T, but surrounded by a quite different environment within a running simulation. The generation of lens masses at the beginning of the run, and the ray-shooting itself, use as much parallelism as is available. The generation of the cell tree has not been parallelised, but is executed concurrently with another phase in the simulation. All these operations are logically the same as in the CPU-T version but execute without knowledge that they are running concurrently, or that their data may be contained within huge disk files. Disk files naturally degrade performance and parallel speed-up, but there are various caching schemes which can make up for that. CPU-P has been built to have the same accuracy as CPU-T, because the same numerical algorithms are used, however due to the reimplementation in a different programming language (FORTRAN for CPU-T, C++ for CPU-P), with different random number generators, and a different precision of integers (arbitrary precision has been implemented), there will naturally be minor differences.

GPU-D
The GPU parallel implementation of direct ray-shooting is achieved by splitting the total number of light rays into N B batches, which are deflected in parallel. For convenience, the deflection calculation of equation (4) is recast as: where T t=1 N t = N * for N T parallel threads, and each processing step considers only a single light ray. The contributions from γ and σ c are determined at the end of the deflection calculation.
The GPU-D code, described in Thompson et al. (2010), was implemented in CUDA according to the following algorithm: 1. N * lens coordinates are obtained (see below) and stored in CPU memory. Coordinates are loaded into GPU device memory as required; 2. N B = 2 17 light ray coordinates are pseudo-randomly generated on the CPU and loaded into GPU device memory; 3. GPU computation is initialised; computation is split into groups of N T = 128 threads; 4. Each thread group loads N T lenses and calculates deflection on N T rays; this is repeated until all lenses and rays are exhausted; 5. Once computation on the GPU is complete, results are copied back from the GPU device memory to system memory; 6. CPU maps the ray locations onto the source pixel grid in order to obtain the magnification map; 7. Steps 2-6 repeated until an average of N rays rays per source pixel is reached.
The only change in approach from that of Thompson et al. (2010) is the way the N * lenses coordinates are generated: for this work, we first run a tree code and output lens positions to a data file. These coordinates are then read from the data file by the GPU-D code, and subsequent processing is carried out as described above. Without this modification, we would not be able to make pixel-by-pixel comparisons of the magnification maps, and would be left with just statistical comparisons based on, for example, magnification probability distributions. For full details on choice of optimisation of GPU parameters, especially N B and N T , see Thompson et al. (2010), as these are somewhat hardware dependent.

Comparing the methods
In this section, we consider the relative accuracies of CPU-T, CPU-P and GPU-D over a much wider range of parameter space than has been considered previously, and we compare processing times between CPU-T and GPU-D, noting regions of parameter space where each code might currently be preferred. The former aspect is important in assessing the validity of all microlensing results that have relied on the tree code to date, the latter in motivating a choice of code for future applications. The choice of model parameters for comparison is somewhat arbitrary. We attempt to cover the majority of σ and γ parameter space typically probed by realistic lens models, as well as typical map sizes and pixel scales used in those analyses (see for example: Witt et al. 1995;Kochanek 2004;Morgan et al. 2006;Bate et al. 2008;Mediavilla et al. 2009).

Accuracy
Since the tree code uses an approximation to the lens locations, in some sense the direct approach produces a magnification map that is more accurate for a given configuration of lenses. An unanswered question is whether the hierarchical approach maintains sufficient accuracy as the number of calculations, Φ, is increased beyond approximately 2 × 10 14 . This restriction is largely historical, due to the impracticality of using the direct ray-shooting approach on a single-core CPU beyond this limit.
We conducted simulations for a range of parameters, summarised below: • N pix = 1024 2 , 2048 2 and 4096 2 ; • N rays = 100, 200, 500, 1000 and 2000; and These parameter choices allow us to probe the range: where we assume that N FLOP = 10. We choose all lenses to have the same lens mass as this does not influence the run-time.
The CPU-T/CPU-P accuracy parameter was set to 0.60 for all of our simulations, as this was found to reliably produce acceptable magnification maps. CPU-P maps were generated using a minimum of 2 compute nodes, to ensure that the parallelisation was not producing inaccurate results. There is a subtle difference between the way N rays is defined for CPU-T/CPU-P and GPU-D. In both codes, the same number of rays are shot across the entire shooting grid. The CPU-T/CPU-P rays are shot on a regular grid, thus in the absence of lensing exactly N rays are received at each CPU-T/CPU-P pixel. The GPU-D rays are shot from random positions in the shooting grid, and so on average N rays are received at each GPU-D pixel. This will introduce Poisson noise into the GPU-D magnification maps, however the total number of rays shot is so large (10 8 to 4 × 10 10 in our simulations) that any error introduced by shooting rays randomly will be insignificant.
We use two methods for assessing accuracy: the root mean square difference between corresponding magnification maps, and a Kolmogorov-Smirnov (K-S) test on probability distributions for change in magnitude ∆m. The K-S test results in a measurement of the significance with which we can reject the null hypothesis, that the two samples were drawn from the same distribution. Change in magnitude is calculated from the (per pixel) map magnifications using the following relationship: Throughout, we define "excellent agreement" as a root mean square difference between maps of less than 0.2 magnitudes (often much less than 0.2 magnitudes), and a K-S test probability that the null hypothesis is not rejected of greater than 99.7%.

No shear
To begin with, we set both the external shear, γ, and the smooth matter fraction, s = σ c /σ, to zero. This allows us to consider just the effect of compact objects, which are directly related to the number of lenses via equation (5). Later, we allow both γ and s to vary.
As indicative results, we present the calculated RMS values and K-S probabilities, P(KS ), for the N pix = 4096 2 case in Table 1. For σ ≥ 0.3, we find excellent agreement between the magnification maps from CPU-T and GPU-D when N rays ≥ 500. The σ = 0.1 and 0.2 cases differ slightly. Here, a low number of lenses (N * = 65 and 163 respectively) mean that there is very little variation in magnification across the maps. The K-S test probabilities are therefore strongly affected by small differences in magnification probability distributions. For N pix = 1024 2 and 2048 2 , N rays = 1000 are sufficient to obtain agreement between CPU-T and GPU-D. For N pix = 4096 2 even N rays = 2000 are not enough. In other words, a larger number of rays per pixel are required for agreement between CPU-T and GPU-D when σ is low (≤ 0.2).
We note also that for a given σ and N rays , the RMS error tends to increase with increasing N pix . For example, when σ = 0.4 and N rays = 500, we obtain the following RMS differences: 0.039 mag (N pix = 1024 2 ), 0.081 mag (N pix = 2048 2 ) and 0.168 mag (N pix = 4096 2 ). In other words, for a fixed accuracy we must increase N rays by a factor of ∼ 2 if we increase N pix by a factor of 4.
For the comparisons between CPU-P and GPU-D, we find the same results. In fact, the CPU-T and CPU-P maps are virtually indistinguishable down to the accuracy of our comparison. They therefore behave identically when compared with GPU-D.

Shear and convergence
We then fixed the convergence at σ = 0.45, and examined the impact of varying shear over the range 0.0 ≤ γ ≤ 1.0 in increments of ∆γ = 0.1. The number of pixels was fixed at N pix = 2048 2 . Once again, the smooth matter fraction s was set to zero. We find very similar results to the no-shear case discussed above: for N rays ≥ 500, the results of CPU-T and GPU-D are indistinguishable (RMS difference < 0.2 magnitudes and K-S probability > 99.7%), and similarly for comparisons between CPU-P and GPU-D.

Smooth matter fraction
Finally, we test the accuracy of CPU-T and CPU-P when the convergence is split into a continuously distributed component σ c and a compact stellar component σ * . We set γ = 0.5 and σ = σ c + σ * = 0.45, and vary s in the range 0.0 ≤ s ≤ 0.9 in increments of ∆s = 0.1. Since there is no microlensing if s = 1, we also simulated the s = 0.99 case. For convenience, we choose N pix = 2048 2 . Table 1: Comparison between the CPU-T, CPU-P and GPU-D codes for magnification maps with N pix = 4096 2 , γ = 0.0 and s = 0.0 (corresponding to σ c = 0.0). RMS is the root mean square difference between magnification maps (compared pixel by pixel), and P(KS ) is the K-S test probability. The CPU-T and CPU-P comparisons are virtually indistinguishable.
CPU-T and GPU-D CPU-P and GPU-D N rays σ RMS P For s ≤ 0.8, we find the same results as the two cases discussed above: for N rays ≥ 500 the results of CPU-T and GPU-D are in excellent agreement. The s = 0.9 case has a slightly higher RMS difference at N rays = 500 of 0.4, although the K-S probability is still > 99.7%. For the s = 0.99 case, however, no agreement is reached between CPU-T and GPU-D. Visual inspection of the maps shows that too few rays have been shot in the CPU-T case. This suggests that more rays should be shot in CPU-T for higher smooth matter fractions; N rays = 1000 should be sufficient for s = 0.90, however more than N rays > 2000 is required in the s = 0.99 case. Once again, CPU-P performs identically to CPU-T.
By comparing the accuracy of the two tree codes, which rely on a series of approximations to the lens distribution in order to speed up calculation time, with GPU-D, which conducts the entire brute-force calculation, we have verified the accuracy of the standard tool for high optical depth microlensing simulations. Indeed, we have considered a wider range of parameters, as represented through the quantity Φ, than previous work (Φ > 2 × 10 14 ). While not unexpected, our results confirm that the tree code can indeed be used with confidence, provided a sufficiently large number of rays per pixel are shot (i.e. typically N rays ≥ 1000 for the parameter space explored here).

Timing Tests
Having established accuracy of the three codes under consideration, we now turn to timing tests. Specifically, we determine the time taken to generate magnification maps with a given set of model parameters, allowing us to determine the regions of (γ, σ, N pix ) parameter space where one of the implementations may be favoured.
We do not consider timing comparisons between CPU-T and CPU-P, as these are reported in Garsden & Lewis (2010). Since CPU-P run-times scale with the number of nodes available, it is difficult to establish which configuration provides the most appropriate comparison between CPU-P and GPU-D. Using the hardware cost as a factor in choosing between software alternatives, we suggest that the cost of purchasing a suitable high-end graphics card plus computer is much lower than the set-up costs for a modest CPU-based computing cluster, once networking overheads and cluster management are considered. We reflect further on the implications this has for the future of microlensing computations in Section 4.
Care must be taken when quoting absolute speeds for different hardware -we use our results as a guide to the relative processing speeds of CPU-T and GPU-D, but note that alternative hardware can change timing by factors of a few. Our benchmarks were conducted on the following equipment: • CPU-T Intel Xeon Dual Core 3060 CPU (2.4 GHz) with 1 GB RAM.
• GPU-D NVIDIA S1070 Tesla unit connected via dual PCIe x8, or higher, buses to an Intel Q6600 Quad Core CPU (2.4 GHz) with GB RAM. The Tesla S1070 comprises four GPUs, each with 240 stream processors, running at 1.296 GHz. Each GPU is able to access 4 GB of RAM. The quoted peak processing rate is 2.488 TFLOP/s, however, for GPU-D, a peak of 1.28 TFLOP/s was obtained using all four GPUs (Thompson et al. 2010). We note that the majority of our timing benchmarks in the present work were conducted using only two of the four GPUs available on the Tesla S1070 4 .
Since we have determined that the magnification maps produced by CPU-T are accurate for N rays ≥ 1000 across the majority of parameter space, we conduct our timing tests by comparing the total time taken by CPU-T and GPU-D to generate identical maps with N rays = 1000, for each combination of parameters. The results of these tests are discussed below.

No shear
Timing results for N pix = 1024 2 , 2048 2 and 4096 2 , with γ = 0 and s = 0, are plotted in Figure 2. Each panel shows both σ (along the top axis), and the logarithm of N * (along the bottom axis). Note that as s = 0 for these simulations, σ = σ * .
The performance of CPU-T is essentially constant as σ increases. The same is not true for GPU-D; it performs better for lower σ (lower N * ) values. This result is as expected; as the number of lenses is increased, the direct code is required to perform correspondingly more calculations, whereas the approximations in the tree code (specifically the depth to which the cell tree is accessed) mean it performs roughly the same number of calculations at each step.
The shape of the GPU-D timing curves is the same as those presented in Figure 2 of Thompson et al. (2010). For low numbers of lenses (N * 10 4 ), overheads induced by the need to transfer data between GPU and CPU affect the performance of GPU-D. For N * 10 4 , near the limit of our comparison, the performance of the code scales linearly with N * . While this linear growth suggests that there is a cross-over in run-time between CPU-T and GPU-D at around N * = 10 5 , it must be recalled that this is approaching the maximum lens limit of CPU-T, due to the memory requirements for storing the cell tree.
We note that although the performance of GPU-D is decreasing as σ (or equivalently, N * ) is increased, it is still faster than CPU-T for all of our γ = 0.0 simulations. Using all four GPUs on the Tesla S1070 unit, instead of the two we use for this detailed timing comparison, results in factor of two speed-up.

Shear and convergence
Real lensing galaxies are unlikely to be perfectly spherically symmetric, and so realistic lens models do not contain zero external shear. We therefore also conducted timing tests for varying γ. Equation 4 shows that the actual calculation of the effect of external shear should not significantly add to computation time -it simply adds one extra operation to the calculation of the deflection of each ray. However, the side lengths L ± of the rectangular shooting grid that is projected on to a square magnification map with side length L map depends on the shear in the following way: When 1 − σ − γ ≈ 0 or 1 − σ + γ ≈ 0, one side of the shooting grid can be very large. This leads to a large area for the lens plane, and hence a correspondingly large number of lenses (see Equation 5 and the discussion following it). We conducted timing tests for the case where N rays = 1000, N pix = 2048 2 , s = 0 and σ = σ * = 0.45, with 0 ≤ γ ≤ 1.0 and ∆γ = 0.1. The results of these tests are presented in Figure  3. The figure shows γ versus the logarithm of the time taken to generate the map. We also plot time against the theoretically predicted average magnification µ av in that figure, where µ av is given by: A negative magnification is interpreted as a parity flip in the image. For low magnifications (when 1 − σ − γ and 1 − σ + γ are very different from zero), GPU-D is faster than CPU-T. There is a crossover point at a magnification of ∼ 5 where CPU-T begins to out-perform GPU-D. For high magnifications, CPU-T is significantly faster than GPU-D. We emphasise again that these timing tests were conducted using only two of the available four GPUs. Running on all four GPUs should improve the performance of GPU-D by approximately a factor of two.
The timing behaviour of GPU-D with the addition of a nonzero shear is exactly the same as the zero shear case. Initially, overheads dominate the code's performance, but as the number of lenses is increased to ∼ 10 4 , the code begins to scale linearly with N * .
The same is not true for CPU-T -the timing behaviour in the zero shear and varying shear cases seems to be quite different. This is due to the way in which the shooting grid is stretched by the addition of a non-zero shear. Equation 11 shows that for certain combinations of σ and γ the shooting grid will be significantly stretched in one direction. A ray shot from one corner of this stretched grid is close to fewer lenses than a ray shot from a corner of a square with the same area and thus a greater number of lenses can be grouped together into pseudolenses. Although N * goes up for these cases, the number of lenses that need to be included directly in the calculation goes down for the tree code.
We expect this behaviour to be generic as σ is varied -CPU-T will be faster than GPU-D in the high magnification regime on current-generation hardware, whereas GPU-D will out-perform CPU-T in the low magnification regime. However, as we discuss in Section 4, the former situation may only be a limit in the short-term.

Smooth matter fraction
Finally, we present the results of the timing tests for varying smooth matter fraction in Figure 4. Tests were conducted for N pix = 2048 2 , σ = 0.45, γ = 0.5, and N rays = 1000. This combination of σ and γ corresponds to the "worst case" magnification for GPU-D timing, and the "best case" magnification for CPU-T timing (see the previous section).
Increasing the smooth matter fraction for fixed σ corresponds to decreasing N * , as mass is transferred from the compact stellar component to the continuously distributed component. We therefore see exactly the timing behaviour we expect: CPU-T takes roughly the same amount of time for each s, whereas GPU-D is faster for fewer lenses.
At the distances from the centre of lensing galaxies where lensed images are typically observed (∼ 2 − 10 kpc), we expect the smooth matter fraction to be quite high. Indeed, a few rough microlensing measurements of smooth matter percentage at lensed image positions exist: 80% to 85% at the location of the A and D images in SDSS J0924+0219 (Keeton et al., 2006); ∼ 80% in HE 1104+080 (Chartas et al., 2009); ∼ 90% in PG 1115+080 (Pooley et al., 2009); and ∼ 70% in RXJ 1131-1231 (Dai et al., 2010). For smooth matter fractions in this range, we find that GPU-D already performs comparably to CPU-T in high magnification systems, and out-performs it in low magnification systems.
We note that our accuracy tests in Section 3.1.3 showed that even N rays = 2000 were not sufficient for CPU-T to obtain results indistinguishable from GPU-D at s = 0.99, at least with the CPU-T accuracy parameter used throughout this paper. The s = 0.99 point in the Figure 4 CPU-T results is therefore likely to be lower than the actual run-time required to produce a useable map.

Discussion
Our intention in this work is not to imply that the tree code is excessively inaccurate. However, its reliance on an approximation to the lens distribution to reduce the number of calculations per magnification map warrants further investigation that has not been feasible until now. The standard tree code (Wambsganss 1990;1999) has been in widespread use throughout the lensing community since the early 1990s. Though it is unsurprising, it is nevertheless reassuring to know that CPU-T remains accurate over an extended range of Φ-space. Our comparison here demonstrates that researchers can, and should, continue to use the tree code for microlensing simulations with confidence. Additionally, we have demonstrated that this accuracy is also maintained by the CPU-P code, which enables exploration of problems orders of magnitude beyond those that can be achieved on a single-core CPU.
However, we stress that the direct inverse ray-shooting code is more accurate than its tree code counterparts by simple virtue that it includes contributions from all lenses without approximations -the entire calculation is done by brute force, taking advantage of the massively parallel GPU architecture to complete the computation in a reasonable timeframe. While CPU-T can indeed be used for a wide range of astrophysically relevant scenarios, GPU-D is now a highly competitive alternative.
Given that we have verified the accuracy of the tree codes, though, it is reasonable to ask in which areas of parameter space each code should be used. The larger the number of lenses, the better the tree codes perform relative to GPU-D. For the values of N pix and N rays we examined in our benchmarking, the crossover point where CPU-T begins to out-perform its GPUbased counterpart is N * ∼ 50, 000 lenses. This can occur either when (σ − γ) or (σ + γ) are close to 1. For scenarios with N pix > 4096 2 and N * > 10 5 , there are advantages in moving to a distributed CPU-based computing cluster, although this can be a financially costly approach. One of the perceived benefits of GPU computing is the heavily reduced dollar per Tflop/s they offer, which makes them an attractive alternative as "desktop" supercomputers (see Thompson et al. 2010). The downside of GPUs, at this point in time, is the paucity of astronomy codes that have been adapted to run on them, whereas CPUbased computing clusters can support a much wider range of existing codes.
Since about 2005, single-core CPU processing speeds have plateaued.
Additional computing performance has been achieved, and indeed Moore's Law growth in processing power has been maintained, by moving to multi-core CPUs. A growing range of applications across a variety of scientific disciplines are now taking advantage of the "many-core" processing architecture of GPUs, which are providing a preview of the future of CPU architectures. Yet if GPUs do indeed point towards a future where many-core processing is mainstream, the need to convert existing trusted, tested, accepted codes will become more pressing. A question that many researchers will have to answer is which code should be adapted?
We suggest that decisions to move other standard astronomy computations to GPU may benefit from a similar approach to Thompson et al. (2010). Rather than attempting to directly port an existing advanced code (i.e. CPU-T) to GPU, they considered a brute force approach that was not feasible for singlecore CPU, but which exhibited the type of fine-grained, massive parallelism that is ideal for GPGPU. This resulted in much faster code development time, despite the learning curve associated with adopting CUDA and the specifics of GPU parallel programming. Moreover, we have shown that in most astrophysically-motivated cases, we are already no worse off using GPU-D compared with CPU-T, with the added benefit of greater intrinsic accuracy.
The N * -crossover in timing tests that we have identified does depend on the hardwared used. For example, our benchmarking was performed using two of the four GPUs available on the Tesla S1070, while the CPUs used for CPU-T were a few years old. The key point, however, is that the GPU-D is already parallel -GPU performance will increase, and more GPU cores will be placed in individual machines, yet no further code development is required. The single-core tree code, however, is more heavily restricted -CPU-T can no longer make significant processing gains by simply relying on Moore's Law to make single-core CPU hardware faster.

Conclusion
In this work, we have examined the accuracy of three highly competive alternatives for generating microlensing magnification maps: a widely-used hierachical tree code for single-core CPUs; a parallel, large data implementation for a distributed CPU-based computing cluster; and a brute force solution that utilises the massive parallelism of modern graphics processing units. We have demonstrated that these three codes do indeed present comparable accuracy, and choice of approach will depend on considerations relating to the scale and nature of the problem to be investigated.
While considered impractical as a calculating approach even a few years ago, our timing tests show that direct, inverse rayshooting is now a viable, and indeed more accurate, option. Moving a simple algorithm to GPU, rather than attempting to directly implement a more advanced solution -originally developed to overcome a hardware restriction of single-core CPUsto achieve run-times that are already similar, suggests an approach that other early adopters of GPUs for astronomy computation may wish to consider.