Daniël M J Peeters
email: d.m.j.peeters@tudelft.nl

Franc ¸ois-Xavier Irisarri

Mostafa M Abdalla

OPTIMISING THE PLY DROPPING ORDER IN VARIABLE STIFFNESS, VARIABLE THICKNESS LAMINATES USING STACKING SEQUENCE TABLES

Keywords: optimisation, stacking sequence table, variable thickness, variable stiffness

Modern composite structures offer two avenues of optimising performance. One, optimising a single stacking sequence over the structure leading to constant stiffness designs. Two, varying the stiffness over the structure. This may be achieved by dropping plies, changing the thickness, or by steering the fibres, changing the fibre angles. Optimising ply drops involves two decisions: the selection of ply drop boundaries, and of the ply drop order. Previous work considered the problems of optimising the fibre angle distribution and ply drop boundaries but ply drop order was pre-specified. This paper extends the work to simultaneously optimise the fibre angle distribution, the ply drop boundaries, and order. The optimisation of fibre angle distribution lends itself to gradient-based methods. The ply drop boundary optimisation is formulated using topology optimisation techniques and is thus also solvable using gradient-based methods. The ply drop order optimisation requires discrete variables and is hence approached using an evolutionary algorithm based on stacking sequence tables. In this paper an efficient multi-step algorithm is developed combining the optimisation of all aspects of variable stiffness laminates. The results indicate that significant improvements may be obtained by including the ply drop order in the optimisation at a relatively modest computational cost.

Introduction

Composite materials are attractive due to their high stiffness-to-weight and strength-to-weight ratio. It has been shown that by spatially varying the stiffness, even better performance can be obtained without adding extra weight. Varying the stiffness can be done in two ways: either by changing the fibre angles, by steering the fibres, or by changing the number of plies from one point to the next by dropping plies. To develop constant thickness, steered laminates, a three-step optimisation approach has been developed. [START_REF] Ijsselmuiden | Optimal design of variable stffness composite structures using lamination parameters[END_REF][START_REF] Samuel | Optimization of variable-stiffness panels for maximum buckling load using lamination parameters[END_REF] In the first step the optimal stiffness distribution in terms of lamination parameters is found, in the second step the optimal fibre angles are obtained, and in the third step the optimal fibre paths are retrieved. A thickness variation has already been implemented in the first step of the optimisation, and showed that significant improvements could be obtained, however, no information about fibre angle or number of plies is available at the first step, so the physical construction of the laminate remains to be found. A method to optimise the fibre angle and ply drop locations has recently been developed, but the ply drop order was pre-specified. [START_REF] Peeters | Optimization of ply drop locations in variable-stiffness composites[END_REF] ECCM17 -17 th European Conference on Composite Materials Munich, Germany, 26-30 th June 2016 2

Thickness variation in a laminate is described using two variables: the ply drop location and the ply drop order. The most popular approach is to use an evolutionary algorithm, typically a genetic algorithm, to optimise the number of layers per 'patch', while also optimising the ply drop order and stacking sequence, limited to a discrete set of angles (e.g., 0 • , ±45 • , and 90 • ). This area of research is referred to as laminate blending [START_REF] Kim | Optimum design of composite structures with ply drop using genetic algorithm and expert system shell[END_REF][START_REF]Genetic algorithm optimization and blending of composite laminates by locally reducing laminate thickness[END_REF][6][START_REF] Sørenn | Global blending optimization of laminated composites with discrete material candidate selection and thickness variation[END_REF][START_REF] Liu | Bilevel optimization of blended composite wing panels[END_REF] and assumes that potential ply drop locations (i.e., patch boundaries) are pre-specified by the user. A technique where the fibre angle is not restricted to a discrete set has also been developed, however, no manufacturing constraints, limiting the change in fibre angle from one element to the next, are posed. [START_REF] De Leon | Aeroelastic tailoring using fiber orientation and topology optimization[END_REF] Other techniques where the ply drop locations are not pre-specified use continuous optimisation. Shape optimisation [START_REF] Delgado | Optimization of composite structures: A shape and topology sensitivity analysis[END_REF] is used to determine the shape and hence ply coverage and ply drop locations of the different layers. The optimisation is performed using a level-set approach with fibre angles limited to a discrete set. Another continuous method is the discrete material and thickness optimisation method, where the fibre angles belong to a discrete set and fictitious density variables are used to select the ply angles at any given location. This has been done for compliance and buckling optimisation. [START_REF] Søren | Dmto a method for discrete material and thickness optimization of laminated composite structures[END_REF][START_REF] Lund | Buckling topology optimization of laminated multi-material composite shell structures[END_REF] For this method, also a thickness filter has been implemented to get to physically feasible designs. [START_REF] Sørensen | Thickness filters for gradient based multi-material and thickness optimization of laminated composite structures[END_REF] Thickness optimisation for buckling load under uni-and bi-axial compression has also been performed. This work showed that large improvement in buckling load could be made without affecting the in-plane stiffness. [START_REF] Joshi | Thickness optimization for maximum buckling loads in composite laminated plates[END_REF] The easiest ply drop orders are inner or outer blending, where the layers are dropped from the symmetry plane, or from the outside respectively. To determine the optimal ply drop order, guide-based designs can be used: [START_REF]Genetic algorithm optimization and blending of composite laminates by locally reducing laminate thickness[END_REF] a stacking sequence for the thickest laminate, called the guide laminate, is defined and the number of layers per patch. The stacking sequence is then derived by dropping layers from the inside or outside, depending whether inner or outer blending is used, from the guide laminate. A method that offers more possible ply drop orders and takes into account industrial guidelines is using stacking sequence tables. [START_REF] Irisarri | Optimal design of laminated composite structures with ply drops using stacking sequence tables[END_REF] A ply drop order and guide laminate are optimised.

This paper aims to combine the ideas behind the optimisation of optimising ply drop locations and the stacking sequence tables. [START_REF] Peeters | Optimization of ply drop locations in variable-stiffness composites[END_REF][START_REF] Irisarri | Optimal design of laminated composite structures with ply drops using stacking sequence tables[END_REF] The outcome will be an optimal steered, variable thickness laminate: both the ply drop location and ply drop order are optimised. The remainder of the paper is organised as follows: first the general approach is explained in section 2, next a description of the stacking sequence tables is given in section 3, some numerical examples in section 4 and finally a conclusion in section 5.

Approach

The overall optimisation strategy is based on the three step optimisation approach developped by IJsselmuiden [START_REF] Ijsselmuiden | Optimal design of variable stffness composite structures using lamination parameters[END_REF]. The present work focuses on the second step of the optimisation and combines the optimisation of the spatial distribution of the fiber angles and the optimisation of the ply drop location and ply drop order. Figure 1 gives an overview of the proposed optimisation approach.

Step 1 (see [START_REF] Samuel | Optimization of variable-stiffness panels for maximum buckling load using lamination parameters[END_REF][START_REF]Design of variable-stiffness composite panels for maximum buckling load[END_REF][START_REF] Khani | Design of variable stiffness panels for maximum strength using lamination parameters[END_REF]) returns an idealized design defined as an optimal stiffness and thickness distribution over the structure.

The idealized design gives an upper bound of the performance of the structure.

In the proposed approach, Step 2 is subdivided into two successive optimisation phases referred to as

Step 2.1 and and Step 2.2 respectively. Both phases combine an evolutionary optimiser and a gradientbased optimizer for their complementary efficiency in solving combinatorial problems and continuous problems respectively.

Step 2.1 aims at providing a relevant initial guess of the fibre angle distributions per ply, total thickness distribution and ply drop order for the subsequent optimisation phase.

Step 2.1 takes as input the ide- Step 2.2 aims at converting the laminate thickness distributions into properly defined ply drop locations. The optimisation is initialized with the non-dominated front issued from Step 2.1. An evolutionary algorithm (EA) specialized for ply drop order optimisation is hybridized with a gradient-based method devised for ply drop location optimisation. For each ply drop order generated by the EA, a topology-like optimisation is performed using a fictitious density distribution for each ply. The densities are forced to converge to either one or zero which defines the ply coverage and the ply drop locations. The gradient based optimisation alternates between density optimisation and fiber angle optimisation in each ply. The optimisation of the ply drop locations is not discussed in this paper, the interested reader is referred to Peeters and Abdalla [START_REF] Peeters | Optimization of ply drop locations in variable-stiffness composites[END_REF].

Stacking sequence table optimisation

The Evolutionary Algorithms (EA) devised in the present work for Step 2.1 and Step 2.2 (see Figure 1) share the same architecture based on the algorithm proposed in [START_REF] Irisarri | Optimal design of laminated composite structures with ply drops using stacking sequence tables[END_REF]. The main features of the algorithm that are preserved in the present work are briefly summarized in the following. The algorithm uses two populations: a regular population and an archive population. The archive population contains a fixedsized subset of the best solutions encountered during the evolutionary process. The optimization starts with an initial population and an empty archive. The following steps are then iteratively repeated. • The solutions composing the current populations are evaluated.

• The archive is updated using the environmental selection method proposed in [START_REF] Zitzler | SPEA2: Improving the Strength Pareto Evolutionary Algorithm[END_REF].

• Solution fitness values are assigned to the union of the current population and the archive population. • The fitness function aggregates a density estimate and the Pareto-rank of the solution to enable multi-objective optimization. • Optimization constraints formulated on the mechanical performances of the solutions are enforced through a modified binary tournament selection scheme [START_REF] Deb | A constraint handling strategy for robust multi-criterion optimization[END_REF].

• The current population is replaced after mutation and recombination of the selected solutions.

• Design guidelines whose evaluation require no mechanical analysis of the solutions, for example symmetry and balance of the laminates, are enforced by construction of the solutions, using specific variation operators.

Step 2.1. SST-based stiffness matching

A Stacking Sequence Table (SST) describes the sequence of ply insertions ensuring the transition between a thin laminate and a thicker one. The SST is encoded within the EA using two integer vectors.

The vector S S T lam represents the stacking sequence of the thickest laminate. The S S T ins corresponds to the rank of insertion of the plies, and reciprocally defines the ply drop sequence. Rank 0 is attributed to the plies that form the thinner laminate for the SST. The first ply inserted is given rank 1, the second ply inserted is given rank 2 and so on. For example, if S S T ins = [0 3 6 1 5 7 2 0 4] then the ply drop order is [6/3/5/9/2/7/4/8/1]. In the present work, laminates are symmetrical and composed of balanced pairs of plies which allows to encode only one fourth of the SST. Additionally, the covering guideline and internal continuity guideline are enforced during the optimization (see [START_REF] Irisarri | Optimal design of laminated composite structures with ply drops using stacking sequence tables[END_REF]) for strength concerns. The covering guideline imposes that the plies on the surface of the laminates are never dropped in the SST. The internal continuity guideline states that a continuous ply (rank = rank min = 0 in S S T ins ) should be kept every three consecutive dropped plies (rank > 0) . Here, these guidelines are evaluated on the vectors S S T ins , thus the guidelines are transposed to pairs of plies. The internal continuity guideline implies a dependency of the maximal number of plies n max to the minimal number of plies n min within the SST:

