M Ourak

B Tamadazte

N Andreff

C Pieralli

O Gaiffe

B Wacogne

Towards registration of multimodal images of vocal folds based on mutual information

This paper deals with mutual information-based registration of multimodal images for laser phonomicrosurgery of the vocal folds. The images to be registered are white light images (white light camera) versus fluorescence images. This work is carried out within the framework of the European project µRALP which involves the use of microrobotic system for endoluminal laser phonosurgery. The designed system includes two fiber bundles connected to a high speed camera and one fiber bundle used for fluorescence image. Using the mutual information based registration method, it will be possible to represent the visible information in the fluorescence image and use it in the other image.

I. INTRODUCTION

A phonosurgery system uses a stereo-microscope and a semi-automatic system to control the laser source positioned at 400 mm outside the patient (see Fig. 1). This technique relies completely on the skills of surgeons who must operate through a microscope, control the laser directly by hand, and deals with the associated poor ergonomics of the operating setup [START_REF] Burns | Endoscopic Laser Resection of Laryngeal Cancer: Is It Oncologically Safe? Position Statement From the American Broncho-Esophagological Association[END_REF], [START_REF] Jackel | Twenty-five Years Experience with Laser Surgery for Head and Neck Tumors[END_REF]. The µRALP project will enable to overcome such limitations by developing a more intuitive system (inserted through the mouth into the laryngo-pharyngeal cavity) that allows to position the laser close to the vocal chords (ie. 20 mm). Mainly, the proposed system includes a piezoelectric actuated mirror (two degrees of freedom), two fiber bundles for high frequency imaging, one fiber bundle for fluorescence image, and a Helium-Neon (HeNe) laser (see Fig. 2).

The developed system will work as follows: the surgeon defines a trajectory around the tumor to incise, in the fluorescence image. Thereby, the objective of µRALP system is to control automatically the displacement of the laser beam along the pre-defined trajectory. In order to ensure a better accuracy during automatic resection, it is required to close the control loop over an exteroceptive sensor (high speed camera in our case). Namely, this so-called visual servoing (vision-based robot control) is robust to robot calibration errors and to robot environment variations [START_REF] Chaumette | Visual Servo Control, Part I: Basic Approaches[END_REF].

In general, visual servoing requires defining a reference features (points, lines, image, 3D pose, etc.) that allows the system to evolve from its current position to its target position. In our case, the reference can be defined (by the surgeon) in the fluorescence image frame R f and the system will move in the high speed white light image frame R s . Thus, it is necessary to compute the transformation ( T , ρ) between both frames (ie. R f and R s ). The mutual informationbased method seems very appropriate for the computation of this transformation, especially when the images are of different modalities.

In the literature, the mutual information-based registration has been widely discussed. Zitova et al. [START_REF] Zitova | Image Registration Methods: a Survey[END_REF] has classified the registration techniques for medical applications into two main categories: area-based and feature-based methods. Generally, a registration technique follows mainly these four steps: features detection, features matching, transformation estimation, and image resampling.

This paper focuses in the adaptation of the mutual information technique in the case of fluorescence images vs. high speed white light images. Thereby, it contains some basic definitions about the mutual information, the context of use, some primary results and the future investigations of this work.

II. MUTUAL INFORMATION

The mutual information is inspired by the information theory domain. It is based on the measure of information, commonly called entropy, in a message. In the case of image processing, this allows to measure the similarity between two images I 1 and I 2 .

A. Entropy

In information theory domain, the entropy was first introduced by Hartley [START_REF] Hartley | Transmission of Information[END_REF]. There, he used it to measure the uncertainty of a signal s between a transmitter and a receiver. The first Hartley entropy function, is given by:

H 0 (X) = log b |X| (1)
where X is a vector containing random values x i .

Move over, the mean information quantity associated to each source value is the expectation of the own information of each event X = x i . According to Shannon, adding a weight to the output of the Hartley entropy allows to rewrite the equation [START_REF] Burns | Endoscopic Laser Resection of Laryngeal Cancer: Is It Oncologically Safe? Position Statement From the American Broncho-Esophagological Association[END_REF] as:

H = i p i log 2 ( 1 p i ) = - i p i log 2 (p i ) (2)
where p i is the associated probability of each event x i . Applied to the images, the measure of the entropy allows us to compute the degree of similarity information between two images (higher is the similarity between images, lower is the entropy).

B. Joint histogram

To estimate the joint probability distribution of the gray values between two images, it is possible to use the notion of joint histogram. This joint histogram is a 3D map (N 2 in N), whose 2D inputs are the grayscale value of each pixel of I 1 (i, j) and of the corresponding pixel in I 2 (i, j) and the output is the amount of such couple (I 1 (i, j), I 2 (i, j)) (see Fig. 3 and algorithm 1).

Thereby, the joint histogram allows illustrating the similarity between two images. In practice, lower the similarity between images, more the joint histogram is uniform (its for j from 0 to v do 7:

h(I 1 (i, j), I 2 (i, j)) = h(I 1 (i, j), I 2 (i, j)) + 1;

8:

end for 9: end for 10: return h graph is a plane). Conversely, if two images are perfectly similar, the joint histogram is represented by a diagonal line.

C. Mutual information

Based on the entropy and the joint histogram principles, the mutual information can be defined following the three methods described in [START_REF] Josien | Mutual Information based Registration of Medical Images: a Survey[END_REF], [START_REF] Dame | Second Order Optimization of Mutual Information for Real-time image Registration[END_REF]. However, for our work, we use the method using the marginal and joint entropy based on the value of the normalized joint histogram, given by (3).

M I(I 1 , I 2 ) = H(I 1 ) + H(I 2 ) -H(I 1 , I 2 ) (3)
where H(I 1 ) and H(I 2 ) are the marginal entropy of I 1 and I 2 , respectively, and H(I 1 , I 2 ) is joint entropy computed from the joint histogram. Let us assume the rigid transformation T SE(3) between the images I 1 and I 2 . Thereby, the transformation can be estimated by maximizing the mutual information (M I) i .

T = arg max t∈SE(3) M I[I 1 , t(I 2 )] ( 4
)
where t is a possible rigid transformation. If we add a zoom factor, then we need to multiply the transformation matrix with the scaling factor:

zoom =   ρ 0 0 0 ρ 0 0 0 1   ( 5
)
where ρ is the scale value. Hence, the equation ( 4) can evolve using the zoom parameter as:

( T , ρ) = arg max t,ρ M I[I 1 , ρ(t(I 2 ))] (6)
III. OPTIMIZATION Once the basics of the mutual information are addressed, it is necessary to establish an effective optimization method able to compute the maximum (respectively the minimum) of the joint information between both images. In this work, we opted to use a non-gradient-based optimization technique commonly known as the simplex method [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF]. This choice is justified by the fact that the simplex is better fitted to an optimization problem with limited variables number [START_REF] Zhu | Influence of Implementation Parameters on Registration of MR and SPECT Brain Images by Maximization of Mutual Information[END_REF]. The different steps of the presented mutual information based images registration are shown in Fig. 4.