n max ≤ 3 × n min -2.
When reached, the limitation is relaxed by incrementing the value of rank min . For instance, Figure 2 shows a unidirectional thickness transition between a 12-ply laminate and a 36-ply laminate compatible with the following encoding: S S T lam = [45 30 60 0 -75 -15 -60 45 90] and S S T ins = [0 3 6 1 5 7 2 0 4].

For stiffness matching the thickness distribution in the structure is set to the thickness distribution of the idealized design rounded to the nearest upper number of plies. Subsequently, the thickness distribution is kept unchanged. The stiffness distribution of the idealized design is targeted. Proximity in stiffness space is evaluated using the distance proposed in [START_REF] Irisarri | Improved Shepard's Method for the Optimization of Composite Structures[END_REF]. The optimization aims at minimizing both the A-distance in membrane space and the D-distance in bending space as defined by Equation 1:

d(C 1 , C 2 ) = C -1 1 : C 2 + C 1 : C -1 2 -6. ( 1
)
where the operator : stands for the trace of the product of two square matrices and C 1 and C 2 are either two membrane stiffness matrices or two bending stiffness matrices. The objective functions correspond to the sum of the distances computed at the nodes and weighted by the nodal areas. The optimization returns a set of non-dominated solutions. Each solution represent a trade-off between proximity to the membrane stiffness distribution of the idealized design and its bending stiffness distribution.

Step 2.2. Ply-drop order optimization

At Step 2.2, the design variables are distributed between the overall evolutionary algorithm and the gradient-based solver used during the ply drop location optimiser. [START_REF] Peeters | Optimization of ply drop locations in variable-stiffness composites[END_REF] The EA is used to solve the combinatorial optimization problem of the ply drop order and calls for each solution the gradient-based solver to perform local improvement of the topology and fiber steering of the plies. Thus, the EA operates on the ply drop sequence only, which is encoded by the vector S S T ins . The permutation operator for the vector S S T ins proposed in [START_REF] Irisarri | Optimal design of laminated composite structures with ply drops using stacking sequence tables[END_REF] is reused in the present work. The operator permutes the insertion ranks of two plies within the vector S S T ins for plies with non-zero ranks. A specific crossover operator is devised to operate on insertion rank vectors S S T ins , inspired by the ordered crossover operator [START_REF] Goldberg | Genetic algorithms in search, optimization, and machine learning[END_REF]:

1. Let S S T 01 ins and S S T 02 ins be two insertion rank vectors. 2. A random subset of ranks strictly superior to rank min is selected. 3. In the proposed example rank min = 1, and the selected ranks are underlined: S S T 01 ins = [0 2 3 5 1 6 7 0 4] and S S T 02 ins = [0 6 7 1 4 5 2 0 3]. 4. To create the offspring vector S S T 11 ins , vector S S T 01 ins is copied and the selected subset is reordered according to S S T 02 ins . Vector S S T 12 ins similarly is a copy of S S T 02 ins with the ordering of S S T 01 ins for the selected subset. As it does not operate on ranks lower than rank min , the proposed crossover operator preserves the internal continuity guideline, as defined in Section 3.1.

Results

The proposed method is applied to the buckling load maximisation of a simply supported square plate under longitudinal compression. The first two buckling modes are considered in the optimisation. Furthermore, a stiffness constraint is used: all solutions are required to present higher or equal longitudinal stiffness than the 24-ply quasi-isotropic constant thickness laminate. The side length of the plate is 500mm. The base ply properties are the following: E 11 = 177GPa, E 22 = 10.8GPa, G 12 = 7.6GPa and ν 12 = 0.27. For the optimisation, the number of plies can locally vary from 8 plies to 36 plies. The maximum volume allowed corresponds to a 24-ply constant thickness laminate. A global steering constraint is enforced that corresponds to a minimal steering radius of 333 mm. The performances of the best solution found are detailed in Table 1. The corresponding ply insertion chromosome S S T ins is [0 5 4 1 7 2 3 0 6]. The fibre paths are presented in Figure 3 for each design layer composed of a ply and its adjacent balanced counterpart. The solution forms large stiffeners along the longitudinal sides of the plate to increase the in-plane stiffness.

In order to quantify the influence of the ply drop order on the buckling performance of the plate, a buckling minimization was run at Step 2.2 with the same initialization of the EA. The minimal critical buckling load obtained 2 is 4.449 (i.e., 117.8% improvement with respect to the best constant thickness design). Hence 15% lower than the best buckling load found.

Conclusion

In this paper, a method to optimise the fibre angle distribution, ply drop location and ply drop order is presented. The optimisation is a combination of an evolutionary and gradient based optimiser. The ply drop order is optimised using a permutation genetic algorithm, while the ply drop locations and fibre angle distributions are optimised using the method of successive approximations. The ply drop locations are optimised using topology-like ideas: each layer is given a fictitious density distribution and this is optimised in two steps: in the first step the overall thickness is optimised, in the second step the exact ply drop locations are determined at places where the fictitious density changes from near one to near zero.

Initial results indicate that large improvements in performance are possible by also optimising the drop order. A flat, square plate under uni-axial compression has been optimised for buckling with stiffness constraint. The stiffness is constrained to be at least the stiffness of a quasi-isotropic laminate with the same weight without ply drops. The influence of the dropping order on the buckling load is large: up to 32.8% of the buckling load of the best steered constant stiffness laminate. These results show that checking just one possible drop order during the ply drop location optimisation can have a significant influence on the final result.

While for the fibre angle distribution some manufacturing constraints, such as a maximum change in fibre angle between adjacent nodes, is taken into account, the effect of gaps and overlaps occurring due to the steering is not taken into account. Furthermore, the minimal distance between consecutive ply drops is not constrained, meaning multiple ply drops can occur at the same point, or within a short distance from each other. While this may be the subject of future research, a post-processing step can also be implemented to adhere to this minimal distance.

ECCM17 - 17 thFigure 1 .

 171 Figure 1. Overview of the optimisation strategy.

Figure 2 .

 2 Figure 2. Schematic of a unidirectional symmetrical thickness transition.

ECCM17 - 17 thFigure 3 .

 173 Figure 3. Outcome of Step 2.2. Ply fiber orientation shown with their balance counterpart.

Table 1

 1 shows shows the highest buckling loads obtained after optimisation Step 1, Step 2.1 and Step 2.2 of the proposed method. The results are compared with two constant thickness solutions: the constant thickness with straight fibers, and the steered fibre laminate. The steered fiber variable thickness solution was obtained in[START_REF] Peeters | Optimization of ply drop locations in variable-stiffness composites[END_REF] with S S T ins defined as [0 6 4 2 1 3 7 5 0]. All buckling loads are normalized with respect to the critical buckling load of the 24-ply quasi-isotropic "black metal" solution (i.e., all lamination parameters equal to zero). The reference solution to assess the performance improvement is the steered fibre constant thickness solution. The idealized design obtained at Step 1 is defined as

	ECCM17 -17 th European Conference on Composite Materials
	Munich, Germany, 26-30 th June 2016

Table 1 .

 1 Overview of the highest buckling load obtained after each optimisation step, compared with results presented in[START_REF] Peeters | Optimization of ply drop locations in variable-stiffness composites[END_REF] and constant thickness solutions. a combined stiffness (i.e., lamination parameters) and thickness distribution and gives an upper bound on the achievable performances. The constant thickness with straight fibers solution is the following: [±41 ± 41 ± 40 ± 35 ± 19 ± 6] s .

Daniël M.J. Peeters, Franc ¸ois-Xavier Irisarri and Mostafa M. Abdalla

Acknowledgements

This work is supported by the CANAL (CreAting Non-conventionAl Laminates) project, part of the European Union Seventh Framework Program.


Massimiliano Fasi
email: massimiliano.fasi@postgrad.manchester.ac.uk

Julien Langou
email: julien.langou@ucdenver.edu

Yves Robert
email: yves.robert@inria.fr

Bora Uçar
email: bora.ucar@ens-lyon.fr

A Backward/Forward Recovery Approach for the Preconditioned Conjugate Gradient Method $

Keywords: Fault-tolerance, Silent errors, Algorithm-based fault tolerance, Checkpointing, Sparse matrix-vector multiplication, Preconditioned conjugate gradient method. 2010 MSC: 68M15, 94C12, 65F08

Several recent papers have introduced a periodic verification mechanism to detect silent errors in iterative solvers. Chen [PPoPP'13, has shown how to combine such a verification mechanism (a stability test checking the orthogonality of two vectors and recomputing the residual) with checkpointing: the idea is to verify every d iterations, and to checkpoint every c × d iterations. When a silent error is detected by the verification mechanism, one can rollback to and re-execute from the last checkpoint. In this paper, we also propose to combine checkpointing and verification, but we use algorithm-based fault tolerance (ABFT) rather than stability tests. ABFT can be used for error detection, but also for error detection and correction, allowing a forward recovery (and no rollback nor re-execution) when a single error is detected. We introduce an abstract performance model to compute the performance of all schemes, and we instantiate it using the preconditioned conjugate gradient algorithm. Finally, we validate our new approach through a set of simulations.

Introduction

Silent errors (or silent data corruptions) have become a significant concern in HPC environments [START_REF] Moody | Design, modeling, and evaluation of a scalable multi-level checkpointing system[END_REF]. There are many sources of silent errors, from bit flips in cache caused by cosmic radiations, to wrong results produced by the arithmetic logic unit. The latter source becomes relevant when the computation is performed in the low voltage mode to reduce the energy consumption in large-scale computations. But the low levels of voltage dramatically reduces the dependability of the system.

The key problem with silent errors is the detection latency: when a silent error strikes, the corrupted data is not identified immediately, but instead only when some numerical anomaly is detected in the application behavior. It is clear that this detection can occur with an arbitrary delay. As a consequence, the de facto standard method for resilience, namely checkpointing and recovery, cannot be used directly. Indeed, the method of checkpointing and recovery applies to fail-stop errors (e.g., hardware crashes): such errors are detected immediately, and one can safely recover from the last saved snapshot of the application state. On the contrary, because of the detection latency induced by silent errors, it is often impossible to know when the error struck, and hence to determine which checkpoint (if any) is valid to safely restore the application state. Even if an unlimited number of checkpoints could be kept in memory, there would remain the problem of identifying a valid one.

In the absence of a resilience method, the only known remedy to silent errors is to re-execute the application from scratch as soon as a silent error is detected. On large-scale systems, the silent error rate grows linearly with the number of components, and several silent errors are expected to strike during the execution of a typical large-scale HPC application [START_REF] Cappello | Fault tolerance in petascale/exascale systems: Current knowledge, challenges and research opportunities[END_REF][START_REF] Cappello | Toward exascale resilience[END_REF][START_REF] Cappello | Toward exascale resilience: 2014 update[END_REF][START_REF] Schroeder | Understanding failures in petascale computers[END_REF]. The cost of re-executing the application one or more times becomes prohibitive, and other approaches need to be considered.

Several recent papers have proposed to introduce a verification mechanism to be applied periodically in order to detect silent errors. These papers mostly target iterative methods to solve sparse linear systems, which are natural candidates to periodic detection. If we apply the verification mechanism every, say, d iterations, then we have the opportunity to detect the error earlier, namely at most d -1 iterations after the actual faulty iteration, thereby stopping the progress of a flawed execution much earlier than without detection. However, the cost of the verification may be non-negligible in front of the cost of one iteration of the application, hence the need to trade off for an adequate value of d. Verification can consist in testing the orthogonality of two vectors (cheap) or recomputing the residual (cost of a sparse matrix-vector product, more expensive). We survey several verification mechanisms in Section 2. Note that in all these approaches a selective reliability model is enforced, where the parts of the application that are not protected are assumed to execute in a reliable mode.

While verification mechanisms speed up the detection of silent errors, they cannot provide correction, and thus they cannot avoid the re-execution of the application from scratch. A solution is to combine checkpointing with verification. If we apply the verification mechanism every d iterations, we can checkpoint every c × d iterations, thereby limiting the amount of re-execution considerably. A checkpoint is always valid because it is being preceded by a verification. If an error occurs, it will be detected by one of the c verifications performed before the next checkpoint. This is exactly the approach proposed by Chen [START_REF] Chen | Online-ABFT: An Online Algorithm Based Fault Tolerance Scheme for Soft Error Detection in Iterative Methods[END_REF] for a variety of methods based on Krylov subspaces, including the widely used conjugate Gradient (CG) algorithm. Chen [START_REF] Chen | Online-ABFT: An Online Algorithm Based Fault Tolerance Scheme for Soft Error Detection in Iterative Methods[END_REF] gives an equation for the overhead incurred by checkpointing and verification, and determines the best values of c and d by finding a numerical solution of the equation. In fact, computing the optimal verification and checkpoint intervals is a hard problem. In the case of pure periodic checkpointing, closed-form approximations of the optimal period have been given by Young [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF] and Daly [START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF]. However, when combining checkpointing and verification, the complexity of the problem grows. To the best of our knowledge, there is no known closed-form formula, although a dynamic programming algorithm to compute the optimal repartition of checkpoints and verifications is available [START_REF] Benoit | Assessing general-purpose algorithms to cope with fail-stop and silent errors[END_REF].

For linear algebra kernels, another widely used technique for silent error detection is algorithm-based fault tolerance (ABFT). The pioneering paper of Huang and Abraham [START_REF] Huang | Algorithm-Based Fault Tolerance for Matrix Operations[END_REF] describes an algorithm capable of detecting and correcting a single silent error striking a dense matrix-matrix multiplication by means of row and column checksums. ABFT protection has been successfully applied to dense LU [START_REF] Du | Algorithm-based fault tolerance for dense matrix factorizations[END_REF], LU with partial pivoting [START_REF] Yao | Detection of soft errors in LU decomposition with partial pivoting using algorithm-based fault tolerance[END_REF], Cholesky [START_REF] Hakkarinen | Fail-stop failure algorithm-based fault tolerance for cholesky decomposition, Parallel and Distributed Systems[END_REF] and QR [START_REF] Du | Soft error resilient QR factorization for hybrid system with GPGPU[END_REF] factorizations, and more recently to sparse kernels like SpMxV (matrix-vector product) and triangular solve [START_REF] Shantharam | Fault Tolerant Preconditioned Conjugate Gradient for Sparse Linear System Solution[END_REF]. The overhead induced by ABFT is usually small, which makes it a good candidate for error detection at each iteration of the CG algorithm.

The beauty of ABFT is that it can correct errors in addition to detecting them. This comes at the price of an increased overhead, because several checksums are needed to detect and correct, while a single checksum is enough when just detection is required. Still, being able to correct a silent error on the fly allows for forward recovery. No rollback, recovery nor re-execution are needed when a single silent error is detected at some iteration, because ABFT can correct it, and the execution can be safely resumed from that very same iteration. Only when two or more silent errors strike within an iteration we do need to rollback to the last checkpoint. In many practical situations, it is reasonable to expect no more than one error per iteration, which means that most roll-back operations can be avoided. In turn, this leads to less frequent checkpoints, and hence less overhead.

The major contributions of this paper are an ABFT framework to detect multiple errors striking the computation and a performance model that allows to compare methods that combine verification and checkpointing. The verification mechanism is capable of error detection, or of both error detection and correction. The model tries to determine the optimal intervals for verification and checkpointing, given the cost of an iteration, the overhead associated to verification, checkpoint and recovery, and the rate of silent errors. Our abstract model provides the optimal answer to this question, as a function of the cost of all application and resilience parameters.

We instantiate the model using a CG kernel, preconditioned with a sparse approximate inverse [START_REF] Benzi | A sparse approximate inverse preconditioner for nonsymmetric linear systems[END_REF], and compare the performance of two ABFT-based verification mechanisms. We call the first scheme, capable of error detection only, ABFT-Detection and the second scheme, which enhances the first one by providing single error correction as well, ABFT-Correction. Through numerical simulations, we compare the performance of both schemes with Online-Detection, the approach of Chen [START_REF] Chen | Online-ABFT: An Online Algorithm Based Fault Tolerance Scheme for Soft Error Detection in Iterative Methods[END_REF] (which we extend to recover from memory errors by checkpointing the sparse matrix in addition to the current iteration vectors). These simulations show that ABFT-Correction outperforms both Online-Detection and ABFT-Detection for a wide range of fault rates, thereby demonstrating that combining checkpointing with ABFT correcting techniques is more efficient than pure checkpointing for most practical situations.

Our discussion focuses on the sequential execution of iterative methods. Yet, all our techniques extend to parallel implementation based on the message passing paradigm (with using, e.g., MPI). In an implementation of SpMxV in such a setting, the processing elements (or processors) hold a part of the matrix and the input vector, and hold a part of the output vector at the end. A recent exposition of different algorithms can be found elsewhere [START_REF] Kaya | Analysis of partitioning models and metrics in parallel sparse matrix-vector multiplication[END_REF]. Typically, the processors perform scalar multiply operations on the local matrix and the input vector elements, when all required vector elements have been received from other processors. The implementations of the MPI standard guarantees correct message delivery, i.e., checksums are incorporated into the message so as to prevent transmission errors (the receives can be done in-place and hence are protected). However, the receiver will obviously get corrupted data if the sender sends corrupted data. Silent errors can indeed strike at a given processor during local scalar multiply operations. Performing error detection and correction locally implies global error detection and correction for the SpMxV. Note that, in this case, the local matrix elements can form a matrix which cannot be assumed to be square in general (for some iterative solvers they can be). Furthermore, the mean time between failures (MTBF) reduces linearly with the number of processors. This is well-known for memoryless distributions of fault inter-arrival times and remains true for arbitrary continuous distributions of finite mean [START_REF] Aupy | Checkpointing algorithms and fault prediction[END_REF]. Therefore, resilient local matrix vector multiplies are required for resiliency in a parallel setting.

The rest of the paper is organized as follows. Section 2 provides an overview of related work. Section 3 provides background on ABFT techniques for the PCG algorithm, and presents both the ABFT-Detection and ABFT-Correction approaches. Section 4 discusses the proposed ABFT method for the SpMxV. Section 5 is devoted to the abstract performance model. Section 6 reports numerical simulations comparing the performance of ABFT-Detection, ABFT-Correction and Online-Detection. Finally, we outline main conclusions and directions for future work in Section 7.

Related work

We classify related work along the following topics: silent errors in general, verification mechanisms for iterative methods, and ABFT techniques.

Silent errors

Considerable efforts have been directed at error-checking to reveal silent errors. Error detection is usually very costly. Hardware mechanisms, such as ECC memory, can detect and even correct a fraction of errors, but in practice they are complemented with software techniques. The simplest technique is triple modular redundancy and voting [START_REF] Lyons | The use of triple-modular redundancy to improve computer reliability[END_REF], which induces a costly verification. For high-performance scientific applications, process replication (each process is equipped with a replica, and messages are quadruplicated) is proposed in the RedMPI library [START_REF] Fiala | Detection and correction of silent data corruption for large-scale highperformance computing[END_REF]. Elliot et al. [START_REF] Elliott | Combining partial redundancy and checkpointing for HPC[END_REF] combine partial redundancy and checkpointing, and confirm the benefit of dual and triple redundancy. The drawback is that twice the number of processing resources is required (for dual redundancy). A comprehensive list of general-purpose techniques and references is provided by Lu et al. [START_REF] Lu | When is multi-version checkpointing needed[END_REF].

Application-specific information can be very useful to enable ad-hoc solutions, which dramatically decrease the cost of detection. Many techniques have been advocated. They include memory scrubbing [START_REF] Hwang | Cosmic rays don't strike twice: understanding the nature of dram errors and the implications for system design[END_REF] and ABFT techniques (see below).

As already stated, most papers assume on a selective reliability setting [START_REF] Hoemmen | Fault-tolerant iterative methods via selective reliability[END_REF][START_REF] Hoemmen | Fault-Tolerant Iterative Methods via Selective Reliability[END_REF][START_REF] Bridges | Fault-tolerant linear solvers via selective reliability[END_REF][START_REF] Sao | Self-stabilizing iterative solvers[END_REF]. It essentially means that one can choose to perform any operation in reliable or unreliable mode, assuming the former to be error-free but energy consuming and the latter to be error-prone but preferable from an energy consumption point of view.

Iterative methods

Iterative methods offer a wide range of ad-hoc approaches. For instance, instead of duplicating the computation, Benson et al. [START_REF] Benson | Silent error detection in numerical time-stepping schemes[END_REF] suggest coupling a higher-order with a lower-order scheme for PDEs. Their method only detects an error but does not correct it. Self-stabilizing corrections after error detection in the CG method are investigated by Sao and Vuduc [START_REF] Sao | Self-stabilizing iterative solvers[END_REF]. Heroux and Hoemmen [START_REF] Heroux | Fault-tolerant iterative methods via selective reliability[END_REF] design a fault-tolerant GMRES capable of converging despite silent errors. Bronevetsky and de Supinski [START_REF] Bronevetsky | Soft error vulnerability of iterative linear algebra methods[END_REF] provide a comparative study of detection costs for iterative methods.

As already mentioned, a nice instantiation of the checkpoint and verification mechanism that we study in this paper is provided by Chen [START_REF] Chen | Online-ABFT: An Online Algorithm Based Fault Tolerance Scheme for Soft Error Detection in Iterative Methods[END_REF], who deals with sparse iterative solvers. For PCG, the verification amounts to checking the orthogonality of two vectors and to recomputing and checking the residual (see Section 3 for further details).

As already mentioned, our abstract performance model is agnostic of the underlying error-detection technique and takes the cost of verification as an input parameter to the model.

ABFT

The very first idea of algorithm-based fault tolerance for linear algebra kernels is given by Huang and Abraham [START_REF] Huang | Algorithm-Based Fault Tolerance for Matrix Operations[END_REF]. They describe an algorithm capable of detecting and correcting a single silent error striking a matrix-matrix multiplication by means of row and column checksums. This germinal idea is then elaborated by Anfinson and Luk [START_REF] Anfinson | A Linear Algebraic Model of Algorithm-Based Fault Tolerance[END_REF], who propose a method to detect and correct up to two errors in a matrix representation using just four column checksums. Despite its theoretical merit, the idea presented in their paper is actually applicable only to relatively small matrices, and is hence out of our scope. Bosilca et al. [START_REF] Bosilca | Algorithm-based fault tolerance applied to high performance computing[END_REF] and Du et al. [START_REF] Du | Algorithm-based fault tolerance for dense matrix factorizations[END_REF] present two relatively recent survey.

The problem of algorithm-based fault-tolerance for sparse matrices is investigated by Shantharam et al. [START_REF] Shantharam | Fault Tolerant Preconditioned Conjugate Gradient for Sparse Linear System Solution[END_REF], who suggest a way to detect a single error in an SpMxV at the cost of a few additional dot products. Sloan et al. [START_REF] Sloan | Algorithmic Approaches to Low Overhead Fault Detection for Sparse Linear Algebra[END_REF] suggest that this approach can be relaxed using randomization schemes, and propose several checksumming techniques for sparse matrices. These techniques are less effective than the previous ones, not being able to protect the computation from faults striking the memory, but provide an interesting theoretical insight.

A preliminary version of this work appeared [START_REF] Fasi | Combining Backward and Forward Recovery to Cope with Silent Errors in Iterative Solvers[END_REF]. In the preliminary version, we dealt with CG only, without preconditioning. In practice, CG is most often used with preconditioning to accelerate the convergence. In the current work, we focus on preconditioned CG (PCG for short). We provide new results and experiments, together with several extensions on single error correction and multiple error detection.

CG-ABFT

We streamline our discussion on the CG method, however, the techniques that we describe are applicable to any iterative solver that use sparse matrix vector multiplies and vector operations. This list includes many of the nonstationary iterative solvers such as CGNE, BiCG, BiCGstab where sparse matrix transpose vector multiply operations also take place. In particular, we consider a PCG variant where the application of the preconditioner reduces to the computation of two SpMxV with triangular matrices [START_REF] Benzi | A sparse approximate inverse preconditioner for nonsymmetric linear systems[END_REF], which are a sparse factorization of an approximate inverse of the coefficient matrix. In fact, the model discussed in this paper can be profitably employed for any sparse inverse preconditioner that can be applied by means of one or more SpMxV.

We first provide a background on the CG method and overview both Chen's stability tests [START_REF] Chen | Online-ABFT: An Online Algorithm Based Fault Tolerance Scheme for Soft Error Detection in Iterative Methods[END_REF] and our ABFT protection schemes.

The code for the variant of the PCG method we use is shown in Algorithm 1. The main loop features three sparse matrix-vector multiply, two inner products (for p i q and r i+1

2 ), and three vector operations of the form axpy. Two of the sparse matrix-vector multiplications are carried out at line 10 where we compute z i+1 ← M Mr i+1 . Here the matrix M is the upper triangular factor of the approximate inverse M M, where M M well approximates A -1 [START_REF] Benzi | A sparse approximate inverse preconditioner for nonsymmetric linear systems[END_REF].

Chen's stability tests [START_REF] Chen | Online-ABFT: An Online Algorithm Based Fault Tolerance Scheme for Soft Error Detection in Iterative Methods[END_REF] amount to checking the orthogonality of vectors p i+1 and q, at the price of computing p i+1 q p i+1 q i , and to checking the residual at the price of an additional SpMxV operation Ax i -b. The dominant cost of these verifications is the additional SpMxV operation.

Algorithm 1 The PCG algorithm.

Require: A, M ∈ R n×n , b, x 0 ∈ R n , ε ∈ R Ensure: x ∈ R n : Ax -b ≤ ε 1: r 0 ← b -Ax 0 ; 2: z 0 ← M Mr 0 ; 3: p 0 ← z 0 ; 4: i ← 0; 5: while r i > ε ( A • r 0 + b ) do 6:
q ← Ap i ;

7:

α i ← r i 2 /p i q;

8:

x i+1 ← x i + α p i ; 9:
r i+1 ← r i -α q;

10:

z i+1 ← M Mr i+1 ; 11:
β ← r i+1 2 / r i 2 ;

12:

p i+1 ← z i+1 + β p i ;
13:

i ← i + 1; 14: end whilereturn x i ;

Our only modification to Chen's original approach is that we also save the sparse matrix A in addition to the current iteration vectors. This is needed when a silent error is detected: if this error comes for a corruption in data memory, we need to recover with a valid copy of the data matrix A. This holds for the three methods under study, Online-Detection, ABFT-Detection and ABFT-Correction, which have exactly the same checkpoint cost.

We now give an overview of our own protection and verification mechanisms. We use ABFT techniques to protect the SpMxV, its computations (hence the vector q), the matrix A, and the input vector p i . Since ABFT protection for vector operations is as costly as repeated computation, we use triple modular redundancy (TMR) for them for simplicity.

Although theoretically possible, constructing an ABFT mechanism to detect up to k errors is practically not feasible for k > 2. The same mechanism can be used to correct up to k/2 . Therefore, we focus on detecting up to two errors and correcting the error if there was only one. That is, we detect up to two errors in the computation q ← Ap i (two entries in q are faulty), or in p i , or in the sparse representation of the matrix A. With TMR, we assume that the errors in the computation are not overly frequent so that two out of three are correct (we assume errors do not strike the vector data here). Our fault-tolerant PCG versions thus have the following ingredients: ABFT to detect up to two errors in the SpMxV and correct the error, if there was only one; TMR for vector operations; and checkpoint and roll-back in case errors are not correctable.

We assume the selective reliability model in which all checksums and checksum related operations are non-faulty, also the tests for the orthogonality checks are non-faulty.

ABFT-SpMxV

Here, we discuss the proposed ABFT method for the SpMxV (combining ABFT with checkpointing is described later in Section 5.2). The proposed ABFT mechanisms are described for detecting single errors (Section 4.1), multiple errors (Section 4.2), and correcting a single error (Section 4.3).

Single error detection

The overhead of the standard single error correcting ABFT technique is too high for the sparse matrix-vector product case. Shantharam et al. [START_REF] Shantharam | Fault Tolerant Preconditioned Conjugate Gradient for Sparse Linear System Solution[END_REF] propose a cheaper ABFT-SpMxV algorithm that guarantees the detection of a single error striking either the computation or the memory representation of the two input operands (matrix and vector). Because their results depend on the sparse storage format adopted, throughout the paper we will assume that sparse matrices are stored in the compressed storage format by rows (CSR), that is by means of three distinct arrays, namely Colid ∈ N nnz(A) , Val ∈ R nnz(A) and Rowidx ∈ N n+1 ; these arrays correspond to the arrays JA, AA, and IA in the more standard notation [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]Sec. 3.4]). Here nnz(A) is the number of non-zero entries in A.