IV. FIRST RESULTS

The experimental validation is performed with the objective to determine the estimated transformation ( T , ρ) applied to image I 2 to be confounded with the image I 1 . This transformation is here decomposed in translations (t x and t y ), rotation θ z and zoom ρ.

The validation tests were carried out using different images, between monomodal images acquired in various conditions (point of view, lighting illumination, adding objects onto the scene and scale changing), and between fluorescence images vs. white light image of vocal chords. This paper show only the results using fluorescence vs. white light images. In this experiment, we added a fictive incision mark on the fluorescence image where the tumor is visible. To display the difference between image I 1 and the transformed image ρ(t(I 2 )), we combine them as their checkerboard masked sum (Fig. 5). Thus a quality criterion of the registration is the continuity of the combined image a long the check-board lines. By analyzing the Fig. 6, it can be noticed the continuity of the combination of I 1 in I 2 as shown in Fig. 6(b).

Otherwise, the obtained values of the different parameters of ( T , ρ) for the test shown in Fig. 6 are: t x = -13.92 pixels, t y = -1.31 pixels, θ z = -0.0141 • , and ρ = 0.158 pixels. This result is obtained in 110 seconds for a 170×150 image size using a 2.5 GHz PC.

V. CONCLUSION

This paper presented the initial results using the proposed mutual information registration method. The considered images are multimodal (ie. fluorescence image vs. white light image), and such image pairs will be used to detect the cancerous tumor in the vocal chords. The mutual information based registration method gives promising results and shows interesting robustness with respect to different kind and quality of the images.

The next step is to develop a new metric for visual servoing using the proposed technique. A new control law, which does not require any matching nor tracking step (direct visual servoing), based on mutual information will be designed. This is to ensure that the surgeon can plan a trajectory (ie. incision) on the fluorescence image and the endoluminal laser phonosurgery microrobot system will be able to follow this reference in the high speed white light images.

Fig. 1 .

 1 Fig. 1. Current laser phonosurgery setup.

Fig. 2 .

 2 Fig. 2. New laser phonosurgery system proposed by the µRALP project.

Fig. 3 .

 3 Fig. 3. Representation of the joint histogram.

Fig. 4 .

 4 Fig. 4. Key steps of the proposed mutual information based registration.

2 Fig. 5 .

 25 Fig. 5. Illustration of the construction of the result image (⊕ operator to combine the two images I 1 and I 2 ).

Fig. 6 .

 6 Fig. 6. Illustration of the registration results between I 1 and I 2 . Column (a) represents the test images (ie. fluorescence image vs. white light image), and column (b) represents the obtained combined image with a zoom in to show the details.

  Algorithm 1 Joint histogram algorithm Require: I 1 , I 2 two images 1: init : binSize ← 255;

2: init : h(binSize, binSize); 3: u ← height(I 1 ); 4: v ← width(I 1 ); 5: for i from 0 to u do 6:

ACKNOWLEDGMENT This work was supported by µRALP, the EC FP7 ICT Collaborative Project no. 288663 (http://www.microralp.eu), and by ACTION, the French ANR Labex no. ANR-11-LABX-01-01 (http://www.labex-action.fr).


Simplice Donfack
email: sdonfack@eecs.utk.edu

Jack Dongarra
email: dongarra@eecs.utk.edu

Mathieu Faverge
email: mathieu.faverge@inria.fr

Mark Gates
email: mgates3@eecs.utk.edu

Jakub Kurzak
email: kurzak@eecs.utk.edu

Piotr Luszczek
email: luszczek@eecs.utk.edu

Ichitaro Yamazaki
email: iyamazak@eecs.utk.edu

On Algorithmic Variants of Parallel Gaussian Elimination: Comparison of Implementations in Terms of Performance and Numerical Properties *

Keywords: D.1.3 [Concurrent Programming]: Parallel programming, F.2.1 [Numerical Algorithms and Problems]: Computations on matrices, G.4 [MATHEMATICAL SOFTWARE]: Parallel and vector implementations Algorithms, Implementation Gaussian elimination, LU factorization, parallel, shared memory, multicore

Gaussian elimination is a canonical linear algebra procedure for solving linear systems of equations. In the last few years, the algorithm received a lot of attention in an attempt to improve its parallel performance. This article surveys recent developments in parallel implementations of the Gaussian elimination. Five different flavors are investigated. Three of them are based on different strategies for pivoting: partial pivoting, incremental pivoting, and tournament pivoting. The fourth one replaces pivoting with the Random Butterfly Transformation, and finally, an implementation without pivoting is used as a performance baseline. The technique of iterative refinement is applied to recover numerical accuracy when necessary. All parallel implementations are produced using dynamic, superscalar, runtime scheduling and tile matrix layout. Results on two multisocket multicore systems are presented. Performance and numerical accuracy is analyzed.

INTRODUCTION

Gaussian elimination has a long history that can be traced some 2000 years back [START_REF] Grcar | Mathematicians of Gaussian elimination[END_REF]. Ever since, dense systems of linear equations have been a critical cornerstone for some of the most compute intensive applications. Any improvement in the time to solution for these dense linear systems has a direct impact on the execution time of numerous applications. A short list of domains directly using dense linear equations to solve some of the most challenging problems our society faces are: airplane wing design, radar cross-section studies, flow around ships and other off-shore constructions, diffusion of solid bodies in a liquid, noise reduction, and diffusion of light by small particles.

Computational experiments of self-sustaining fusion reactions could give us an informed perspective on how to build a device capable of producing and controlling the high performance [START_REF] Aymar | Overview of ITER-FEAT -the future international burning plasma experiment[END_REF]. Modeling the heating response of plasma due to radio frequency waves in the fast wave time scale leads to solving the generalized Helmholtz equation. The time harmonic terms of effective approximations of the electric field, magnetic field, and distribution function as a time-averaged equilibrium satisfy the equation. The Scientific Discovery through Advanced Computing project (SciDAC) Numerical Computation of Wave Plasma-Interactions in Multi-dimensional Systems developed and implemented a simulation code that gives insight into how electromagnetic waves can be used for driving current flow, heating and controlling instabilities in the plasma. The code is called AORSA [START_REF] Jaeger | Advances in full-wave modeling of radio frequency heated multidimensional plasmas[END_REF][START_REF] Jaeger | Sheared poloidal flow driven by mode conversion in Tokamak plasmas[END_REF][START_REF] Jaeger | Global-wave solutions with self-consistent velocity distributions in ion cyclotron heated plasmas[END_REF] and stands for All ORders Spectral Algorithm. The resulting computation requires a solution of a sys-tem of linear equations exceeding half a million unknowns and the fastest method is LU factorization through Gaussian elimination with partial pivoting [START_REF] Barrett | Complex version of high performance computing LINPACK benchmark (HPL)[END_REF].

The electromagnetic community is a major user of dense linear systems solvers. Of particular interest to this community is the solution of the so-called radar cross-section problem -a signal of fixed frequency bounces off an object; the goal is to determine the intensity of the reflected signal in all possible directions. The underlying differential equation may vary, depending on the specific problem. In the design of stealth aircraft, the principal equation is the Helmholtz equation. To solve this equation, researchers use the method of moments [START_REF] Harrington | Origin and development of the method of moments for field computation[END_REF][START_REF] Wang | Generalized Moment Methods in Electromagnetics[END_REF]. In the case of fluid flow, the problem often involves solving the Laplace or Poisson equation. Here, the boundary integral solution is known as the panel methods [START_REF] Hess | Panel methods in computational fluid dynamics[END_REF][START_REF] Hess | Calculation of potential flows about arbitrary bodies[END_REF], so named from the quadrilaterals that discretize and approximate a structure such as an airplane. Generally, these methods are called boundary element methods. Use of these methods produces a dense linear system of size O(N) by O(N), where N is the number of boundary points (or panels) being used. It is not unusual to see size 3N by 3N, because of three physical quantities of interest at every boundary element. A typical approach to solving such systems is to use LU factorization. Each entry of the matrix is computed as an interaction of two boundary elements. Often, many integrals must be computed. In many instances, the time required to compute the matrix is considerably larger than the time for solution. The builders of stealth technology who are interested in radar cross-sections are using direct Gaussian elimination methods for solving dense linear systems. These systems are always symmetric and complex, but not Hermitian. Another major source of large dense linear systems is problems involving the solution of boundary integral equations [START_REF] Edelman | Large dense numerical linear algebra in 1993: the parallel computing influence[END_REF]. These are integral equations defined on the boundary of a region of interest. All examples of practical interest compute some intermediate quantity on a two-dimensional boundary and then use this information to compute the final desired quantity in three-dimensional space. The price one pays for replacing three dimensions with two is that what started as a sparse problem in O(n 3 ) variables is replaced by a dense problem in O(n 2 ). This article is organized as follows: at the beginning the motivation for this survey is given in Section 2; then the related work is discussed in Section 3; the original contribution of this work is described in Section 4; Section 5 presents the algorithms evaluated here and Section 6 gives implementation details; Section 7 contains the performance and numerical results of the experiments. Finally, Section 8 concludes the article and Section 9 shows potential future extensions of this work.

MOTIVATION

The aim of this article is to give exhaustive treatment to both performance and numerical stability of various pivoting strategies that have emerged over the past few years to cope with the need for increased parallelism in the face of the paradigm shifting switch to multicore hardware [START_REF] Buttari | The Impact of Multicore on Math Software[END_REF]. Rather than focusing on multiple factorizations and performance across multiple hardware architectures [START_REF] Agullo | Comparative study of one-sided factorizations with multiple software packages on multi-core hardware[END_REF], we switch our focus to multiple pivoting strategies and provide uniform treatment for each, while maintaining sufficient hardware variability to increase the meaningful impact of our results Ever since the probabilistic properties of partial pivoting [START_REF] Trefethen | Average case analysis of Gaussian elimination[END_REF] were established, the community has fully embraced the method despite the high upper bound on pivot growth that it theoretically can incur [START_REF] Foster | Gaussian elimination with partial pivoting can fail in practice[END_REF]: O(2 n ) versus a much more acceptable O(n 3 ) provided by complete pivoting [START_REF] Golub | Matrix Computations[END_REF]. Probabilistically, similar results have been proven for no-pivoting LU [START_REF] Yeung | Probabilistic analysis of Gaussian elimination without pivoting[END_REF], which lead to include it in this survey, but complete lack of any pivoting strategy is discouraged for practical applications.

RELATED WORK

Incremental pivoting [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] has its origins in pair-wise pivoting [START_REF] Sorensen | Analysis of pairwise pivoting in Gaussian elimination[END_REF], both of which have numerical issues relating to the pivot growth factor [START_REF] Agullo | LU factorization for acceleratorbased systems[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]. But they offer a reduction of dependencies between tasks that aids parallelism that has become so important with the introduction and current proliferation of multicore hardware.

The partial pivoting code that we use for evaluation in this article is based on the parallel panel factorization that uses a recursive formulation of LU [START_REF] Dongarra | Recursive approach in sparse matrix LU factorization[END_REF][START_REF] Gustavson | Recursion leads to automatic variable blocking for dense linear-algebra algorithms[END_REF]. This recent implementation was introduced to address the bottleneck of the panel computation [START_REF] Dongarra | Exploiting fine-grain parallelism in recursive LU factorization[END_REF][START_REF] Dongarra | Exploiting fine-grain parallelism in recursive LU factorization[END_REF], and has been successfully applied to matrix inversion [START_REF] Dongarra | High performance matrix inversion based on lu factorization for multicore architectures[END_REF]. However, it should not be confused with a globally recursive implementation based on column-cyclic distribution [START_REF] Luszczek | Anatomy of a globally recursive embedded LINPACK benchmark[END_REF]. Neither it is similar to a non-recursive parallelization effort [START_REF] Castaldo | Scaling LAPACK panel operations using parallel cache assignment[END_REF] that focuses only on cacheefficient implementation of the existing Level 1 BLAS kernels rather than using the recursive formulation that was used in this article's experiments.

With incremental pivoting, the numerical quality of the solution does not match that of the standard partial pivoting scheme in the panel factorization because it has to be replaced by a form of pair-wise pivoting [START_REF] Sorensen | Analysis of pairwise pivoting in Gaussian elimination[END_REF], which is related to an updating-LU for out-of-core computations [START_REF] Yip | FORTRAN Subroutines for Out-of-Core Solutions of Large Complex Linear Systems[END_REF] when blocking of computations has to be done for better performance. It has resurfaced in the form of what is now called the pivoting incremental strategy [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] that allows pairwise treatment of tiles and, as a consequence, reduces dependencies between tasks and aids parallelism. This causes the magnitude of pivot elements to increase substantially, which is called the pivot growth factor, and results in rendering the factorization numerically unstable [START_REF] Agullo | LU factorization for acceleratorbased systems[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].

A probabilistic technique called the Partial Random Butterfly Transformation (PRBT) is an alternative to pivoting and may, in a sense, be considered a preconditioning step that renders pivoting unnecessary. It was originally proposed by Parker [START_REF] Parker | A randomizing butterfly transformation useful in block matrix computations[END_REF] and then applied in practice with adaptations by Baboulin et al [START_REF] Baboulin | Accelerating linear system solutions using randomization techniques[END_REF], that limited the recursion depth without compromising the numerical properties of the method.

The tournament pivoting originated in CALU [START_REF] Demmel | Implementing communication-optimal parallel and sequential QR factorizations[END_REF][START_REF] Donfack | Adapting communication-avoiding LU and QR factorizations to multicore architectures[END_REF][START_REF] Grigori | Communication avoiding Gaussian elimination[END_REF] -Communication-Avoiding LU. The main design goal for this new pivoting scheme was to attain the minimum bounds on the amount of data communicated and the number of messages exchanged between the computing processors. The way to achieve these bounds was to minimize the communication that occurs during the factorization of the panel by performing redundant computations. The extra operations do not hurt the scheme in the global sense because they only account for an increase in the constant for a lower term: O(n 2 ), while the highest order term -O(n 3 ) -remains unchanged. In terms of stability, CALU could potentially cause a greater pivot growth, which increases exponentially with the amount of parallelism. Unlike the partial pivoting, however, the bounds might not be attainable in practice.

Similar study was published before [START_REF] Agullo | Comparative study of one-sided factorizations with multiple software packages on multi-core hardware[END_REF] but it was mostly focused on performance across a wide range of hardware architectures. It featured results for the main three of the one-sided factorization schemes: Cholesky, QR, and LU. No algorithmic variants for a particular method were considered.

ORIGINAL CONTRIBUTION

The unique contribution of this survey is in implementing all the algorithms, being compared using the same framework, the same data layout, and the same set of parallel layout translation routines, as well as the same runtime scheduling system. This allows for gaining a level of insight into the trade-offs of the different methods that one could not reach by comparing published data for different implementations in different environments.

ALGORITHMS 5.1 Partial Pivoting

The LAPACK block LU factorization is the main point of reference here, and LAPACK naming convention is followed. The LU factorization of a matrix A has the form

PA = LU,
where L is a unit lower triangular matrix, U is an upper triangular matrix and P is a permutation matrix. The LAPACK algorithm proceeds in the following steps: Initially, a set of NB columns (the panel) is factored and a pivoting pattern is produced (implemented by the DGETF2 routine). Then the elementary transformations, resulting from the panel factorization, are applied in a block fashion to the remaining part of the matrix (the trailing submatrix). This involves swapping of up to NB rows of the trailing submatrix (DLASWP), according to the pivoting pattern, application of a triangular solve with multiple right-hand-sides to the top NB rows of the trailing submatrix (DTRSM), and finally application of matrix multiplication of the form A i j ← A i j -A ik × A k j (DGEMM), where A ik is the panel without the top NB rows, A k j is the top NB rows of the trailing submatrix and A i j is the trailing submatrix without the top NB rows. Then the procedure is applied repeatedly, descending down the diagonal of the matrix (Figure 1). The block algorithm is described in detail in section 2.6.3 of the book by Demmel [START_REF] Demmel | Applied Numerical Linear Algebra[END_REF] U (done)

L (done) A ij A kj A ik

Incremental Pivoting

The worst performance-limiting aspect of Gaussian elimination with partial pivoting is the panel factorization operation. First, it is an inefficient operation, usually based on a sequence of calls to Level 2 BLAS. Second, it introduces synchronization, by locking the entire panel of the matrix at a time. Therefore, it is desirable to split the panel factorization into a number of smaller, finer-granularity operations, which is the basic premise of the incremental pivoting implementation, also known in literature as the tile LU factorization.

In this algorithm, instead of factoring the panel one column at a time, the panel is factored one tile at a time. The operation proceeds as follows: First the diagonal tile is factored, using the standard LU factorization procedure. Then the factored tile is combined with the tile directly below it, and factored. Then the re-factored diagonal tile is combined with the next tile, and factored again. The algorithm descends down the panel until the bottom of the matrix is reached. At each step, the standard partial pivoting procedure is applied to the tiles being factored. Also, at each step, all the tiles to the right of the panel are updated with the elementary transformations resulting from the panel operations. This way of pivoting is basically the idea of pairwise pivoting applied at the level of tiles, rather than individual elements (Figure 2). The main benefit comes from the fact that updates of the trailing submatrix can proceed alongside panel factorizations, leading to a very efficient parallel execution, where multiple steps of the algorithm are smoothly pipelined.

Tournament Pivoting

The panel factorization is one of the most important tasks, because it creates parallelism for the update of the trailing submatrices. Hence, its ineffective execution suffices to reduce considerably the performance of the overall algorithm. Classic approaches that implement partial pivoting algorithm spend more time to perform communication during the panel factorization and hence are not optimal. This is because pivoting forces the algorithm to factor the panel column by column, and then this leads to an algorithm which communicates asymptotically more than the established lower bounds [START_REF] Demmel | Implementing communication-optimal parallel and sequential QR factorizations[END_REF].

The basic idea of communication avoiding algorithms, initially introduced for distributed memories [START_REF] Demmel | Implementing communication-optimal parallel and sequential QR factorizations[END_REF][START_REF] Grigori | Communication avoiding Gaussian elimination[END_REF], and later adapted to shared memories [START_REF] Donfack | Adapting communication-avoiding LU and QR factorizations to multicore architectures[END_REF], is to replace the search for maximum, performed at each column, by a single reduction of the maximums altogether. This is done thanks to a new pivoting strategy referred to as tournament pivoting (TSLU), which performs redundant computations and is shown to be stable in practice. TSLU reduces the bottleneck introduced by the pivoting operation through a block reduction operation to factor the panel. It factors the panel in two steps. The first one identifies rows, which can be used as good pivots for the factorization of the whole panel, with a tournament selection. The second one swaps the selected pivot to the top of the panel, and then factors the entire panel without pivoting in a tiled Cholesky-like operation. With this strategy, the panel is efficiently parallelized and the communication is provably minimized.

Figure 3 presents the first step of TSLU for a panel W using a binary tree for the reduction operation. First, the panel is partitioned into T r = 4 blocks, that is, W = [W 00 ,W 10 , ...,W T r -1,0 ], where T r represents the number of threads participating in the operation. At the first step of the reduction operation, each thread I applies Gaussian elimination with partial pivoting to its block W I0 , then the resulting permutation matrix P I0 is applied to the original block W I0 , and the first b rows of the permuted block P I0 W I0 are selected as pivot candidates. These first pivot candidates represent the leaves of the reduction tree. In the next steps, W IJ represents the block owned by thread I at step J of the reduction operation. At each node of the tree, the pivot candidates of its descending children are merged one on top of another, and then Gaussian elimination with partial pivoting is applied on the merged block, the resulting permutation matrix is again applied on the original merged block and then, the first b rows are selected as a new pivot candidate. By using a binary tree, this step is repeated log T r times. The pivots obtained at the root of the tree are then considered as the good pivots for the whole panel. Once these pivots are permuted at the top of the panel, each thread I applies Gaussian elimination without partial pivoting to its block W I0 .

The example presented in Figure 3 uses a binary tree with two tiles reduced together at each level, but any reduction tree can be used depending on the underlying architecture. The TSLU implementation in PLASMA, used for experiments in this paper, reduces tiles four by four at each level. The number of 4 tiles has been chosen because it gave a good ratio of kernel efficiency over one single core relative to the time spent to perform the factorization of the subset.

Random Butterfly Transform

As an alternative to pivoting, the Partial Random Butterfly Transformation (PRBT) preconditions the matrix as A r = W ⊤ AV , such that, with probability close to 1, pivoting is unnecessary. This technique was proposed by Parker [START_REF] Parker | A randomizing butterfly transformation useful in block matrix computations[END_REF] and later adapted by Baboulin et al [START_REF] Baboulin | Accelerating linear system solutions using randomization techniques[END_REF].

An n × n butterfly matrix is defined as

B (n) = 1 √ 2 R S R -S ,
where R and S are random diagonal, nonsingular matrices. W and V are recursive butterfly matrices of depth d, defined by

W (n,d) =     B (n/2 d-1 ) 1 . . . 0 . . . . . . . . . 0 . . . B (n/2 d-1 ) 2 d-1     × • • • × B (n) .
We use a depth d = 2, previously found to be sufficient in most cases [START_REF] Baboulin | Accelerating linear system solutions using randomization techniques[END_REF]. Since each R and S is diagonal, W and V can be stored as n × d arrays. Due to the regular sparsity pattern, multiplying an m × n matrix by W and V is an efficient, O(mn) operation.

After applying the PRBT, Gaussian elimination without pivoting is used to obtain the solution, as indicated in Algorithm 1.

Algorithm 1 Solving Ax = b using RBT.

1: A r = W ⊤ AV 2: factor A r = LU without pivoting 3: solve LUy = W ⊤ b for y 4: x = V y
While the randomization reduces the need for pivoting, the lack of pivoting can still be unstable, so we use iterative refinement to reduce the potential for instability. The cost of pivoting is thus avoided, at the expense of applying the PRBT and iterative refinement.

No Pivoting

This implementation of Gaussian elimination completely abandons pivoting. This can be done very rarely in practice without risking serious numerical consequences, or even a complete break-down of the algorithm if a zero is encountered on the diagonal. Here the implementation serves only as a performance baseline. Dropping pivoting increases performance for two reasons. First, the overhead of swapping matrix rows disappears. Second, the level of parallelism dramatically increases, since the panel operations now become parallel, and can also be pipelined with the updates to the trailing submatrix.

Iterative Refinement

Iterative refinement is an iterative method proposed by James Wilkinson to improve the accuracy of numerical solutions to systems of linear equations. When solving a linear system Ax = b, due to the presence of roundoff errors, the computed solution may deviate from the exact solution. Starting with the initial solution, iterative refinement computes a sequence of solutions that converges to the exact solution when certain assumptions are met (Algorithm 2).

Algorithm 2 Iterative refinement using MATLAB TM backslash notation. 1:

repeat 2: r = b -Ax 3: z = L\(U\Pr) 4:
x = x + z 5: until x is "accurate enough" As Demmel points out [13, pp.60], the iterative refinement process is equivalent to Newton's method applied to f (x) = b -Ax. If the computation was done exactly, the exact solution would be produced in one step. Iterative refinement was studied by Wilkinson [START_REF] Wilkinson | The Algebraic Eigenvalue Problem[END_REF], Moler [START_REF] Moler | Iterative refinement in floating point[END_REF], and Demmel [START_REF] Demmel | Error bounds from extra-precise iterative refinement[END_REF], and is covered in the books by Higham [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] and Stewart [START_REF] Stewart | Introduction to Matrix Computations[END_REF].

Iterative refinement introduces a memory overhead. Normally, in the process of factorization, the original matrix A is overwritten with the L and U factors. However, the original matrix is required, in the refinement process, to compute the residual. Therefore, application of iterative refinement doubles the memory requirement of the algorithm.

IMPLEMENTATION 6.1 Tile Layout

It is always beneficial for performance to couple the algorithm with a data layout that matches the processing pattern. For tile algorithms, the corresponding layout is the tile layout, developed by Gustavson et al. [START_REF] Gustavson | Parallel and cache-efficient in-place matrix storage format conversion[END_REF] and shown in Figure 4. The matrix is arranged in square submatrices, called tiles, where each tile occupies a continuous region of memory. The particular type of layout used here is referred to as Column-Column Rectangular Block (CCRB). In this flavor of the tile layout, tiles follow the column-major order and elements within tiles follow the column-major order. The same applies to the blocks A 11 , A 21 , A 12 , and A 22 . Form the standpoint of serial execution, tile layout minimizes conflict cache misses, because two different memory locations within the same tile cannot be mapped to the same set of a set-associative cache. The same applies to the Translation Look-aside Buffer (TLB) misses. In the context of parallel execution, tile layout minimizes the probability of false sharing, which is only possible at the beginning and end of the continuous memory region occupied by each tile, and can easily be eliminated altogether, if the matrix is aligned to cache lines and tiles are divisible by the cache line size. Tile layout is also beneficial for prefetching, which in the case of strided memory access is likely to generate useless memory traffic.

It is only fair to assume that most users of shared memory systems assemble their matrices in the standard column-major layout, common to FORTRAN 77 and LAPACK, Therefore, the overhead of translating the matrix from the column-major layout to the CCRB layout and back is always included in the timing. Because the entire matrix occupies a contiguous region of memory, translation between the tile layout and the legacy FORTRAN 77 layout can be done in place, without changing the memory footprint. Gustavson et al. [START_REF] Gustavson | Parallel and cache-efficient in-place matrix storage format conversion[END_REF] devised a collection of routines for performing this translation in a parallel and cache efficient manner. It is important to observe that the layout translation routines have a broader impact in forming the basis for a fast transposition operation. The codes are distributed as part of the PLASMA library.

Dynamic Scheduling

In order to exploit the fine-grained parallelism to its fullest, efficient multithreading mechanisms have to be designed, where data dependencies are preserved, i.e., data hazards are prevented. This has been done for both the simpler single-sided factorizations, such as Cholesky, LU and QR [2, 3, 8, 9, 16-18, 29, 37], as well as the more complicated two-sided factorizations, such as the reductions to band bi-diagonal and band tri-diagonal form [START_REF] Haidar | Parallel reduction to condensed forms for symmetric eigenvalue problems using aggregated fine-grained and memory-aware kernels[END_REF][START_REF] Haidar | A comprehensive study of task coalescing for selecting parallelism granularity in a two-stage bidiagonal reduction[END_REF][START_REF] Ltaief | Parallel band two-sided matrix bidiagonalization for multicore architectures[END_REF][START_REF] Ltaief | High Performance Bidiagonal Reduction using Tile Algorithms on Homogeneous Multicore Architectures[END_REF][START_REF] Ltaief | Enhancing parallelism of tile bidiagonal transformation on multicore architectures using tree reduction[END_REF][START_REF] Luszczek | Two-stage tridiagonal reduction for dense symmetric matrices using tile algorithms on multicore architectures[END_REF]. The process of constructing such schedules through manipulation of loop indexes and enforcing them by progress tables is tedious and error-prone. Using a runtime dataflow scheduler is a good alternative. A superscalar scheduler is used here. Superscalar schedulers exploit multithreaded parallelism in a similar way as superscalar processors exploit Instruction Level Parallelism (ILP). Scheduling proceeds under the constraints of data hazards: Read after Write (RaW), Write after Read (RaW) and Write after Write (RaW). In the context of multithreading, superscalar scheduling is a way of automatically parallelizing serial code. The programmer is responsible for encapsulating the work in side-effectfree functions (parallel tasks) and providing directionality of their parameters (input, output, input-and-output), and the scheduling is left to the runtime. Scheduling is done by conceptually exploring the Directed Acyclic Graph (DAG), or task graph, of the problem. In practice the DAG is never built entirely, and instead explored in a sliding window fashion. The superscalar scheduler used here is the QUeuing And Runtime for Kernels (QUARK) [START_REF] Yarkhan | QUARK users' guide: QUeueing And Runtime for Kernels[END_REF] system, developed at the University of Tennessee.

EXPERIMENTAL RESULTS

Hardware and Software

The experiments were run on an Intel system with 16 cores and an AMD system with 48 cores. The Intel system has two sockets of eight-core Intel Sandy Bridge CPUs clocked at 2.6 GHz, with a theoretical peak of 16 cores × 2.6 GHz × 8 ops per cycle ≃ 333 G f lop/s in double precision arithmetic. The AMD system has eight sockets of six-core AMD Istanbul CPUs clocked at 2.8 GHz, with a theoretical peak of 48 cores × 2.8 GHz × 4 ops per cycle ≃ 538 G f lop/s in double precision arithmetic.

All presented LU codes were built using the PLASMA framework, relying on the CCRB tile layout and QUARK dynamic scheduler. The GCC compiler was used for compiling the software stack and Intel MKL (Composer XE 2013) was used to provide an optimized implementation of serial BLAS.

Performance

We now study the performance of our implementations on square random matrices in double real precision, and compare their performance with that of MKL. Since each of our implementations uses a different pivoting strategy, for a fair performance comparison, we ran iterative refinement with all the algorithms. Namely, for MKL, we used the MKL iterative refinement routine DGERFS, while for the tile algorithms, we implemented a tiled iterative refinement with the same stopping criteria as that in DGERFS, i.e., the iteration terminates when one of the following three criteria is satisfied: (1) the component-wise backward error, max i |r i |/(|A|| x| + |b|) i , is less than or equal to ((n + 1) * sfmin)/eps, where r is the residual vector (i.e, r = A xb with the computed solution x), eps is the relative machine precision, and sfmin is the smallest value such that 1/sfmin does not overflow, (2) the component-wise backward error is not reduced by half, or (3) the number of iterations is equal to ten. For all of our experiments, the iterative refinements converged in less than ten iterations. We observed that even with partial pivoting, it requires a couple of iterations to satisfy this stopping criteria (see Section 7.3). Furthermore, in many cases, DGERFS of MKL did not scale as well as our implementation. We suspect this is due to the need to compute |A|| x| + |b| at each iteration.

The performance of our implementations is sensitive to the tile size nb. Hence, for each matrix dimension n on a different number of cores, we studied the performance of each algorithm with the tile sizes of nb = 80, 160, 240, 320, and 400. We observed that on 48 cores of our AMD machine, the performance is especially sensitive to the tile size, and we tried the additional tile sizes of nb = 340, 360 and 380. In addition, the performance of our incremental pivoting is sensitive to the inner blocking size, and we tried using the block sizes of ib = 10, 20, and 40 for both nb = 80 and 160, ib = 10, 20, 30, 40, 60, and 80 for nb = 240, and ib = 10, 20, 40, and 80 for both nb = 320 and 400. Figure 5 shows the performance obtained using the tile and block sizes that obtained the highest performance of factorization. 1 For the tile algorithms, we included the data layout conversion time as a part of the solution time. We summarize our finding below:

• RBT used the default transformation depth of two and added only a small overhead over no-pivoting in all the test cases.

• In comparison to other pivoting strategies, incremental pivoting could exploit a large number of cores more effectively. As a result, when the performance is not dominated by the trailing submatrix updates (e.g., for a small matrix dimension on 48 cores), it obtained performance that is close to that of nopivoting. However, for a large matrix dimension, due to the special kernels required for the incremental pivoting to update the trailing submatrix, its performance was lower than that of the partial or tournament pivoting LU that uses BLAS-3 DGEMM of MKL for their trailing submatrix updates.

• In comparison to MKL, partial and tournament pivoting LU may reduce the communication and could effectively utilize a larger number of cores. As a result, when the communication impacts the performance (e.g., for medium sized matrices on multiple sockets), they outperformed MKL. We can clearly see this on 48 cores of our AMD machine.

• MKL performed well especially when the trailing submatrix update dominated the performance. Moreover, on the single socket, MKL outperformed no-pivoting for a large enough matrix dimension. This could be because a tiled implementation looses efficiency due to its smaller BLAS kernels used during the trailing submatrix updates.

Clearly, for a large enough matrix size, the computational kernels used for the trailing submatrix update (e.g., DGEMM) dominate the performance. On the other hand, when the communication has a significant impact on the performance (e.g., on a distributed-memory machine), the algorithms that reduce the communication (e.g., RBT and tournament-pivoting) can outperform the standard algorithm. This is implied by our results on 48 cores, and the difference could be greater when the communication becomes more dominant. In our experiments on a shared-memory machine, partial pivoting LU outperformed tournament pivoting in all the cases. However, when the communication latency becomes a significant part of the performance, the tournament pivoting may be favored.

Finally, Figure 6 shows the strong scaling of our implementations. 2In these experiments, for each algorithm, we used the tile size that obtained the best performance on the 48 or 16 cores of the AMD or Intel machine, respectively. Since the matrix dimensions were relatively small in this study, in order to obtain high performance, it becomes imperative to reduce communication and exploit parallelism as much as possible. In particular, RBT and incremental pivoting obtained excellent parallel efficiency.

Accuracy

To study the numerical behavior of our implementations, we used the synthetic matrices from two recent papers [START_REF] Baboulin | Accelerating linear system solutions using randomization techniques[END_REF][START_REF] Grigori | CALU: A communication optimal LU factorization algorithm[END_REF]. In our experiments, all the five pivoting strategies performed similarly on most of the test matrices. In this section, we only present the results of the test matrices which demonstrate the performance characteristics of the pivoting strategies. We also conducted the numerical experiments using all the matrix dimensions used in Section 7.2, but here, we only show the results of n = 30000, which represent the performance trends of the pivoting strategies for all other matrix dimensions. Table 1 shows some properties of these test matrices and the stability results of using partial pivoting. In the table, the second through the sixth test matrices are from the paper by Grigori et al. [START_REF] Grigori | CALU: A communication optimal LU factorization algorithm[END_REF], where the first two have relatively small condition numbers, while the rest are more ill-conditioned. The last three matrices are from the paper by Baboulin et al. [START_REF] Baboulin | Accelerating linear system solutions using randomization techniques[END_REF], where the last test matrix gfpp is one of the pathological matrices that exhibit exponential growth factor using partial pivoting. Since the condition number of the gfpp matrix increases rapidly with the matrix dimension, we used the matrix dimension of m = 1000 for our study. Finally, our incremental and tournament pivoting exhibit different numerical behavior using different tile sizes. For the numerical results presented here, we used the tile and block sizes that obtain the best performance on the 16 core Intel SandyBridge. All the results are in double real precision. • For all the test matrices, tournament pivoting obtained initial backward errors comparable to those of partial pivoting.

• No-pivoting was unstable for five of the test matrices (i.e., ris, fiedler, orthog, {-1 1}, and gfpp). For the rest of the test matrices, the initial backward errors of no-pivoting were significantly greater than those of partial pivoting, but were improved after a few iterative refinements.

• Our incremental pivoting failed for the fiedler, orthog, and gfpp matrices. For other test matrices, its backward errors were greater than those of partial pivoting, but were improved to be in the same order as those of partial pivoting after a few iterative refinements. The only exception was with the ris matrix, where the refinements stagnated before reaching a similar accuracy as that of partial pivoting. In each column of the ris matrix, entries with smaller magnitudes are closer to the diagonal (i.e., A(i, j) = 0.5/(nij + 1.5)). As a result, the permutation matrix P of partial pivoting has ones on the anti-diagonal.

• When no-pivoting was successful, its backward errors were similar to those of RBT. On the other hand, RBT was more stable than no-pivoting, being able to obtain small backward errors for the fiedler, {-1,1}, and gfpp matrices.

• Partial pivoting was not stable for the pathological matrix gfpp. On the other hand, RBT randomizes the original structure of the matrix and was able to compute the solution of reasonable accuracy. It is also possible to construct pathological test matrices, where partial pivoting is unstable while tournament pivoting is stable, and vice versa [START_REF] Grigori | CALU: A communication optimal LU factorization algorithm[END_REF].

Figure 7(b) shows the relative forward error norms of our implementations, which were computed as xx ∞ / x ∞ . We observed similar trends in the convergence of the forward error norms as in that of the backward errors. One difference was with the orthog test matrix, where iterative refinements could not adequately improve the forward errors of incremental pivoting and RBT. Also, even though the backward errors of the ris test matrix were in the order of machine epsilon with partial and tournament pivoting, their relative forward errors were in the order of O(1) due to the large condition number.

CONCLUSIONS

When implemented well, using a fast, recursive panel factorization, tile data layout, and dynamic scheduling, the canonical LU factorization with partial (row) pivoting, is a fast and numerically robust method for solving dense linear systems of equations. On a shared-memory multicore system, its asymptotic performance is very close to the performance of LU factorization without pivoting.

In our experiments on synthetic matrices, tournament pivoting turned out to be as stable as partial pivoting, which has been theoretically proven by its inventors in the first place. It also proved to be fairly fast. However, it failed to deliver on its promise of outperforming partial pivoting, which can be attributed to the shared-memory environment. The method has much more potential for distributed memory systems, where communication matters much more.

Incremental pivoting showed the worst asymptotic performance due to the use of exotic kernels, instead of the GEMM kernel. On the other hand, it showed strong scaling properties almost as good as RBT and no-pivoting. It is harder to make strong claims about its numerical properties. Its initial residual is usually worse than that of partial and tournament pivoting, but in most cases the accuracy is quickly recovered in iterative refinement. It can fail in some situations, when partial and tournament pivoting prevail.

RBT is the fastest method, both asymptotically, and in terms of strong scaling, because it only adds a small overhead of preprocessing and postprocessing to the time of factorization without pivoting. Similarly to incremental pivoting, it produces a high initial residual, but the accuracy can be recovered in iterative refinement. Similarly to incremental pivoting, it can fail in situations, when partial pivoting and tournament pivoting prevail.

And finally, it can be observed that iterative refinement is a powerful mechanism of minimizing the backward error, which in most cases translates to minimizing the forward error.

FUTURE DIRECTIONS

Although we believe that the wide range of synthetic matrices with different properties gives a good understanding of numerical properties of the different flavors of the LU factorization, ultimately it would be invaluable to make such comparison using matrices coming from real world applications, such as plasma burning or radar cross section.

Although we believe that the comparison using relatively large shared memory systems, by todays standards, gives a good insight into the performance properties of the different LU factorization algorithms, we acknowledge that the picture can be very different in a distributed memory environment. Ultimately we would like to produce a similar comparison using distributed memory systems.

Figure 1 :

 1 Figure 1: The block LU factorization (Level 3 BLAS algorithm of LAPACK).

Figure 2 :

 2 Figure 2: Incremental LU factorization.

Figure 3 :

 3 Figure 3: Example of panel factorization using TSLU with P = 4 processors. A binary tree is used for the reduction operation

Figure 4 :

 4 Figure 4: The Column-Column Rectangular Block (CCRB) matrix layout.

Figure 5 :

 5 Figure 5: Asymptotic performance comparison of LU factorization algorithms.

  AMD Opteron (m = n = 4000).

  Intel SandyBridge (m = n = 2000).

Figure 6 :

 6 Figure 6: Strong scaling comparison of LU factorization algorithms.

Figure 7 (

 7 Figure 7(a) shows the component-wise backward errors, max i |r i |/(|A|| x| + |b|) i , at each step of the iterative refinements. For these experiments, the right-hand-side b is chosen such that the entries of the exact solution x are uniformly distributed random numbers in the range of [-0.5, 0.5]. Below, we summarize our findings:

random circul riemann ris compan fiedler orthog {- 1

 1

Figure 7 :

 7 Figure 7: Numerical accuracy of LU factorization algorithms, showing error before refinement (top point of each line) and after each refinement iteration (subsequent points).

Table 1 :

 1 Properties of test matrices and stability results of using partial pivoting (n = m = 30000).

	name	description	A 1	cond(A, 2)	L 1	L -1 1	max |U(i, j)| max |U(i, i)|	U 1	cond(U, 1)
	random	dlarnv(2)	7.59e+03	4.78e+05	1.50e+04 8.60e+03	1.54e+02	1.19e+02	2.96e+05	2.43e+09
	circul	gallery('circul', 1 : n)	2.43e+04	6.97e+02	6.66e+03 8.64e+03	5.08e+03	3.87e+03	1.50e+06	4.23e+07
	riemann gallery('riemann', n)	1.42e+05	3.15e+05	3.00e+04 3.50e+00	3.00e+04	3.00e+04	2.24e+05	1.25e+08
	ris	gallery('ris', n)	1.16e+01	3.34e+15	2.09e+04 3.43e+02	7.34e+00	3.30e+00	3.46e+02	1.42e+21
	compan	compan(dlarnv(3))	4.39e+00	1.98e+04	2.00e+00 1.01e+01	3.39e+00	1.85e+00	1.90e+01	8.60e+01
	fiedler gallery('fiedler', 1 : n)	1.50e+04	1.92e+09	1.50e+04 1.37e+04	2.00e+00	1.99e+00	2.71e+04	9.33e+09
	orthog	gallery('orthog', n)	1.56e+02	1.00e+00	1.91e+04 1.70e+03	1.57e+03	1.57e+02	2.81e+03	3.84e+08
	{-1,1} gfpp †	3.00e+04 gfpp(triu(rand(n)),1e-4) 1.00e+03 A(i, j) = -1 or 1	1.81e+05 1.42e+19	3.00e+04 8.67e+03 9.02e+02 2.10e+02	5.47e+03 4.98e+00	3.78e+02 2.55e+00	1.00e+06 4.28e+02	8.35e+08 5.96e+90
	† For gfpp, n = 1000.

We computed the Gflop/s as the ratio of the number of flops required for LU factorization and a pair of forward and backward substitutions

,

[START_REF] Agullo | Numerical linear algebra on emerging architectures: The PLASMA and MAGMA projects[END_REF] 3 n 3 + 2n 2 , over the total time required for the factorization and iterative refinement.

Using the default transformation depth of two, our current implementation of RBT assumes the matrix dimension to be a multiple of four times the tile size. Hence, for these experiments, in order to use the same block sizes as those used for no-pivoting, we set the matrix dimensions to be m = 1920 and m

= 3840 on the AMD and Intel machines, respectively, for RBT.

* This work was supported in part by the U. S. Department of Energy, the National Science Foundation, and the Intel Corporation.


Efi Fogel

Monique Teillaud

The Computational Geometry Algorithms Library CGAL *

is an open source software library that provides industrial and academic users with easy access to reliable implementations of efficient geometric algorithms.

Usage. CGAL is used in a diverse range of domains requiring geometric computation such as computer graphics, scientific visualization, computer aided design and modeling, geographic information systems, molecular biology, medical imaging, and many more. Since CGAL provides a wide range of components, we restrict ourselves to mentioning just a few here. As an example application of CGAL, a series of packages are provided which are useful in robotics and automation: Minkowski sums, offset polygons, Boolean operations on curved regions. The high precision of CGAL allows users to solve geometric problems involving motion in restricted environments, such as those arising in assembly planning. The robustness and efficiency of components such as the Delaunay triangulation and mesh construction and manipulation packages makes CGAL attractive for simulations, in particular those involving proteins, particle physics, fluid dynamics, medical modeling, biophysics, geophysics, and astronomy. Indeed, the aforementioned components are largely used in these areas. Some support for manipulations of polynomials and for solving univariate polynomial equations and bivariate polynomial systems is also provided, as well as handling for convex quadratic programs.

History of the CGAL Open Source Project. Several European research groups started to develop their own small geometry libraries in the early 90's. In 1996, a consortium of eight sites was created to gather the work of these groups into a single software library, namely CGAL. Their main goal was to promote research in computational geometry and to translate the results into robust software suitable for industrial applications. Around this time the Computational Geometry Impact Task Force Report [C + 96, C + 99] made a series of recommendations. Amongst these recommendations, the production and distribution of usable (and useful) geometric software, and the need to establish a reward structure for software implementations in academia, were key. On November 2003, when version 3.0 was released, CGAL officially became an Open Source project, allowing new contributors to join the project.

License. CGAL is distributed under the GPL license (apart from a few basic parts, which are distributed under the LGPL license). In particular, it is publicly and freely available for academic use. Commercial licenses are offered by Geometry Factory, a company founded in 2003 mainly for this purpose.

Editorial board. The CGAL editorial board was created in 2001. It currently consists of thirteen members. The main task of the editorial board is to assure the quality of CGAL. It is also responsible for making decisions about technical matters and coordinating communication and promotion of CGAL. All new packages must be submitted to the Editorial Board to be reviewed before they can be accepted and integrated into the library, in a process that is very similar to the standard review process for papers published in conference proceedings or journals. More information about the submission process is available at http://www.cgal.org/review_process_rules.html.

Style and Techniques.

CGAL is a unique library both in general and within the field of computational geometry in particular, as it consists of a large number of components with a homogeneous API (Application Programming Interface). Careful choices in design and programming style have made CGAL the de facto standard in the field of applied computational geometry. Its development started whilst the standardization process of C++ and the STL (Standard Template Library) was taking place. Indeed, the programming style is very close to the programming style of STL; it rigorously adheres to the generic programming paradigm-a discipline that consists of the gradual lifting of concrete algorithms abstracting over details, while retaining the algorithm semantics and efficiency. The programming style of CGAL also facilitates the process of interfacing with third party software. Each package comes with header files consisting not only of the interface, but also the generic implementation of the package code, comprehensive and didactic on-line documentation, a set of non-interactive standalone example programs, and an optional interactive demo with a graphical user interface.

Robustness. CGAL follows the exact geometric-computation paradigm, which simply amounts to ensuring that errors in predicate evaluations do not occur; it guarantees robustness of the applied algorithms. We additionally remark that every package also includes a collection of function and regression test. 1 The tests provided by each package are combined into one place to form the CGAL test suite. This test suite is run daily and its results are automatically assembled, analyzed, and reported.

Impact. Measuring the impact of software is a difficult task, especially in the Open Source software community. Even if some hard numbers can be found, they can be difficult to interpret. The following facts may shed some light on the impact of CGAL:

• There are roughly 1000 downloads per month from http://gforge.inria.fr/

• CGAL is included in various software distribution channels, such as Fedora, Debian/Ubuntu, and Macports.

• The range of uses of CGAL is very broad, as shown by the sample list of projects using CGAL, which is available at http://www.cgal.org/projects.html. In addition, many projects shown in http://acg.cs.tau.ac.il/projects use CGAL or even describe the development of a CGAL component.

• The CGAL triangulation packages were integrated in Matlab 2009a.2

• Springer has published a book entitled "CGAL Arrangements and Their Applications" authored by some of the developers of the 2D Arrangements package and its derivatives.

• Concerning the public mailing lists, there are currently -4000 subscribers to the announcement list cgal-announce@lists-sop.inria.fr -1500 to the public discussion list cgal-discuss@lists-sop.inria.fr, with high traffic: users are free to ask questions, which are often rapidly answered by the developers or other users.

However these are not distributed as part of the public releases.

Watch the video at http://www.mathworks.com/products/demos/shipping/matlab/ New-MATLAB-Mathematics-Features-in-R2009a.html

Acknowledgments

We thank Ross Hemsley for helping us distilling the text.

-sop.inria.fr/members/Monique.Teillaud/