Shantharam et al. can protect y ← Ax, where A ∈ R n×n and x, y ∈ R n . To perform error detection, they rely on a column checksum vector c defined by

c j = n i=0 a i,j (1)
and an auxiliary copy x of the x vector. After having performed the actual SpMxV, to validate the result, it suffices to compute n i=1 y i , c x and c x , and to compare their values. It can be shown [START_REF] Shantharam | Fault Tolerant Preconditioned Conjugate Gradient for Sparse Linear System Solution[END_REF] that, in case of no errors, these three quantities carry the same value, whereas if a single error strikes either the memory or the computation, one of them must differ from the other two. Nevertheless, this method requires A to be strictly diagonally dominant, a condition that seems to restrict too much the practical applicability of their method. Shantharam et al. need this condition to ensure detection of errors striking an entry of x corresponding to a zero checksum column of A. We further analyze that case and show how to overcome the issue without imposing any restriction on A.

A nice way to characterize the problem is expressing it in geometrical terms. Consider the computation of a single entry of the checksum as

(w A) j = n i=1 w i a i,j = w A j ,
where w ∈ R n denotes the weight vector and A j the j-th column of A. Let us now interpret such an operation as the result of the scalar product •, • :

R n × R n → R defined by u, v → u v.
It is clear that a checksum entry is zero if and only if the corresponding column of the matrix is orthogonal to the weight vector. In (1), we have chosen w to be such that w i = 1 for 1 ≤ i ≤ n, in order to make the computation easier. Let us see now what happens without this restriction.

The problem reduces to finding a vector w ∈ R n that is not orthogonal to any vector out of a basis B = {b 1 , . . . , b n } of R n -the rows of the input matrix. Each of these n vectors is perpendicular to a hyperplane h i of R n , and w does not verify the condition

w, b i = 0 , (2)
for any i, if and only if it lies on h i . Since the Lebesgue measure in R n of an hyperplane of R n itself is zero, the union of these hyperplanes is measurable and has measure 0. Therefore, the probability that a vector w randomly picked in R n does not satisfy condition (2) for any i is zero. Nevertheless, there are many reasons to consider zero checksum columns. First of all, when working with finite precision, the number of elements in R n one can have is finite, and the probability of randomly picking a vector that is orthogonal to a given one could be larger than zero. Moreover, a coefficient matrix usually comes from the discretization of a physical problem, and the distribution of its columns cannot be considered as random. Finally, using a randomly chosen vector instead of (1, . . . , 1) increases the number of required floating point operations, causing a growth of both the execution time and the number of rounding errors (see Section 6). Therefore, we would like to keep w = (1, . . . , 1) as the vector of choice, in which case we need to protect SpMxV with matrices having zero column sums. There are many matrices with this property, for example the Laplacian matrices of graphs [START_REF] Chung | Spectral Graph Theory[END_REF]Ch. 1].

In Algorithm 2, we propose an ABFT SpMxV method that uses weighted checksums and does not require the matrix to be strictly diagonally dominant. The idea is to compute the checksum vector and then shift it by adding to all entries a constant value chosen so that all elements of the new vector are different from zero. We give the generalized result in Theorem 1.

Theorem 1 (Correctness of Algorithm 2). Let A ∈ R n×n be a square matrix, let x, y ∈ R n be the input and output vector respectively, and let x = x. Let us assume that the algorithm performs the computation

y ← A x , (3)
where A ∈ R n×n and x ∈ R n are the possibly faulty representations of A and x respectively, while y ∈ R n is the possibly erroneous result of the sparse matrixvector product. Let us also assume that the encoding scheme relies on 1. an auxiliary checksum vector

c = n i=1 a i,1 + m, . . . , n i=1 a i,n + m ,
where m is such that

c j = n i=1 a i,j + m = 0 , (4)
Algorithm 2 Shifting checksum algorithm.

Require: A ∈ R n×n , x ∈ R n Ensure: y ∈ R n such that y = Ax or the detection of a single error 1: x ← x; 2: [w, c, m, c r ] = computeChecksum(A); return SpMxV(A, x, x , w, c, m, c r ); 3: function computeChecksum(A) 4:
Generate w ∈ R n+1 ;

5:

w ← w 1:n ;

6: c ← w A; 7: if min(| c |) = 0 then; 8:
Find m that verifies (4); error ("Soft error is detected"); We consider these three possible cases in the following discussion.

d x ← c x ; 28: d r ← c r -s r ; 29: if d x = 0 ∧ d x = 0 ∧ d r = 0 then return y 1:n ; 30:
Case a. Let us assume, without loss of generality, that the error has struck at the pth position of y, which implies y i = y i for 1 ≤ i ≤ n with i = p and y p = y p +ε, where ε ∈ R \{0} represents the value of the error that has occurred. Summing up the elements of y gives n+1 i=1

y i = n i=1 n j=1 a i,j x j + m n j=1 x j + ε = n j=1 c j x j + ε = c x + ε , that violates condition (i).
Case b. A single error in the A matrix can strike one of the three vectors that constitute its sparse representation:

• a fault in Val that alters the value of an element a i,j implies an error in the computation of y i , which leads to the violation of the safety condition (i) because of (a),

• a variation in Colid can zero out an element in position a i,j shifting its value in position a i,j , leading again to an erroneous computation of y i ,

• a transient fault in Rowindex entails an incorrect value of s r and hence a violation of condition (iii).

Case c. Let us assume, without loss of generality, an error in position p of x. Hence we have that x i = x i for 1 ≤ i ≤ n with i = p and x p = x p + ε, for some ε ∈ R \ {0}. Noting that x = x , the sum of the elements of y gives n+1 i=1 Let us remark that computeChecksum in Algorithm 2 does not require the input vector x of SpMxV as an argument. Therefore, in the common scenario of many SpMxV with the same matrix, it is enough to invoke it once to protect several matrix-vector multiplications. This observation will be crucial when discussing the performance of these checksumming techniques.

y i = n i=1 n j=1 a i,j x j + m n j=1 x j = n i=1 n j=1 a i,j x j + m n j=1 x j + ε n i=1 a i,p + εm = n j=1 c j x j + ε n i=1 a i,p + m = c x + ε n i=1 a i,p + m ,
Shifting the sum checksum vector by an amount is probably the simplest deterministic approach to relax the strictly diagonal dominance hypothesis, but it is not the only one. An alternative solution is described in Algorithm 3, which basically exploits the distributive property of matrix multiplication over matrix addition. The idea is to split the original matrix A into two matrices of the same size, A and A, such that no column of either matrix has a zero checksum. Two standard ABFT multiplications, namely y ← Ax and y ← A x, are then performed. If no error occurs neither in the first nor in the second computation, the sum of y and y is computed in reliable mode and then returned. Let us note that, as we expect the number of non-zeros of A to be much smaller than n, we store sparsely both the checksum vector of A and the y vector.

We do not write down an extended proof of the correctness of this algorithm, and limit ourselves to a short sketch. We consider the same three cases as in the proof of Theorem 1, without introducing any new idea. An error in the computation of y or y can be detected using the dot product between the corresponding column checksum and the x error. An error in A can be detected by either c r or an erroneous entry in y or y, as the matrix loop structure of the sparse multiplication algorithm has not been changed. Finally, an error in the pth component of x would make the sum of the entries of y and y differ from c x and c x , respectively. The evaluation of the performance of the two algorithms, though straightforward from the point of view of the computational cost, has to be carefully assessed in order to devise a valid and practical trade-off between Algorithm 2 and Algorithm 3.

In both cases computeChecksum introduces an overhead of O (nnz(A)), but the shift version should in general be faster containing less assignments than its counterpart, and this changes the constant factor hidden by the asymptotic notation. Nevertheless, as we are interested in performing many SpMxV with the same matrix, this pre-processing overhead becomes negligible.

The function SpMxV has to be invoked once for each multiplication, and hence more care is needed. Copying x and initializing y both require n operations, and the multiplication is performed in time O (nnz(A)), but the split version pays an n more to read the values of the sparse vector b. The cost of the verification step depends instead on the number of zeroes of the original checksum vector, that is also the number of non-zero elements of the sparse vector c. Let us call this quantity n . Then the overhead is 4n for the shifting and 3n + 3n for the splitting, that requires also the sum of two sparse vectors to return the result. Hence, as summarized in Table 1, the two methods bring different overheads into the computation. Comparing them, it is immediate to see that the shifting method is cheaper when

n > n 5 ,
while it has more operations to do when the opposite inequality holds. For the equality case, we can just choose to use the first method, because of the cheaper preprocessing phase. In view of this observation, it is possible to devise a simple algorithm that exploits this trade-off to achieve better performance. It suffices to compute the checksum vector of the input matrix, count the number of non-zeros and choose which detection method to use accordingly.

We also note that by splitting the matrix A into say pieces and checksumming each piece separately, we can possibly protect A from up to errors, by Algorithm 3 Splitting checksum algorithm. c, m ← 0;

4:

for i ← 1 to nnz(A) do 5:
ind ← Colid i ;

6:

c ind ← c ind + Val i ; 7: m ind ← i; 8:
end for

9:
k ← 0;

10:

for i ← 1 to n do 11:

if c i = 0 ∧ m i = 0 then 12:
b m i ← true;

13:

c i ← Val m i ;
14:

c i ← c i -c i ; 15: end if 16:
end for

y i ← y i + Val j • x ind ; 31:
y i ← y i + Val j •
d x ← c x -c y ; d x ← c x -c y ;
38: error ("Soft error is detected");

d x ← c x -c y ; d x ← c x -c y ; 39: d r ← c r -s r ; 40: if d x = 0 ∧ d x = 0 ∧ d x = 0 ∧ d x = 0 ∧ d r =

43:

end if 44: end function protecting each piece against a single error (obviously the multiple errors should hit different pieces).

Multiple error detection

With some effort, the shifting idea in Algorithm 2 can be extended to detect errors striking a single SpMxV. Let us consider the problem of detecting up to k errors in the computation of y ← Ax introducing an overhead of O (kn). Let k weight vectors w (1) , . . . , w (k) ∈ R n be such that any sub-matrix of W = w (1) w (2) . . . w (k) has full rank. For k = 2 and k = 3, we give a possible choice for the W matrix explicitly. For larger values of k, we suggest the use of rectangular random matrices, which besides having full rank sub-matrices are also well conditioned with high probability [START_REF] Edelman | Eigenvalues and condition numbers of random matrices[END_REF][START_REF] Chen | Condition numbers of Gaussian random matrices[END_REF].

To build our ABFT scheme let us note that, if no error occurs, for each weight vector w ( ) it holds that

w ( ) A = n i=1 w ( ) i a i,1 , . . . , n i=1 w ( ) i a i,n ,
and hence that

w ( ) Ax = n i=1 w ( ) i a i,1 x 1 + • • • + n i=1 w ( ) i a i,n x n = n i=1 n j=1 w ( ) i a i,j x j .
Similarly, the sum of the entries of y weighted with the same w ( )

is n i=1 w ( ) i y i = w ( ) 1 y 1 + • • • + w ( ) n y n = w ( ) 1 n j=1 a 1,j x j + • • • + w ( ) n n j=1 a n,j x j = n i=1 n j=1 w ( ) i a i,j x j ,
and we can conclude that n i=1 w ( )

i y i = w ( ) A x ,
for any w ( ) with 1 ≤ ≤ k.

To convince ourself that with these checksums it is actually possible to detect up to k errors, let us suppose that k errors, with k ≤ k, occur in positions p 1 , . . . , p k , and let us denote by y the faulty vector where y p i = y p i + ε p i for ε p i ∈ R \ {0} and 1 ≤ i ≤ k and y i = y i otherwise. Then for each weight vector we have

n i=1 w ( ) i y i - n i=1 w ( ) i y i = k j=1 w ( ) p j ε p j .
Said otherwise, the occurrence of the k errors is not detected if and only if, for 1 ≤ ≤ k, all the ε p i respect k j=1 w ( )

p j ε p j = 0 . (5)
We claim that there cannot exist a vector (ε p 1 , . . . , ε p k ) ∈ R k \ {0} such that all the conditions in ( 5) are simultaneously verified. These conditions can be expressed in a more compact way as a linear system

    w (1) p 1 • • • w (1) p k . . . . . . . . . w (k) p 1 • • • w (k) p k         ε p 1 . . . ε p k     =    0 . . . 0    .
Denoting by W * the coefficient matrix of this system, it is clear that the errors cannot be detected if only if (ε p 1 , . . . ε p k ) ∈ ker(W * ) \ {0}. Because of the properties of W, we have that rk(W * ) = k. Moreover, it is clear that the rank of the augmented matrix

    w (1) p 1 • • • w (1) p k . . . . . . . . . w (k) p 1 • • • w (k) p k 0 . . . 0    
is k as well. Hence, by means of the Rouché-Capelli theorem, the solution of the system is unique and the null space of W * is trivial. Therefore, this construction can detect the occurrence of k errors during the computation of y by comparing the values of the weighted sums y w ( ) with the result of the dot product (w ( ) A)x, for 1 ≤ ≤ k. However, to get a true extension of the algorithm described in the previous section, we also need to make it able to detect errors that strike the sparse representation of A and that of x. The first case is simple, as the k errors can strike the Val or Colid arrays, leading to at most k errors in y, or in Rowindex , where they can be caught using k weighted checksums of the Rowindex vector.

Detection in x is much trickier, since neither the algorithm just described nor a direct generalization of Algorithm 2 can manage this case. Nevertheless, a proper extension of the shifting technique is still possible. Let us note that there exists a matrix M ∈ R k×n such that

W A + M = W .
The elements of such an M can be easily computed, once that the checksum rows are known. Let x ∈ R n be the faulty vector, defined by

x i =    x i + ε p i , 1 ≤ i ≤ k , x i otherwise.
for some k ≤ k, and let us define y = A x. Now, let us consider a checksum vector x ∈ R n such that x = x and let assume that it cannot be modified by a transient error. For 1 ≤ ≤ k, it holds that n i=1 w ( )

i y i + n j=1 m ,j x j = n i=1 n j=1 w ( ) i a i,j x j + k i=1 ε p i   n j=1 w ( ) j a j,p i   + n j=1 m ,j x j + k i=1 ε p i m l,p i = n j=1 n i=1 w ( ) i a i,j x j + n j=1 m ,j x j + k i=1 ε p i   n j=1 w ( ) j a j,p i + m l,p i   = n j=1 n i=1 w ( ) i a i,j x j + n j=1 m ,j x j + k i=1 ε p i w ( ) p i = n j=1 n i=1 w ( ) i a i,j + m ,j x j + k i=1 ε p i w ( ) p i = w ( ) x + k i=1 ε p i w ( ) p i = w ( ) x + k i=1 ε p i w ( ) p i .
Therefore, an error is not detected if and only if the linear system

    w (1) p 1 • • • w (1) p k . . . . . . . . . w (k) p 1 • • • w (k) p k         ε p 1 . . . ε p k     =    0 . . . 0   
has a non-trivial solution. But we have already seen that such a situation can never happen, and we can thus conclude that our method, whose pseudocode we give in Algorithm 4, can also detect up to k errors occurring in x. Therefore, we have proven the following theorem.

Theorem 2 (Correctness of Algorithm 4). Let us consider the same notation as in Theorem 1. Let W ∈ R n+1×n be a matrix such that any square submatrix is full rank, and let us denote by W ∈ R n×n the matrix of its first n rows. Let us assume an encoding scheme that relies on Then, up to k errors striking the computation of y or the memory locations that store A or x, cause one of the following conditions to fail:

i. W y = C x, ii. W (x -y),
iii. s Rowindex = c Rowindex .

Let us note that we have just shown that our algorithm can detect up to k errors striking only A, or only x or only the computation. Nevertheless, this result holds even when the errors are distributed among the possible cases, as long as at most k errors rely on the same checkpoint.

It is clear that the execution time of the algorithm depends on both nnz(A) and k. For the computeChecksum function, the cost is, assuming that the weight matrix W is already known, O (k nnz(A)) for the computation of C, and O (kn) for the computation of M and c Rowindex . Hence the number of performed operations is O (k nnz(A)). The overhead added to the SpMxV depends instead on the computation of four checksum matrices that lead to a number of operations that grows asymptotically as kn.

Single error correction

We now discuss single error correction, using Algorithm 4 as a reference. We describe how a single error striking either memory or computation can be not only detected but also corrected at line 27. We use only two checksum vectors, that is, we describe correction of single errors assuming that two errors cannot strike the same SpMxV. By the end of the section, we will generalize this approach and discuss how single error correction and double error detection can be performed concurrently by exploiting three linearly independent checksum vectors. Whenever a single error is detected, regardless of its location (computation or memory), it is corrected by means of a succession of various steps. When one or more errors are detected, the correction mechanism tries to determine their multiplicity and, in case of a single error, what memory locations have been corrupted or what computation has been miscarried. Errors are then corrected using the values of the checksums and, if need be, partial recomputations of the result are performed.

As we did for multiple error detection, we require that any 2 × 2 submatrix of W ∈ R n×2 has full rank. The simplest example of weight matrix having this property is probably

W =        1 1 1 2 1 3 . . . . . . 1 n        .
To detect errors striking Rowidx , we compute the ratio ρ of the second component of d r to the first one, and check whether its distance from an integer is smaller than a certain threshold parameter ε. If this distance is smaller, the algorithm concludes that the σth element, where σ = Round(ρ) is the nearest integer to ρ, of Rowidx is faulty, performs the correction by subtracting the first component of d r from Rowidx σ , and recomputes y σ and y σ-1 , if the error in Algorithm 4 Shifting checksum algorithm for k errors detection.

Require: A ∈ R n×n , x ∈ R n Ensure: y ∈ R n such that y = Ax or the detection of up to k errors 1: x ← x 2: [W, C, M, c] = computeChecksums(A, k); return SpMxV(A, x, x , W, C, M, k, c); 3: function computeChecksums(A, k) 4:
Generate W = w (1) . . . w (k) ∈ R n+1×n ;

5:

W ← W 1:n, * ∈ R n×k 6: C ← W A; 7:
M ← W -C; s ← [0, . . . , 0]; • if z C > 2, then errors can be detected but not corrected, and an error is emitted.

d x ← W y -C x; 24: d x ← W (x -y) -M x; 25: d r ← c -s; 26: if d x = 0 ∧ d x = 0 ∧ d r = 0 then return y; 27:
To correct errors striking x, the algorithm computes ρ, that is the ratio of the second component of d x to the first one, and checks that the distance between d and the nearest integer σ is smaller than ε. Provided that this condition is verified, the algorithm computes the value of the error τ = n i=1 x i -cx σ and corrects x σ = x σ -τ . The result is updated by subtracting from y the vector y τ = Ax τ , where x τ ∈ R n×n is such that x τ σ = τ and x τ i = 0 otherwise. Let us now investigate how detection and correction can be combined and let us give some details about the implementation of ABFT-Correction as defined in Section 3. Indeed, note that double errors could be shadowed when using Algorithm 2, although the probability of such an event is negligible.

Let us restrict ourselves to an easy case, without considering errors in x. As usual, we compute the column checksums matrix C = (W A) , and then compare the two entries of C x ∈ R 2 with the weighted sums

y c 1 = n i=1 y i and y c 2 = n i=1 i y i
where y is the possibly faulty vector computed by the algorithm. It is clear that if no error occurs, the computation verifies the condition δ = y cc = 0. Furthermore, if exactly one error occurs, we have δ 1 , δ 2 = 0 and δ 2 δ 1 ∈ N, and if two strike the vectors protected by the checksum c, the algorithm is able to detect them by verifying that δ = 0.

At this point it is natural to ask whether this information is enough to build a working algorithm or some border cases can bias its behavior. In particular, when δ 2 δ 1 = p ∈ N, it is not clear how to discern between single and double errors. Let ε 1 , ε 2 ∈ R \ {0} be the value of two errors occurring at position p 1 and p 2 respectively, and let y ∈ R n be such that

y i =    y i , 1 ≤ i ≤ n, i = p 1 , p 2 y i + ε 1 , i = p 1 y i + ε 2 , i = p 2 .
Then the conditions

δ 1 = ε 1 + ε 2 , (6)
δ 2 = p 1 ε 1 + p 2 ε 2 , (7)
hold. Therefore, if ε 1 and ε 2 are such that

p (ε 1 + ε 2 ) = p 1 ε 1 + p 2 ε 2 , (8)
it is not possible to distinguish these two errors from a single error of value ε 1 + ε 2 occurring in position p. This issue can be solved by introducing a new set of weights and hence a new row of column checksums. Let us consider a weight matrix W ∈ R n×3 that includes a quadratic weight vector

W =        1 1 1 1 2 4 1 3 9 . . . . . . . . . 1 n n 2       
, and the tridimensional vector

δ = W A x -W y ,
whose components can be expressed as

δ 1 = ε 1 + ε 2 δ 2 = p 1 ε 1 + p 2 ε 2 δ 3 = p 2 1 ε 1 + p 2 2 ε 2 .
To be confused with a single error in position p, ε 1 and ε 2 have to be such that

p (ε 1 + ε 2 ) = p 1 ε 1 + p 2 ε 2 and p 2 (ε 1 + ε 2 ) = p 2 1 ε 1 + p 2 2
ε 2 hold simultaneously for some p ∈ N. In other possible values of the errors are the solution of the linear system

(p -p 1 ) (p -p 2 ) (p 2 -p 2 1 ) (p 2 -p 2 2 ) ε 1 ε 2 = 0 0 .
It is easy to see that the determinant of the coefficient matrix is

(p -p 1 ) (p -p 2 ) (p 2 -p 1 ) ,
which always differs from zero, as long as p, p 1 and p s differ pairwise. Thus, the matrix is invertible, and the solution space of the linear system is the trivial kernel (ε 1 , ε 2 ) = (0, 0). Thus using W as weight matrix guarantees that it is always possible to distinguish a single error from double errors.

Performance model

In Section 5.1, we introduce the general performance model. Then in Section 5.2 we instantiate it for the three methods that we are considering, namely Online-Detection, ABFT-Detection and ABFT-Correction.

General approach

We introduce an abstract performance model to compute the best combination of checkpoints and verifications for iterative methods. We execute T time-units of work followed by a verification, which we call a chunk, and we repeat this scheme s times, i.e., we compute s chunks, before taking a checkpoint. We say that the s chunks constitute a frame. The whole execution is then partitioned into frames. We assume that checkpoint, recovery and verification are error-free operations. Let T cp , T rec and T verif be the respective cost of these operations. Finally, assume an exponential distribution of errors and let q be the probability of successful execution for each chunk: q = e -λT , where λ is the fault rate.

The goal of this section is to compute the expected time E (s, T ) needed to execute a frame composed of s chunks of size T . We derive the best value of s as a function of T and of the resilience parameters T cp , T rec , T verif , and q, the success probability of a chunk. Each frame is preceded by a checkpoint, except maybe the first one (for which we recover by reading initial data again). Following earlier work [START_REF] Bougeret | Checkpointing strategies for parallel jobs[END_REF], we derive the following recursive equation to compute the expected completion time of a single frame:

E (s, T ) = q s (s(T + T verif )) + T cp ) + (1 -q s ) (E (T lost ) + T rec + E (s, T )) . (9)
Indeed, the execution is successful if all chunks are successful, which happens with probability q s , and in this case the execution time simply is the sum of the execution times of each chunk plus the final checkpoint. Otherwise, with probability 1 -q s , we have an error, which we detect after some time E (T lost ), and that forces us to recover (in time T rec ) and restart the frame anew, hence in time E (s, T ). The difficult part is to compute E (T lost ).

For 1 ≤ i ≤ s, let f i be the following conditional probability:

f i = P(error strikes at chunk i|there is an error in the frame) . (

Given the success probability q of a chunk, we obtain that

f i = q i-1 (1 -q) 1 -q s .
Indeed, the first i -1 chunks were successful (probability q i-1 ), the ith one had an error (probability 1 -q), and we condition by the probability of an error within the frame, namely 1 -q s . With probability f i , we detect the error at the end of the ith chunk, and we have lost the time spent executing the first i chunks. We derive that

E (T lost ) = s i=1 f i (i(T + T verif )) .
We have

s i=1 f i = (1-q)h(q) 1-q s where h(q) = 1 + 2q + 3q 2 + • • • + sq s-1 . If m(q) = q +q 2 +• • •+q s = 1-q s+1
1-q -1, we get by differentiation that m (q) = h(q), hence h(q) = -(s+1)q s 1-q + 1-q s+1

(1-q) 2 and finally

E (T lost ) = (T + T verif ) sq s+1 -(s + 1)q s + 1 (1 -q s )(1 -q) .
Plugging the expression of E (T lost ) back into (9), we obtain

E (s, T ) = s(T + T verif )) + T cp + (q -s -1)T rec + T sq s+1 -(s + 1)q s + 1 q s (1 -q) ,
which simplifies into

E (s, T ) = T cp + (q -s -1)T rec + (T + T verif ) 1 -q s q s (1 -q) .
We have to determine the value of s that minimizes the overhead of a frame:

s = argmin s≥1 E (s, T ) sT . ( 11
)
The minimization is complicated and should be conducted numerically (because T , the size of a chunk, is still unknown). Luckily, a dynamic programming algorithm to compute the optimal value of T and s is available [START_REF] Benoit | Assessing general-purpose algorithms to cope with fail-stop and silent errors[END_REF].

Instantiation to PCG

For each of the three methods, Online-Detection, ABFT-Detection and ABFT-Correction, we instantiate the previous model and discuss how to solve [START_REF] Huang | Algorithm-Based Fault Tolerance for Matrix Operations[END_REF].

Online-Detection

For Chen's method [START_REF] Chen | Online-ABFT: An Online Algorithm Based Fault Tolerance Scheme for Soft Error Detection in Iterative Methods[END_REF], we have chunks of d iterations, hence T = dT iter , where T iter is the raw cost of a PCG iteration without any resilience method. The verification time T verif is the cost of the orthogonality check operations performed as described in Section 3. As for silent errors, the application is protected from arithmetic errors in the ALU, as in Chen's original method, but also for corruption in data memory (because we also checkpoint the matrix A). Let λ a be the rate of arithmetic errors, and λ m be the rate of memory errors. For the latter, we have λ m = M λ word if the data memory consists of M words, each susceptible to be corrupted with rate λ word . Altogether, since the two error sources are independent, they have a cumulative rate of λ = λ a + λ m , and the success probability for a chunk is q = e -λT .

Plugging these values in [START_REF] Huang | Algorithm-Based Fault Tolerance for Matrix Operations[END_REF] gives an optimization formula very similar to that of Chen [7, Sec. 5.2], the only difference being that we assume that the verification is error-free, which is needed for the correctness of the approach.

ABFT-Detection

When using ABFT techniques, we detect possible errors every iteration, so a chunk is a single iteration, and T = T iter . For ABFT-Detection, T verif is the overhead due to the checksums and redundant operations to detect a single error in the method.

ABFT-Detection can protect the application from the same silent errors as Online-Detection, and just as before the success probability for a chunk (a single iteration here) is q = e -λT .

ABFT-Correction

In addition to detection, we now correct single errors at every chunk. Just as for ABFT-Detection, a chunk is a single iteration, and T = T iter , but T verif corresponds to a larger overhead, mainly due to the extra checksums needed to detect two errors and correct a single one.

The main difference lies in the error rate. An iteration with ABFT-Correction is successful if zero or one error has struck during that iteration, so that the success probability is much higher than for Online-Detection and ABFT-Detection. We compute that value of the success probability as follows. We have a Poisson process of rate λ, where λ = λ a + λ m as for Online-Detection and ABFT-Detection. The probability of exactly k errors in time T is [START_REF] Mitzenmacher | Probability and Computing: Randomized Algorithms and Probabilistic Analysis[END_REF], hence the probability of no error is e -λT and the probability of exactly one error is λT e -λT , so that q = e -λT + λT e -λT .

(λT ) k k! e -λT

Experiments

Setup

There are two different sources of advantages in combining ABFT and checkpointing. First, the error detection capability lets us perform a cheap validation of the partial result of each PCG step, recovering as soon as an error strikes. Second, being able to correct single errors makes each step more resilient and increases the expected number of consecutive valid iterations. We say an iteration is valid if it is non-faulty, or if it suffers from a single error that is corrected via ABFT.

For our experiments, we use a set of positive definite matrices from the UFL Sparse Matrix Collection [START_REF] Davis | The University of Florida Sparse Matrix Collection[END_REF], with size between 17456 and 74752 and density lower than 10 -2 . These matrices are listed in Table 2. Among these matrices, "Andrews/Andrews" and "GHS psdef/jnlbrng1" are diagonally dominant, others are not. We perform the experiments under Matlab and use the factored approximate inverse preconditioners [START_REF] Benzi | A sparse approximate inverse preconditioner for nonsymmetric linear systems[END_REF][START_REF] Benzi | Stabilized and block approximate inverse preconditioners for problems in solid and structural mechanics[END_REF] in the PCG. The application of these preconditioners requires two SpMxV, which are protected against error using the methods proposed in Section 4 (in all methods Online-Detection, ABFT-Detection, and ABFT-Correction).

At each iteration of PCG, faults are injected during vector and matrix-vector operations but, since we are assuming selective reliability, all the checksums and checksum operations are considered non-faulty. Faults are modeled as bit flips occurring independently at each step, under an exponential distribution of parameter λ, as detailed in Section 5.2. These bit flips can strike either the matrix (the elements of Val , Colid and Rowidx ), or any entry of the PCG vectors r i , z i , q, p i or x i . We chose not to inject errors during the computation explicitly, as they are just a special case of error we are considering. Moreover, to simplify the injection mechanism, T iter is normalized to be one, meaning that each memory location or operation is given the chance to fail just once per iteration [START_REF] Sao | Self-stabilizing iterative solvers[END_REF]. Finally, to get data that are homogeneous among the test matrices, the fault rate λ is chosen to be inversely proportional to M (memory size) with a proportionality constant α ∈ (0, 1); this makes sense as larger the memory used by an application, larger is the chance to have an error. It follows that the expected number of PCG iterations between two distinct fault occurrences does not depend either on the size or on the sparsity ratio of the matrix.

We compare the performance of three algorithms, namely Online-Detection, ABFT-Detection (single detection and rolling back as soon as an error is detected), and ABFT-Correction (correcting single errors during a given step and rolling back only if two errors strike a single operation). We instantiate them by limiting the maximum number of PCG steps to 50 (20 for #924, whose convergence is sublinear) and setting the tolerance parameter at line 5 of Algorithm 1 to 10 -14 . The number of iterations for a non-faulty execution and the achieved accuracy are detailed in Table 2.

Implementing the null checks in Algorithm 2, Algorithm 3 and Algorithm 4 poses a challenge. The comparison d r = 0 is between two integers, and can be correctly evaluated by any programming language using the equality check.

(c A) x) -f l (c (Ax)) | ≤ 2 γ 2n |c | | A | | x | . (12)
We refer the reader to the technical report for the proof [44, Theorem 2]. Let us note that if all of the entries of c are positive, as it is often the case in our setting, the absolute value of c in ( 12) can be safely replaced with c itself. It is also clear that these bounds are not computable, since c | A | | x | is not, in general, a floating point number. This problem can be alleviated by overestimating the bound by means of matrix and vector norms.

Since we are interested in actually computing the bound at runtime, we consider a weaker bound. Recalling that [START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF]Sec. B.7]

A 1 = max 1≤j≤n n i=1 |a i,j | . ( 13
)
we can upper bound the right hand side in so that Though the right hand side of ( 14) is not exactly computable in floating point arithmetic, it requires an amount of operations dramatically smaller than [START_REF] Du | Algorithm-based fault tolerance for dense matrix factorizations[END_REF]; just a few sums for the norm of A. As this norm is usually computed using the identity in [START_REF] Yao | Detection of soft errors in LU decomposition with partial pivoting using algorithm-based fault tolerance[END_REF], any kind of summation yields a relative error of at most n u [43, Sec. 4.6], where n is the maximum number of nonzeros in a column of A, and u is the machine epsilon. Since we are dealing with sparse matrices, we expect n to be very small, and hence the computation of the norm to be accurate. Moreover, since the right hand side in [START_REF] Hakkarinen | Fail-stop failure algorithm-based fault tolerance for cholesky decomposition, Parallel and Distributed Systems[END_REF] does not depend on x, it can be computed just once for a given matrix and weight vector. Using [START_REF] Hakkarinen | Fail-stop failure algorithm-based fault tolerance for cholesky decomposition, Parallel and Distributed Systems[END_REF] as tolerance parameter guarantees no false positives (a computation without any error that is considered as faulty), but allows false negatives (an iteration during which an error occurs without being detected) when the perturbations of the result are small. Nonetheless, this solution works almost perfectly in practice, meaning that though the convergence rate can be slowed down, the algorithms still converges towards the "correct" answer. Though such an outcome could be surprising at first, Elliott et al. [START_REF] Elliott | Quantifying the impact of single bit flips on floating point arithmetic[END_REF][START_REF] Stoyanov | Quantifying the impact of single bit flips on floating point arithmetic[END_REF] showed that bit flips that strike the less significant digits of the floating point representation of vector elements during a dot product create small perturbations of the results, and that the magnitude of this perturbation gets smaller as the size of the vectors increases. Hence, we expect errors that are not detected by our tolerance threshold to be too small to impact the solution of the linear solver.

|f l ((c A) x) -f l (c (Ax)) | ≤ 2 γ 2n n c ∞ A 1 x ∞ . (14)

Simulations

To validate the model, we perform the simulation whose results are illustrated in Table 3. For each matrix, we set λ = 1 16 M and consider the average execution time of 100 repetitions of both ABFT-Detection (columns 5-8) and ABFT-Correction (columns 6-9). In the table we record the checkpointing interval s * i which achieves the shortest execution time E t (s * 1 ), and the checkpointing interval s i which is the best stepsize according to our method, along with its execution time E t ( s i ). Finally, we evaluate the performance of our guess by means of the quantity

l i = E t ( s i ) -E t (s * i ) E t (s * i )
• 100 , that expresses the loss, in terms of execution time, of executing with the checkpointing interval given by our model with respect to the best possible choice.

From the table, we clearly see that the values of s i and s * i are close, since the time loss reaches just above 5% for l 1 and just below 15% for l 2 . This sometimes poor result depends just on the small number of repetitions we are considering, that leads to the presence of outliers, lucky runs in which a small number of errors occur and the computation is carried on in a much quicker way. Similar results hold for other values of λ.

We also compare the execution time of the three algorithms to empirically asses how much their relative performance depend on the fault rate. The results on our test matrices are shown in Fig. 1, where the y-axis is the execution time (in seconds), and the x-axis is the normalized mean time between failure (the reciprocal of α). Here, the larger x = 1 α , the smaller the corresponding value of λ = α M , hence the smaller the expected number of errors. For each value of λ, we draw the average execution time of 50 runs of the three algorithms, using the best checkpointing interval predicted in Section 5.1 for ABFT-Detection and ABFT-Correction, and by Chen [START_REF] Chen | Online-ABFT: An Online Algorithm Based Fault Tolerance Scheme for Soft Error Detection in Iterative Methods[END_REF]Eq. 10] for Online-Detection. In terms of execution time, Chen's method is comparable to ours for middle to high fault rates, since it clearly outperforms ABFT-Detection in five out of nine cases, being slightly faster than ABFT-Correction for lower fault rates.

Intuitively, this behavior is not surprising. When λ is large, many errors occur but, since α is between zero and one, we always have, in expectation, less than one error per iteration. Thus ABFT-Correction requires fewer checkpoints than ABFT-Detection and almost no rollback, and this compensates for the slightly longer execution time of a single step. When the fault rate is very low, instead, the algorithms perform almost the same number of iterations, but ABFT-Correction takes slightly longer due to the additional dot products at each step.

Altogether, the results show that ABFT-Correction outperforms both Online-Detection and ABFT-Detection for a wide range of fault rates, thereby demonstrating that combining checkpointing with ABFT correcting techniques is more efficient than pure checkpointing for most practical situations.

Conclusion

We consider the problem of silent errors in iterative linear systems solvers. At first, we focus our attention on ABFT methods for SpMxV, developing algorithms able to detect and correct errors in both memory and computation using various checksumming techniques. Then, we combine ABFT with replication, in order to develop a resilient PCG kernel that can protect axpy's and dot products as well. We also discuss how to take numerical issues into account when dealing with actual implementations. These methods are a worthy choice for a selective reliability model, since most of the operations can be performed in unreliable mode, whereas only checksum computations need to be performed reliably.

In addition, we examine checkpointing techniques as a tool to improve the resilience of our ABFT PCG and develop a model to trade-off the checkpointing interval so to achieve the shortest execution time in expectation. We implement two of the possible combinations, namely an algorithm that relies on roll back as soon as an error is detected, and one that is able to correct a single error and recovers from a checkpoint just when two errors strike. We validate the model by means of simulations and finally compare our algorithms with Chen's approach, empirically showing that ABFT overhead is usually smaller than Chen's verification cost.

We expect this combined approach to be interesting for other variants of the preconditioned conjugate gradient algorithm [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. Triangular preconditioners seem to be particularly attracting, in that it looks possible to treat them by adapting the techniques described in this paper (Shantharam et al. [START_REF] Shantharam | Fault Tolerant Preconditioned Conjugate Gradient for Sparse Linear System Solution[END_REF] addressed the triangular case).

We discussed the adequacy of the proposed techniques in the single core and distributed memory settings, but did not address shared memory parallel systems with many-cores. In these latter systems, the task-level parallelism necessitates application of accuracy theorems (in the spirit of Theorem 3) to avoid false positives. We leave the resilience of task-parallel executions as a future work.

  r ← w Rowindex ; return w, c, m, c r ; 12: end function 13: function SpMxV(A, x, x , w, c, m, c r ) 14: w ← w 1:n ; 15: s r ← 0; 16: for i ← 1 to n do 17: y i ← 0; 18: s r ← s r + Rowindex i ; 19: for j ← Rowindex i to Rowindex i+1 -1 do 20: ind ← Colid j ; 21:y i ← y i + Val j • x ind ; n+1 ← m w x ;25: c y ← w y; 26: d x ← c x; 27:

  end function for 1 ≤ j ≤ n, 2. an auxiliary checksum y n+1 = m n i=i x i , 3. an auxiliary counter s r initialized to 0 and updated at runtime by adding the value of the hit element each time the Rowindex array is accessed (line 20 of Algorithm 2), 4. an auxiliary checksum c r = n i=1 Rowindex i ∈ N. Then, a single error in the computation of the SpMxV causes one of the following conditions to fail: i. c x = n+1 i=1 y i , ii. c x = n+1 i=1 y i , iii. s r = c r . Proof. There are three possible error sources when computing Equation (3), namely: a. a faulty arithmetic operation during the computation of y, b. a bit flip in the sparse representation of A, c. a bit flip in an element of of x.

  that violates (ii) since n i=1 a i,p + m = 0 by definition of m.

2 :

 2 Require: A ∈ R n×n , x ∈ R n Ensure: y ∈ R n such that y = Ax or the detection of a single error 1: [c, k, c r ] = computeChecksum(A); return SpMxV(A, x, c, c, b, c r ); function computeChecksum(A) 3:

17 :

 17 c r ← n i=1 Rowindex i ; return c, c, b, c r ; 18: end function 19: function SpMxV(A, x, c, c, b, c r ) for i ← 1 to n do 26: s r ← s r + Rowindex i ; 27: for j ← Rowindex i to Rowindex i+1 -1 do 28: ind ← Colid j ; 29: if b j then 30:

  0 then return y + y;

1 .

 1 an auxiliary checksum matrix C = (W A) , 2. an auxiliary checksum matrix M = W -C, 3. a vector of auxiliary counters s Rowindex initialized to the null vector and updated at runtime as in lines 16 -17 of Algorithm 4), 4. an auxiliary checksum vector c Rowindex = W Rowindex .

8 :c

 8 Rowindex ← W Rowindex ; return W, C, M, c Rowindex ; 9: end function 10: function SpMxV(A, x, x , W, C, M, k, c) 11: W ← W 1:n, * ∈ R n×k 12:

  end function Rowindex σ is a decrement; or y σ+1 if it was an increment. Otherwise, it just emits an error.The correction of errors striking Val , Colid and the computation of y are performed together. Let now ρ be the ratio of the second component of d x to the first one. If ρ is near enough to an integer σ, the algorithm computes the checksum C = W A and considers the number z C of non-zero columns of the difference matrix C = |C -C |. At this stage, three cases are possible:• If z C = 0, then the error is in the computation of y σ , and can be corrected by simply recomputing this value.• If z C = 1,then the error has struck an element of Val . Let us call f the index of the non-zero column of C. The algorithm finds the element of Val corresponding to the entry at row σ and column f of A and corrects it by using the column checksums much like as described for Rowidx . Afterwards, y d is recomputed to fix the result. • If z C = 2, then the error concerns an element of Colid . Let us call f 1 and f 2 the index of the two non-zero columns and m 1 , m 2 the first and last elements of Colid corresponding to non-zeros in row σ. It is clear that there exists exactly one index m * between m 1 and m 2 such that either Colid m * = f 1 or Colid m * = f 2 . To correct the error it suffices to switch the current value of Colid m * , i.e., putting Colid m * = f 2 in the former case and Colid m * = f 1 in the latter. Again, y σ has to be recomputed.

10

 10

2

 2

Figure 1 :

 1 Figure 1: Execution time in seconds (y axis) of Online-Detection (dotted), ABFT-Detection (solid line) and ABFT-Correction (dashed) with respect to the normalized MTBF (mean time between failures), x-axis). The matrix number is in the subcaption.

Table 1 :

 1 Overhead comparison for Algorithm 2 and Algorithm 3. Here n denotes the size of the matrix and n the number of null sum columns.

	Algorithm 2 Algorithm 3

  x ind ;

	33:	end if
	34:	end for
	35:	end for
	36:	c y ← n i=1 y i ; c y ← n i=1 y i ;
	37:

  j ← s j + w ij Rowindex i ; Rowindex i to Rowindex i+1 -1 do i ← y i + Val j • x ind ;

	13:	for i ← 1 to n do
	14:	i ← 0;
	15:	for j = 1 to k do
	16:
	17:	end for
	18: for j ← 19: ind ← Colid j ;
	20:
	21:	end for
	22:	end for
	23:

s y

Table 2 :

 2 Test matrices used in the experiments. Name and id are from the University of Florida Sparse Matrix Collection.However, the other two, having floating point operands, are problematic. Since the floating point operations are not associative and the distributive property does not hold, we need a tolerance parameter that takes into account the rounding operations that are performed by each floating point operation. Here, we give an upper bound on the difference between the two floating point checksums, using the standard model [43, Sec. 2.2] to make sure that errors caught by our algorithms really are errors and not merely inaccuracies due to floating point operations (which is tolerable, as non-faulty executions can give rise to the same inaccuracy).

	name	id	size	density steps residual
	Boeing/bcsstk36	341 23052 2.15e-03	50 6.41e-04
	Mulvey/finan512	752 74752 1.07e-04	25 2.19e-14
	Andrews/Andrews	924 60000 2.11e-04	20 1.59e-04
	GHS psdef/wathen100	1288 30401 5.10e-04	50 2.55e-13
	GHS psdef/wathen120	1289 36441 4.26e-04	50 9.16e-14
	GHS psdef/gridgena	1311 48962 2.14e-04	50 5.61e-05
	GHS psdef/jnlbrng1	1312 40000 1.24e-04	50 5.83e-13
	UTEP/Dubcova2	1848 65025 2.44e-04	50 1.16e-05
	JGD Trefethen/Trefethen 20000 2213 20000 1.39e-03	10 6.00e-16

Theorem 3 (Accuracy of the floating point weighted checksums). Let A ∈ R n×n , x ∈ R n , c ∈ R n . If all of the sums involved into the matrix operations are performed using some flavor of recursive summation

[START_REF] Higham | Accuracy and Stability of Numerical Algorithms, 2nd Edition[END_REF] Ch. 4]

, it holds that |f l (

Table 3 :

 3 Experimental validation of the model. Here s i and s * i represent the best checkpointing interval according to our model and to our simulations respectively, whereas E t ( s i ) and E t s * i stand for the execution time of the algorithm using these checkpointing intervals.

	id s 1 E t ( s 1 ) s * 1	E t (s * 1 )	l 1 s 2 E t ( s 2 ) s * 2	E t (s * 2 )	l 2
	341	4 305.42	1 293.22 4.16	4 305.45	1 293.16	4.19
	752 30	13.81 24	13.34 3.57 24	12.17 23	11.93	2.01
	924 23	49.82 30	47.53 4.82 23	44.52 26	42.42	4.96
	1288 22	11.12 19	10.82 2.72 22	11.32 19	11.03	2.58
	1289 16	16.56 13	16.38 1.07 23	13.49 23	13.49	0.00
	1311	4 216.70	1 207.97 4.19	4 220.20	1 208.19	5.77
	1312 25	14.41 22	13.86 3.97 23	12.30 22	12.06	1.96
	1848	4 321.70	1 309.28 4.01	4 366.03	1 314.20 16.49
	2213 19	2.31 12	2.19 5.58 24	2.33 23	2.20	5.94

Acknowledgements

Y. Robert and B. Uçar were partly supported by the French Research Agency (ANR) through the Rescue and SOLHAR (ANR MONU-13-0007) projects. Y. Robert is with Institut Universitaire de France. J. Langou was fully supported by NSF award CCF 1054864.