Nathann Cohen
email: nathann.cohen@gmail.com

Dmitrii V Pasechnik

Implementing Brouwer's database of strongly regular graphs

Keywords: 05E30, strongly regular graphs, association schemes; 68-04, explicit machine computation and programs

Andries Brouwer maintains a public database of existence results for strongly regular graphs on n ≤ 1300 vertices. We have implemented most of the infinite families of graphs listed there in the open-source software Sagemath [18], as well as provided constructions of the "sporadic" cases, to obtain a graph for each set of parameters with known examples. Besides providing a convenient way to verify these existence results from the actual graphs, it also extends the database to higher values of n.

Introduction

Many researchers in algebraic combinatorics or an adjacent field at some point want to get their hands on a list of feasible parameters of strongly regular graphs, and on actual examples of graphs. These graphs are studied and/or used in hundreds of articles; recent highlights in using strongly regular graphs include A. Bondarenko's [START_REF] Bondarenko | On Borsuk's conjecture for two-distance sets[END_REF] and an improvement of the latter by T. Jenrich and A. E. Brouwer [START_REF] Jenrich | A 64-dimensional counterexample to Borsuk's conjecture[END_REF]. While parameters are available from A. E. Brouwer's online database [START_REF] Brouwer | Parameters of strongly regular graphs[END_REF], actually constructing an example can easily take a lot of time and effort. The project described here aims at making these tasks almost trivial by providing the necessary graph constructions, and a way to obtain a strongly regular graph from a tuple of parameters, in the computer algebra system Sagemath [START_REF] Developers | Sage Mathematics Software (Version 7.3)[END_REF] (also known as Sage). It is worth mentioning that a large part of the project makes use of GAP [START_REF]GAP -Groups, Algorithms, and Programming[END_REF] and its packages, in particular L. Soicher's package GRAPE [START_REF] Soicher | Computing with graphs and groups[END_REF].

Exhaustive enumeration of the non-isomorphic strongly regular graphs has been performed for some tuples T of parameters (see T. Spence [START_REF] Spence | Strongly regular graphs on at most 64 vertices[END_REF]). However, the sheer number of non-isomorphic examples (see e.g. D. Fon-Der-Flaass [START_REF] Fon-Der-Flaass | New prolific constructions of strongly regular graphs[END_REF] or M. Muzychuk [START_REF] Muzychuk | A generalization of Wallis-Fon-Der-Flaass construction of strongly regular graphs[END_REF]) makes it hard to expect to be able to generate all of them, for a given T , in reasonable time. Thus we opted for a minimalist approach: for each T we generate an example, provided that one is known. We note, however, that some of constructions implemented can generate many examples with the same T ; e.g. we have implemented the construction to generate the point graph of the generalized quadrangle T * 2 (O) (see [START_REF] Payne | Finite generalized quadrangles[END_REF]) from any hyperoval O ⊂ P G(2, 2 k ). As well, many T 's can be realized by more than one implemented construction, sometimes leading to isomorphic graphs, and sometimes not.

Our desire to take on this project was motivated by the following considerations.

• One wants to double-check that the constructions are correct and their descriptions are complete; indeed, a program is more trustworthy than a proof in some situations, and coding a construction is a good test for completeness of the description provided.

• We wanted to see that the Sage combinatorial, graph-theoretic, and grouptheoretic primitives to deal with such constructions are mature and versatile, so that coding of constructions is relatively easy and quick.

• One learns a lot while working on such a project, both the underlying mathematics, and how the toolset can be improved. In particular, one might come along simplifications of constructions, and this actually happened on couple of occasions, see Sect. 5.3.

• As time goes by, possible gaps in constructions are harder and harder to fill in. Reconstructing omitted proof details becomes a tricky and timeconsuming task.

In particular, as far as the latter item is concerned, we seem to have uncovered at least one substantial gap in constructions (see Sect. 7). Furthermore, a number of constructions needed feedback from their authors or discussions with experts -sometimes quite substantial -to code them.

A large part of the constructions use in a nontrivial way another combinatorial or algebraic object: block design, Hadamard matrix, two-graph, twodistance code, finite group, etc. In particular, at the start of the project some of these were lacking in Sagemath, we needed to implement constructions of certain block designs, regular symmetric Hadamard matrices with constant diagonal (where the gap just mentioned was uncovered), skew-Hadamard matrices, and two-graphs. As well, we created a small database of two-distance codes (see Sect. 4.2).

The remainder of the paper consists of a short introduction to strongly regular graphs, pointing out particular relevant Sagemath features, and a description of our implementations, few of them seemingly novel. We list the constructions that we implemented, and discuss few gaps we discovered in the literature.

Strongly regular graphs and related objects

An undirected regular degree k graph Γ on n vertices (with 0 < k < n -1) is called strongly regular if the vertices u and v of any edge have λ common neighbours, and the vertices u and v of any non-edge have µ common neighbours. One says that Γ has parameters (n, k, λ, µ). Note that the complement of Γ, i.e. the graph with the same set of vertices and edges being precisely the non-edges of Γ, is also a strongly regular whose parameters are related by a simple formula to these of Γ (see e.g. A. E. Brouwer and W. Haemers [START_REF] Brouwer | Spectra of graphs[END_REF] for details).

Example 2.1 Let Γ be the graph with vertices being k-subsets of an m-set, with k ≤ m/2 ; two vertices are adjacent if the corresponding k-subsets intersect in a (k -1)-subset. Such graphs are called Johnson graphs and denoted by J(m, k) (in Sagemath, J(m, k) can be constructed by calling the function graphs.JohnsonGraph(m,k)). Then J(m, 2) is a strongly regular graph, with parameters

( m 2 , 2(m -2), m -2, 4).
Example 2.2 RSHCD -a (n, )-Regular Symmetric Hadamard matrix M with Constant Diagonal is an n × n symmetric ±1-matrix such that: 1) M M T = nI; 2) its rows sums are all equal to δ √ n, where ∈ {-1, +1} and δ is the (constant) diagonal value of M , usually denoted RSHCD -and RSHCD + . These matrices yield regular two-graphs. As well, they yield strongly regular graphs: replacing all the entries equal to the diagonal values by 0, and the remaining entries by 1 gives the adjacency matrix of a strongly regular graph. Some sources further require that both Γ and its complement are connected; in terms of parameters this means 0 < µ < k. This excludes the trivial case of Γ (or its complement) being disjoint union of complete graphs of the same size. Sagemath implementation does not impose this restriction.

A considerable number of techniques ruling out the existence of a strongly regular graph Γ with given parameters (n, k, λ, µ) are known, e.g. based on computing eigenvalues of the adjacency matrix A of Γ. As A generates a dimension 3 commutative subalgebra of C n×n , one sees that there are just 3 distinct eigenvalues of A, and they are determined by the parameters (e.g. the largest eigenvalue is k). Sagemath implements parameter-based techniques to rule out sets of parameters from A. E. Brouwer and J. H. van Lint [START_REF] Brouwer | Strongly regular graphs and partial geometries[END_REF], and from A. E. Brouwer, A. M. Cohen, and A. Neumaier [START_REF] Brouwer | Distance-regular graphs[END_REF].

We use standard terminology for finite permutation groups, finite simple groups, and geometries over finite fields from [START_REF] Brouwer | Distance-regular graphs[END_REF][START_REF] Brouwer | Spectra of graphs[END_REF].

Structure and use of the implementation

The strongly regular graphs are split into two categories: the fixed-size graphs (see Sect.4) and the families of strongly regular graphs (see Sect.5). The parameters (n, k, λ, µ) of fixed-size graphs are hardcoded, while each family of strongly regular graphs has a helper function which takes (n, k, λ, µ) as an INPUT and answers whether the graph family is able to produce a graph with the required parameters. Some families forward their queries to the databases of Balanced Incomplete Block Designs, of Orthogonal Arrays, of Hadamard matrices of various types, and of 2-weight codes.

With this design, it takes ≤ 3 seconds on a modern laptop to know which graphs on < 1300 vertices can be produced by the implemented constructions (i.e. as far as the online database goes).

In Sage, everything is made available to the user through a single function graphs.strongly regular graph that produces a graph matching the provided parameters. Note that µ, the fourth parameter, can be omitted. sage: G=graphs.strongly regular graph(175,30,5) sage: G AS(5)*; GQ [START_REF] Brouwer | Uniqueness and nonexistence of some graphs related to M 22[END_REF][START_REF] Brouwer | Parameters of strongly regular graphs[END_REF]: Graph on 175 vertices One can use the same function to learn whether a set of parameters is realizable, if it is not, or if the existence problem is unsettled: sage: graphs.strongly regular graph(175,30,5,5,existence=True) True sage: graphs.strongly regular graph(57,14,1,existence=True) False sage: graphs.strongly regular graph(3250,57,0,1,existence=True) Unknown Technical details and descriptions of many specific functions is available as a part [START_REF] Cohen | Sagemath: Database of strongly regular graphs[END_REF] of the Sagemath manual, which can be found online [START_REF] Developers | Sage Mathematics Software (Version 7.3)[END_REF].

4 Fixed-size constructions

"Sporadic" examples

Here we did not attempt to give an exhaustive list of references for each graph, for some of them have several papers devoted to them in one or another way.

We identify the corresponding graphs by their parameters, and provide references and some construction details for each of them. [START_REF] Ionin | Combinatorics of symmetric designs[END_REF][START_REF]Handbook of combinatorial designs[END_REF][START_REF] Brouwer | Parameters of strongly regular graphs[END_REF][START_REF] Brouwer | Uniqueness and nonexistence of some graphs related to M 22[END_REF] Hubaut [33, S.9]. Subgraph of common neighbours of a triangle in Suzuki graph. (100, 44, 18, 20) Jørgensen and Klin [START_REF] Jørgensen | Switching of edges in strongly regular graphs. I. A family of partial difference sets on 100 vertices[END_REF]. Built as a Cayley graph.

(100, 45, 20, 20) [START_REF] Jørgensen | Switching of edges in strongly regular graphs. I. A family of partial difference sets on 100 vertices[END_REF]. Built as a Cayley graph.

(105, 32, 4, 12) Goethals and Seidel [START_REF] Goethals | Strongly regular graphs derived from combinatorial designs[END_REF], Coolsaet [START_REF] Coolsaet | The uniqueness of the strongly regular graph srg(105, 32, 4, 12[END_REF].

(120, 63, 30, 36) R. (176, 49, 12, 14) Brouwer [START_REF] Brouwer | Polarities of G. Higman's symmetric design and a strongly regular graph on 176 vertices[END_REF].Built from the symmetric Higman-Sims design.

There exists an involution σ exchanging the points and blocks of the Higman-Sims design, such that each point is mapped onto a block that contains it (i.e. σ is a polarity with all absolute points). The graph is then built by making two vertices u, v adjacent whenever v ∈ σ(u).

( We then build the intersection graph of blocks with intersection size 3.

(196, 91, 42, 42) Ionin and Shrikhande [START_REF] Ionin | Combinatorics of symmetric designs[END_REF].

(210, 99, 48, 45) Klin et al. [START_REF] Klin | Examples of computer experimentation in algebraic combinatorics[END_REF].S 7 acts on the 210 digraphs isomorphic to the disjoint union of K 1 and the circulant 6-vertex digraph which one can obtain using Sagemath as digraphs.Circulant(6, [START_REF] Abel | Existence of five MOLS of orders 18 and 60[END_REF][START_REF] Brouwer | Parameters of strongly regular graphs[END_REF]). This action has 16 orbitals; the package [START_REF] Faradjev | Computer package for computations with coherent configurations[END_REF] found a merging of them, explicitly described in [START_REF] Klin | Examples of computer experimentation in algebraic combinatorics[END_REF], resulting in this graph.

(231, 30, 9, 3) Brouwer [START_REF] Brouwer | Uniqueness and nonexistence of some graphs related to M 22[END_REF]. The Cameron graph.

(243, 110, 37, 60) Goethals and Seidel [START_REF] Goethals | The regular two-graph on 276 vertices[END_REF].Consider the orthogonal complement of the ternary Golay code, which has 243 words. On them we define a graph, with two words adjacent if their Hamming distance is 9.

( (1800,1029,588, 588) Janko and Kharaghani [START_REF] Janko | A block negacyclic Bush-type Hadamard matrix and two strongly regular graphs[END_REF].

Two-weight codes database

The rest of the fixed-size constructions of strongly regular graphs in the database originate from linear d-dimensional two-weight codes of length with weights w 1 and w 2 over F q . We use data shared by Eric Chen [START_REF] Chen | Online database of two-weight codes[END_REF], data by Axel Kohnert [START_REF] Kohnert | Constructing two-weight codes with prescribed groups of automorphisms[END_REF] shared by Alfred Wassermann, data from I. 96 [START_REF] Chen | Online database of two-weight codes[END_REF] Note that some of these codes are members of infinite families; this will be explored and extended in forthcoming work.

Infinite families

These are roughly divided into two parts: graphs related to finite geometries over finite fields (in particular various classical geometries), and graphs obtained by combinatorial constructions.

Graphs from finite geometries

Here q denotes a prime power, and ∈ {-, +}.

• Graphs arising from projective geometry designs are discussed in Sect. 5.2, along with other Steiner graphs.

• Paley graphs. The vertices are the elements of F q , with q ≡ 1 mod 4; two vertices are adjacent if their difference is a nonzero square in F q ; see [8, 9.1.2].

• Polar space graphs. These include polar spaces for orthogonal and unitary groups, see entries O 2d (q), O 2d+1 (q), and U d (q) in [8, Table 9.9]. Sagemath also has an implementation of polar spaces for symplectic groups (entry Sp 2d (q) in [loc.cit.]), but we do not use them in the database, as they have the same parameters as these for orthogonal groups.

• Generalised quadrangle graphs, GQ(s, t) in [8, Table 9.9]. Apart from these appearing as polar space graphs, with s = t = q, s 2 = t = q, and s = q 2 , t = q 3 , we provide other examples, as follows.

-Unitary dual polar graphs. This gives s = q 3 , t = q 2 .

-GQ(q -1, q + 1)-graphs for q odd are constructed following Ahrens and Szekeres, see [49, 3.1.5], and for q even we provide the T * 2 (O) construction, see [49, 3.1.3], from a hyperoval O in PG(2, q).

-GQ(q + 1, q -1) are constructed as line graphs of GQ(q -1, q + 1).

• Affine polar graphs. These are the entry VO 2d (q) in [8, Table 9.9].

• Graphs of non-degenerate hyperplanes of orthogonal polar spaces, with adjacency specified by degenerate intersection; see NO 2d+1 (q) in [8, Table 9.9]. These are constructions by Wilbrink, cf. [10, Sect.7.C]. The implementation in Sagemath simply takes the appropriate orbit and orbital of the orthogonal group acting on the hyperplanes using parameters of the graph, namely v = q d (q d + )/2, k = (q d -)(q d-1 + ).

• Graphs of non-isotropic points of polar spaces, with adjacency specified by orthogonality. These include a number of cases.

-Non-isotropic points of orthogonal polar spaces over F 2 ; see NO 2d [START_REF] Bondarenko | On Borsuk's conjecture for two-distance sets[END_REF] in [8, Table 9.9].

-One class of non-isotropic points of orthogonal polar spaces over F 3 ; see NO 2d (3) in [8, Table 9.9].

-One class of non-isotropic points of orthogonal polar spaces (specified by a non-degenerate quadratic form F ) over F 5 ; see NO ⊥ 2d+1 (5) in [8, Table 9.9]. This is a construction by Wilbrink, cf. [START_REF] Brouwer | Strongly regular graphs and partial geometries[END_REF]Sect.7.D], where the class of points p is described in terms of the type of the quadric specified by p ⊥ ∩ Q, where Q is the set of isotropic points of the space, i.e. Q := {x ∈ P G(2d, 5) | F (x) = 0}, and

p ⊥ := {x ∈ Q | F (p + x) = F (p)}.
The implementation in Sagemath takes {x ∈ P G(2d, 5) | F (x) = ±1} for = +, and the rest of non-isotropic points for = -.

-Non-isotropic points of unitary polar spaces; see NU d(q) in [8, Table 9.9].

• Graphs of Taylor two-graphs, see [8, Table 9.9] and [START_REF] Brouwer | Strongly regular graphs and partial geometries[END_REF]Sect.7E]. Note that we implement an efficient construction that does not need all the triples of the corresponding two-graphs, by first directly constructing the descendant strongly regular graphs on q 3 vertices, and a partition of its vertices into cliques. The latter provides a set to perform Seidel switching on the disjoint union with K 1 , and obtain the strongly regular graph on q 3 + 1 vertices. See Sagemath documentation for graphs.TaylorTwographSRG for details.

• Cossidente-Penttila hemisystems in PG(5, q), for q odd prime power [START_REF] Cossidente | Hemisystems on the Hermitian surface[END_REF], are certain partitions of points of the minus type quadric in PG(5, q) into two parts V , V of equal size. The subgraph Γ of the collinearity graph of the corresponding GQ(q, q 2 ) induced on V has parameters ((q 3 + 1)(q + 1)/2, (q 2 + 1)(q -1)/2, (q -3)/2, (q -1) 2 /2). The way we construct these graphs in Sage is described in Sect. 5.3.1.

Graphs from combinatorics

• Johnson Graphs J(m, 2), see Example 2.1.

• Orthogonal Array block graphs OA(k, n). Sage is able to build a very substantial state-of-the-art collection of orthogonal arrays (often abbreviated as OA), thanks to a large implementation project undertaken in 2013/2014 by the first author in a very productive collaboration with Julian R. Abel and Vincent Delecroix. For the present work no new constructions of OAs were needed, and the link between Sage's OAs and Strongly Regular Graphs databases filled in three new entries in Andries Brouwer's database.

• Steiner Graphs (intersection graphs of BIBD) -Sage can already build several families of Balanced Incomplete Block Designs (when k ≤ 5, or projective planes, or other recursive constructions and fixed-size instances). More constructions from [START_REF]Handbook of combinatorial designs[END_REF] were added to Sage while working on this project.

• Goethals-Seidel graphs, see [START_REF] Goethals | Strongly regular graphs derived from combinatorial designs[END_REF].

• Haemers graphs, see [START_REF] Brouwer | Strongly regular graphs and partial geometries[END_REF]Sect.8.A].

• RSHCD -graphs from (n, )-regular symmetric Hadamard matrices M with constant diagonal, see Example 2.2 for the definition. Several constructions from the literature (and one apparently new one, cf. Sect. 5.3.2)

for this class of Hadamard matrices were implemented in Sage and are available in its Hadamard matrices module.

• Two-graph descendants. Each regular two-graph (a certain class of 3uniform v-vertex hypergraphs having 2µ three-edges on each pair of points, cf. e.g. [START_REF] Brouwer | Spectra of graphs[END_REF]Chap.10]) gives rise to a strongly regular graph with parameters (v -1, 2µ, 3µ -v/2, µ) obtained by descendant construction, see e.g. [START_REF] Brouwer | Spectra of graphs[END_REF]Sect.10.3].

• Switch OA Graphs -these strongly regular graphs are obtained from OA block graphs (see above). From such a graph G obtained from an OA(k, n), the procedure is to (1) add a new isolated vertex v;

(2) perform Seidel switching on the union of {v} and several disjoint n-cocliques of G. Note that a n-coclique in G corresponds to a parallel class of the OA(k, n), and that those are easily obtained from an OA(k + 1, n) (i.e. a resolvable OA(k, n)).

• Polhill Graphs -In [START_REF] Polhill | Negative Latin square type partial difference sets and amorphic association schemes with Galois rings[END_REF], Polhill produced 5 new strongly regular graphs on 1024 vertices as Cayley graphs. His construction is able to produce larger strongly regular graphs of order ≥ 4096, though the current implementation only covers the n = 1024 range.

• Mathon's pseudo-cyclic strongly regular graphs related to symmetric conference matrices, optionally parameterised by a strongly regular graph with parameters of a Paley graph, and a skew-symmetric Latin square [START_REF] Mathon | Symmetric conference matrices of order pq 2 + 1[END_REF][START_REF] Seidel | Two-graphs, a second survey[END_REF].

• Pseudo-Paley and Pasechnik graphs from skew-Hadamard matrices. These are constructions due to Goethals-Seidel [START_REF] Brouwer | Strongly regular graphs and partial geometries[END_REF] and Pasechnik [START_REF] Pasechnik | Skew-symmetric association schemes with two classes and strongly regular graphs of type L 2n-1 (4n -1)[END_REF], constructing graphs on (4m -1) 2 vertices from skew-Hadamard matrices of order 4m. Sage builds the corresponding skew-Hadamard matrices from a small database featuring classical constructions of skew-Hadamard matrices from [START_REF] Hall | Combinatorial theory[END_REF] and small examples from (anti)-circulant matrices [START_REF] Goethals | A skew Hadamard matrix of order 36[END_REF][START_REF] Wallis | A skew-Hadamard matrix of order 92[END_REF].

Novel constructions

Here we collect descriptions of constructions of graphs that are in our view suffciently novel and interesting to mention. Namely, Sect. 5.3.1 describes another construction for a known graph, Sect. 5.3.2 describes a working construction for a graph which was claimed to exist in the literature, although we were unable to verify a number of published constructions (see Sect. 7 for details). Finally, Sect. 5.3.3 discusses an unpublished construction by Ionin and Kharaghani.

Cossidente-Penttila hemisystems

The construction of the hemisystem in [START_REF] Cossidente | Hemisystems on the Hermitian surface[END_REF] requires building GQ(q 2 , q), which is slow. Thus we designed, following a suggestion of T. Penttila, a more efficient approach, working directly in PG(5, q). The partition in question is invariant under the subgroup H = Ω 3 (q 2 ) < O - 6 (q). Without loss in generality H leaves the form B(X, Y, Z) = XY + Z 2 invariant. We pick two orbits of H on the F q -points, one of them B-isotropic, with a representative (1 : 0 : 0), viewed as a point of Π := PG(2, q 2 ), and the other corresponding to points of Π that have all the lines on them intersecting the conic of Π specified by B in zero or two points. We take (1 : 1 : ) as a representative, with ∈ F * q 2 so that 2 + 1 is not a square in F q 2 . Indeed, the conic can be viewed {(0

: 1 : 0)} ∪ {(1 : -t 2 : t) | t ∈ F q 2 }.
The coefficients of a generic line on (1 : 1 : ) are [1 : -1 -b : b], for -1 = b. Thus, to make sure that its intersection with the conic is always even, we need that the discriminant of 1 + (1 + b)t 2 + tb = 0 never vanishes, and this is if and only if 2 + 1 is not a square.

Finally, we need to adjust B, by multiplying it by appropriately chosen ν ∈ F * q 2 , so that (1 : 1 : ) becomes isotropic under the relative trace norm (X : Y : Z) → νB(X, Y, Z) + (νB(X, Y, Z)) q , used to define adjacency in Γ.

Regular symmetric Hadamard matrices of order 324.

We recall the definition of RSHCD + and RSHCD -from Example 2.2. An example M + of RSHCD + order 324 was constructed by Janko, Kharaghani, and Tonchev in [START_REF] Janko | The existence of a Bush-type Hadamard matrix of order 324 and two new infinite classes of symmetric designs[END_REF], and we implemented their construction in Sagemath. See [START_REF] Muzychuk | Symmetric Bush-type Hadamard matrices of order 4m 4 exist for all odd m[END_REF][START_REF] Haemers | Strongly regular graphs with parameters (4m 4 , 2m 4 + m 2 , m 4 + m 2 , m 4 + m 2 ) exist for all m > 1[END_REF] for other examples of RSHCD + of order 324.

We use M + to build an example M -of RSHCD -of order 324, as follows. One is tempted to apply [START_REF] Haemers | Strongly regular graphs with parameters (4m 4 , 2m 4 + m 2 , m 4 + m 2 , m 4 + m 2 ) exist for all m > 1[END_REF]Lemma 11] to M + , which says that for an RSHCD matrix M built from four n × n-blocks M ij , so that

M = M 11 M 12 M 21 M 22 , the matrix T (M ) := M 11 -M 12 -M 21 M 22 (1)
is an RSHCD -, provided that row sums of M 11 and M 22 are 0. However, the latter condition does not hold for M = M + . We are able to "twist" M + so that the resulting matrix is amenable to this Lemma. Namely, it turns our that the matrix

M := M 12 M 11 M 11 M 21 ,
where M = M + , is RSHCD + , its diagonal blocks having row sums 0, as needed by [START_REF] Abel | Existence of five MOLS of orders 18 and 60[END_REF]. Interestingly, the (324, 152, 70, 72)-strongly regular graph corresponding to T (M ) has a vertex-transitive automorphism group of order 2592, twice the order of the (intransitive) automorphism group of the (324, 153, 72, 72)-strongly regular graphcorresponding to M + . As far as we know, this is the only known example of such a vertex-transitive graph. Other graphs with such parameters were constructed in [START_REF] Haemers | Strongly regular graphs with parameters (4m 4 , 2m 4 + m 2 , m 4 + m 2 , m 4 + m 2 ) exist for all m > 1[END_REF].

A (765,192,48,48)-graph

We were unable to implement the construction of a graph with these parameters described, as a part of an infinite family, in Ionin and Kharaghani [START_REF] Ionin | New families of strongly regular graphs[END_REF]. The authors of the latter were very kind to send us an updated construction of the graph in question, which we successfully implemented in Sagemath, see [START_REF] Cohen | graphs.IoninKharaghani765Graph. Sagemath Trac ticket[END_REF]. This construction can be found in the documentation of the Sagemath function graphs.IoninKharaghani765Graph.

They have also posted an update [START_REF] Ionin | New families of strongly regular graphs[END_REF] to [START_REF] Ionin | New families of strongly regular graphs[END_REF]; we have not yet tried to implement the updated version in full generality.

Missing values

Among the 1150 realizable, according to Andries Brouwer's database, parameter sets, our implementation can realize 1142. Up to taking graph complements, the list of currently missing entries is as follows.

(196 90 40 42) RSHCD -(may not exist, cf. Sect. 7) (196 135 94 90) Huang, Huang and Lin [START_REF] Huang | On a class of strongly regular designs and quasi-semisymmetric designs[END_REF] (378 116 34 36) Muzychuk S6(n = 3, d = 3) [START_REF] Muzychuk | A generalization of Wallis-Fon-Der-Flaass construction of strongly regular graphs[END_REF] (512 133 24 38) Godsil (q = 8, r = 3) [START_REF] Godsil | Krein covers of complete graphs[END_REF] The fisrt entry is discussed in Sect. 7. Implementation of the remaining three entries is currently in progress.

Incorrect RSHCD constructions ?

We were unable to reproduce the following two constructions of Regular Symmetric Hadamard matrices with Constant Diagonal (RSHCDs) and thus the corresponding strongly regular graphs (in the sense of [START_REF] Brouwer | Spectra of graphs[END_REF]Sect.10.5]).

• In [8, Sect.10.5.1, (i)], the construction of RSHCD(196, -) is attributed to [START_REF] Ionin | Combinatorics of symmetric designs[END_REF], in which the existence of a (4k If n -1 and n + 1 are odd prime powers, there exists a RSHCD(n 2 , +).

We implemented the construction provided in [58, Corollary 5.12], but that did not lead us to the expected strongly regular graph. We also note that while Corollary 5.12 does not claim that the provided matrices are regular, that claim appears in the theorem on which it relies. The author of [START_REF] Wallis | Combinatorics: Room squares, sum-free sets, Hadamard matrices[END_REF]Corollary 5.12] did not answer our message, and we discarded this construction as broken in our work.

This construction should have been able to produce a RSHCD(676, +) and a RSHCD(900, +). Fortunately in the end it was not required.

The following construction of RSHCDs needed a lot of effort and a number of discussions with Andries Brouwer to correct crucial misprints in several sources and combine them into a working construction.

• In [8, Sect.10.5.1, (iv)] one finds the following claim, attributed to [START_REF] Wallis | Combinatorics: Room squares, sum-free sets, Hadamard matrices[END_REF]Corollary 5.16].

If a + 1 is a prime power and there exists a symmetric conference matrix of order a, then there exists a RSHCD(a 2 , +).

Following [loc.cit.] did not lead us to the expected result; as it turns out, [loc.cit.] has a typo, and the correct formulae should be taken from the original source [57, Corollary 17] by J. Wallis and A. L. Whiteman. An essential ingredient in this construction, referred to in [START_REF] Wallis | Some classes of Hadamard matrices with constant diagonal[END_REF], is a special pair of difference sets due to G. Szekeres [START_REF] Szekeres | Tournaments and Hadamard matrices[END_REF]Theorem 16], defined in (4.1) and (4.2) there. However, (4.2) has a typo (it has -instead of + sign), invalidating the construction. Fortunately, a correct definition for may be found [START_REF] Wallis | Combinatorics: Room squares, sum-free sets, Hadamard matrices[END_REF]Theorem 2.6,p.303]. This construction allowed us to produce RSHCD(676, +) and a RSHCD(900, +).

Entries added to database during this work

By linking Sage's database of Orthogonal Arrays with its database of Strongly Regular Graphs, we were able to fill in the following three values:

• (196, 78, 32, 30) -can be obtained from an OA(6, 14) [START_REF] Todorov | Four mutually orthogonal Latin squares of order 14[END_REF] • (324, 102, 36, 30) -can be obtained from an OA(6, 18) [START_REF] Abel | Existence of five MOLS of orders 18 and 60[END_REF] • (324, 119, 46, 42) -can be obtained from an OA(7, 18) [START_REF] Abel | Existence of five MOLS of orders 18 and 60[END_REF] This can be seen as a by-product of making two mathematical databases, which formerly only existed in printed form, inter-operable. In our implementation, any update of the combinatorial designs databases can be beneficial for the database of strongly regular graphs.

We obtained a (1024, 462, 206, 210)-graph while going through the constructions from [START_REF] Polhill | Negative Latin square type partial difference sets and amorphic association schemes with Galois rings[END_REF], although this value did not appear in the online database at that time.

  Sect.9.1.7 (vi)]. The M 22 -graph.(100,[START_REF]GAP -Groups, Algorithms, and Programming[END_REF] 0,[START_REF] Brouwer | Uniqueness and nonexistence of some graphs related to M 22[END_REF] [8, Sect.9.1.7 (vii)]. The Higman-Sims graph.

  Mathon, cf. [10, Sect.6.A]. The distance-2 graph of J(10, 3).(120, 77, 52, 44) Unique by J. Degraer K. Coolsaet[START_REF] Degraer | Classification of some strongly regular subgraphs of the mclaughlin graph[END_REF].We first build a 2 -(21, 7, 12) design, by removing two points from the Witt design on 23 points.We then build the intersection graph of blocks with intersection size 3. 144, 39, 6, 12) A. A. Ivanov, M. H. Klin, and I. A. Faradjev[START_REF] Ivanov | The primitive representations of the non-abelian simple groups of order less than 10 6 , part II[END_REF] Table 9].An orbital of degree 39 (among 2 such orbitals) of the group P SL 3 (3) acting on the (right) cosets of a subgroup of order 39.

	(126, 25, 8, 4)	R. Mathon, cf. [10, Sect.6.A].The distance-(1 or 4) graph of
	J(9, 4).
	(126, 50, 13, 24) Goethals, cf. [10].

((162, 56, 10, 24) [33, S.12].The complement of the subgraph induced on the neighbours of a vertex in the complement of McLaughlin graph. (175, 72, 20, 36) [10, Sect.10.B (iv)].Obtained from the line graph Λ of Hoffman-Singleton Graph, by setting two vertices to be adjacent if their distance in Λ is exactly 2. For more information, see http://www.win.tue.nl/ ~aeb/graphs/ McL.html.

  253, 140, 87, 65) [33, S.6]; a rank 3 representation of M 23 .We first build the Witt design on 23 points which is a 2 -[START_REF] Godsil | Krein covers of complete graphs[END_REF][START_REF] Brouwer | Distance-regular graphs[END_REF][START_REF] Fon-Der-Flaass | New prolific constructions of strongly regular graphs[END_REF] design. We then build the intersection graph of blocks with intersection size 3. Mathon and Rosa[START_REF] Mathon | A new strongly regular graph[END_REF].The vertices of the graph are all 280 partitions of a set of cardinality 9 into 3-sets, e.g. {{a, b, c}, {d, e, f }, {g, h, i}}. The cross-intersection of two partitions P = {P 1 , P 2 , P 3 } and P = {P 1 , P 2 , P 3 } being defined as {P i ∩P j : 1 ≤ i, j ≤ 3}, two vertices of 'G' are set to be adjacent if the cross-intersection of their respective partitions does not contain exactly 7 nonempty sets. 784, 243, 82, 72) R. Mathon, cf. [10, Sect.6.D].This and the following two are Mathon's graphs from merging classes in the product of pseudo-cyclic association scheme for action of O 3 (8) on elliptic lines in PG[START_REF] Bondarenko | On Borsuk's conjecture for two-distance sets[END_REF][START_REF] Brouwer | Spectra of graphs[END_REF], studied by H. D. L. Brouwer and van Eijl[START_REF] Brouwer | On the p-rank of the adjacency matrices of strongly regular graphs[END_REF].This graph is built on the words of weight 12 in the binary Golay code. Two of them are then made adjacent if their symmetric difference has weight 12.(1782,416, 100, 96)[33, S.15]. Suzuki graph, rank 3 representation of Suz.

	(275, 112, 30, 56) [33, S.13]. The McLaughlin graph.
	(276, 140, 58, 84) Haemers and Tonchev [28].The graph is built from from
	McLaughlin graph, with an added isolated vertex. We then perform Seidel
	switching on a set of 28 disjoint 5-cliques.
	(280, 117, 44, 52) Hollmann [31].
	(784, 270, 98, 90) R. Mathon, cf. [10, Sect.6.D].
	(784, 297, 116, 110) R. Mathon, cf. [10, Sect.6.D].
	(936, 375, 150, 150) Janko and Kharaghani [38].
	(1288,792, 476, 504)

(280, 135, 70, 60)

[37,

Table 9, p.51].This graph is built from the rank 4 action of J 2 on the cosets of a subgroup 3.P GL(2, 9). (324, 152, 70, 72) See Sect. 5.3.2.RSHCD -, see Example 2.2. We build an apparently new example using the (324, 153, 72, 72)-graph; other example may be found in [29]. (324, 153, 72, 72) See Sect. 5.3.2.RSHCD + , see Example 2.2. We build the example from [39]; more examples may be found in [47, 29]. (416, 100, 36, 20) [33, S.14]; (rank 3 representation of G 2 (4)).This graph is isomorphic to the subgraph of the Suzuki graph [33, S.15] induced on the neighbours of a vertex. (560, 208, 72, 80) [37, Table 9, p.45].Obtained as the union of 4 orbitals (among the 13 that exist) of the group Sz(8) in its primitive action on 560 points. (630, 85, 20, 10) W. Haemers [27], see also [10, Sect.10.B.(v)].This graph is the line graph of a pg(5, 18, 2); its point graph is the (175, 72, 20, 36)-srg from this table. One then selects a subset of 630 maximum cliques in the latter to form the set of lines of the pg(5, 18, 2). (765, 192, 48, 48) Ionin and Kharaghani, see Sect. 5.3.3. (

  Bouyukliev and J. Simonis [3, Theorem 4.1], and from L. A. Disset [19].

	Graph parameters		Code parameters	Ref.
	n	k	λ	µ q		d	w 1	w
	81	50	31	30 3	15	4	9	[11]
	243 220 199 200 3	55	5	36	[11]
	256 153	92	90 4	34	4	24	[11]
	256 170 114 110 2	85	8	40	[11]
	256 187 138 132 2	68	8	32	[11]
	512	73	12	10 2 219	9	96 [11]
	512 219 102	84 2	73	9	32	[11]
	512 315 202 180 2	70	9	32	[43]
	625 364 213 210 5	65	4	50	[11]
	625 416 279 272 5	52	4	40	[11]
	625 468 353 342 5	39	4	30	[3]
	729 336 153 156 3 168	6 108 [19]
	729 420 243 240 3 154	6	99 [11]
	729 448 277 272 3 140	6	90	[43]
	729 476 313 306 3 126	6	81	[11]
	729 532 391 380 3	98	6	63	[11]
	729 560 433 420 3	84	6	54	[11]
	729 616 523 506 3	56	6	36	[11]
	1024 363 122 132 4 121	5	88	[19]
	1024 396 148 156 4 132	5	96 [19]
	1024 429 176 182 4 143	5 104 [19]
	1024 825 668 650 2 198 10

  Communication with the authors did not solve the issue, and we are not aware of any other proof of the existence of a (196, 90, 40, 42)-strongly regular graph.

	• In [8, Sect.10.5.1, (iii)] one finds the following claim, attributed to [58,
	Corollary 5.12].

2

, 2k 2 + k, k 2 + k)-strongly regular graph, equivalent to a RSHCD(196, -) for k = 7, is claimed in Theorem 8.2.26.(iii). The latter says that the RSHCD(196, -) can be easily obtained from the RSHCD(196, +) from

[START_REF] Ionin | Combinatorics of symmetric designs[END_REF] Theorem 8.2.26.(ii)

]. While the construction of (ii) was successfully implemented in Sage, following the authors' instructions for (iii) did not lead us to the RSHCD(196, -).

Acknowledgments

The authors thank Andries Brouwer, Eric Chen, Luis Disset, Hadi Kharaghani, Misha Muzychuk, Tim Penttila, John B. Polhill, Leonard Soicher, Vladimir Tonchev, and Alfred Wasserman for many helpful discussions and communications.

The second author was partially supported by the EU Horizon 2020 research and innovation programme, grant agreement OpenDreamKit No 676541.


Guillaume Aupy
email: aupy@cse.psu.edu

Anne Benoit
email: anne.benoit@ens-lyon.fr

Aurélien Cavelan
email: aurelien.cavelan@ens-lyon.fr

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

Yves Robert
email: yves.robert@ens-lyon.fr

Hongyang Sun
email: hongyang.sun@ens-lyon.fr

Bora Uc
email: bora.ucar@ens-lyon.fr

Coping with silent errors in HPC applications

This chapter describes a unified framework for the detection and correction of silent errors, which constitute a major threat for scientific applications at extreme-scale. We first motivate the problem and explain why checkpointing must be combined with some verification mechanism. Then we introduce a generalpurpose technique based upon computational patterns that periodically repeat over time. These patterns interleave verifications and checkpoints, and we show how to determine the pattern minimizing expected execution time. Then we move to application-specific techniques and review dynamic programming algorithms for linear chains of tasks, as well as ABFT-oriented algorithms for iterative methods in sparse linear algebra.

Introduction

For High-Performance Computing (HPC) applications, scale is a major opportunity. Massive parallelism with 100,000+ nodes is the most viable path to achieving sustained petascale performance. Future platforms will exploit even more computing resources to enter the exascale era.

Unfortunately, scale is also a major threat, because resilience becomes a key challenge. Even if each node provides an individual MTBF (Mean Time Between Failures) of, say, one century, a machine with 100,000 such nodes encounters on average a failure every 9 hours, an interval much shorter than the execution time of many HPC applications. Note that (i) a one-century MTBF per node is an optimistic figure, given that each node features several hundreds of cores; and (ii) in some scenarios for the path to exascale computing [START_REF] Dongarra | The international exascale software project: a call to cooperative action by the global high-performance community[END_REF], one envisions platforms including up to one million such nodes, whose MTBF will decrease to 52 minutes.

Several kinds of errors need to be considered when computing at scale. In the recent years, the HPC community has traditionally focused on fail-stop errors, such as hardware failures. The de facto general-purpose technique to recover from failstop errors is checkpoint/restart [START_REF] Chandy | Distributed snapshots: Determining global states of distributed systems[END_REF][START_REF] Elnozahy | A survey of rollback-recovery protocols in message-passing systems[END_REF]. This technique employs checkpoints to periodically save the state of a parallel application, so that when an error strikes some process, the application can be restored into one of its former states. There are several families of checkpointing protocols, but they share a common feature: each checkpoint forms a consistent recovery line, i.e., when an error is detected, one can rollback to the last checkpoint and resume execution, after a downtime and a recovery time. Many models are available to understand the behavior of the checkpointing and restarting techniques [START_REF] Bougeret | Checkpointing strategies for parallel jobs[END_REF][START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF][START_REF] Ozaki | Distribution-free checkpoint placement algorithms based on min-max principle[END_REF][START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF].

While the picture is quite clear for fail-stop errors, the community has yet to devise an efficient approach to cope with silent errors, primary source of silent data corruptions. Such errors must also be accounted for when executing HPC applications [START_REF] Moody | Design, Modeling, and Evaluation of a Scalable Multi-level Checkpointing System[END_REF][START_REF] O'gorman | The effect of cosmic rays on the soft error rate of a DRAM at ground level[END_REF][START_REF] Ziegler | Accelerated testing for cosmic soft-error rate[END_REF][START_REF] Ziegler | Cosmic ray soft error rates of 16-Mb DRAM memory chips[END_REF][START_REF] Ziegler | IBM experiments in soft fails in computer electronics[END_REF]. They may be caused, for instance, by soft errors in L1 cache, arithmetic errors in the ALU (Arithmetic and Logic Unit), or bit flips due to cosmic radiation. The main issue is that the impact of silent errors is not immediate, since they do not manifest themselves until the corrupted data impact the result of the computation (see Figure 1), leading to a failure. If an error striking before the last checkpoint is detected after that checkpoint, then the checkpoint is corrupted, and cannot be used to restore the application. If only fail-stop failures are considered, a checkpoint cannot contain a corrupted state, because a process subject to failure cannot create a checkpoint or participate to the application: failures are naturally contained to failed processes. When dealing with silent errors, however, faults can propagate to other processes and checkpoints, because processes continue to participate and follow the protocol during the interval that separates the occurrence of the error from its detection.

In Figure 1, X s and X d are random variables that represent the time until the next silent error and its detection latency, respectively. We usually assume that silent errors strike according to a Poisson process of parameter λ , so that X s has the distribution of an exponential law of parameter λ and mean 1/λ . On the contrary, it is very hard to make assumptions on the distribution of X d . To alleviate the problem of detection latency, one may envision to keep several checkpoints in memory, and to restore the application from the last valid checkpoint, thereby rolling back to the last correct state of the application [START_REF] Lu | When is multi-version checkpointing needed?[END_REF]. This multiple-checkpoint approach has three major drawbacks. First, it is demanding in terms of storage: each checkpoint typically represents a copy of the entire memory footprint of the application, which may well correspond to several terabytes. The second drawback is the possibility of fatal failures. Indeed, if we keep k checkpoints in memory, the approach requires that the last checkpoint still kept in memory to precede the instant when the error currently detected struck. Otherwise, all live checkpoints would be corrupted, and one would have to re-execute the entire application from scratch. The probability of a fatal failure for various error distribution laws and values of k can be evaluated [START_REF] Aupy | On the combination of silent error detection and checkpointing[END_REF]. The third and most serious drawback of this approach applies even without memory constraints, i.e., if we could store an infinite number of checkpoints in memory. The critical point is to determine which checkpoint is the last valid one, information which is necessary to recover from a valid application state. However, because of the detection latency (which is unknown), we do not know when the silent error has indeed occurred, hence we cannot identify the last valid checkpoint, unless some verification mechanism is enforced.

We introduce such verification mechanisms in this chapter. In Section 2, we discuss several approaches to validation (recomputation, checksums, coherence tests, orthogonalization checks, etc). Then in Section 3 we adopt a general-purpose approach, which is agnostic of the nature of the verification mechanism. We consider a divisible-load application (which means that we can take checkpoints at any instant), and we partition the execution into computational patterns that repeat over time. The simplest pattern is represented by a work chunk followed by a verified checkpoint, which corresponds to performing a verification just before taking each checkpoint. If the verification succeeds, then one can safely store the checkpoint. If the verification fails, then a silent error has struck since the last checkpoint, and one can safely recover from it to resume the execution of the application. We compute the optimal length of the work chunk in the simplest pattern in Section 3.1, which amounts to revisiting Young and Daly's formula [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF][START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF] for silent errors. While taking a checkpoint without verification seems a bad idea (because of the memory cost, and of the risk of saving corrupted data), a validation step not immediately followed by a checkpoint may be interesting. Indeed, if silent errors are frequent enough, verifying the data in between two (verified) checkpoints, will reduce in expectation the detection latency and thus the amount of work to be re-executed due to possible silent errors. The major goal of Section 3 is to determine the best pattern composed of m work chunks, where each chunk is followed by a verification and the last chunk is followed by a verified checkpoint. We show how to determine m and the length of each chunk so as to minimize the makespan, that is the total execution time.

Then we move to application workflows. In Section 4, we consider application workflows that consist of a number of parallel tasks that execute on a platform, and that exchange data at the end of their execution. In other words, the task graph is a linear chain, and each task (except maybe the first and the last one) reads data from its predecessor and produces data for its successor. This scenario corresponds to a high-performance computing application whose workflow is partitioned into a succession of (typically large) tightly-coupled computational kernels, each of them being identified as a task by the model. At the end of each task, we can either perform a verification on the task output, or perform a verification followed by a checkpoint. We provide dynamic programming algorithms to determine the optimal locations of checkpoints and verifications.

The last technique that we illustrate is application-specific. In Section 5, we deal with sparse linear algebra kernels, and we show how to combine ABFT (Algorithm Based Fault Tolerance) with checkpointing. In a nutshell, ABFT consists in adding checksums to application data, and to view them as extended data items. The application performs the same computational updates on the original data and on the checksums, thereby avoiding the need to recompute the checksums after each update. The salient feature of this approach is forward recovery: ABFT is used both as an error verification and error correction mechanism: whenever a single error strikes, it can be corrected via ABFT and there is no need to rollback for recovery. Finally, we outline main conclusions and directions for future work in Section 6.

Verification mechanisms

Considerable efforts have been directed at error-checking to reveal silent errors. Error detection is usually very costly. Hardware mechanisms, such as ECC (Error Correcting Code) memory, can detect and even correct a fraction of errors, but in practice they are complemented with software techniques. General-purpose techniques are based on replication [START_REF] Engelmann | The case for modular redundancy in large-scale highh performance computing systems[END_REF][START_REF] Ferreira | Evaluating the Viability of Process Replication Reliability for Exascale Systems[END_REF][START_REF] Schroeder | Understanding failures in petascale computers[END_REF][START_REF] Zheng | Reliability-aware scalability models for high performance computing[END_REF]. Indeed, performing the operation twice and comparing the results of the replicas makes it possible to detect a single silent error. With Triple Modular Redundancy [START_REF] Lyons | The use of triple-modular redundancy to improve computer reliability[END_REF] (TMR) , errors can also be corrected by means of a voting scheme. Another approach, proposed by Moody et al. [START_REF] Ni | ACR: Automatic checkpoint/restart for soft and hard error protection[END_REF], is based on checkpointing and replication and enables detection and fast recovery of applications from both silent errors and hard errors.

Coming back to verification mechanisms, application-specific information can be helpful in designing ad hoc solutions, which can 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], but also ABFT techniques [START_REF] Bosilca | Algorithm-based fault tolerance applied to high performance computing[END_REF][START_REF] Huang | Algorithm-based fault tolerance for matrix operations[END_REF][START_REF] Shantharam | Fault tolerant preconditioned conjugate gradient for sparse linear system solution[END_REF], such as coding for the SpMxV (Sparse Matrix-Vector multiplication) kernel [START_REF] Shantharam | Fault tolerant preconditioned conjugate gradient for sparse linear system solution[END_REF], and coupling a higher-order with a lower-order scheme for Ordinary Differential Equations [START_REF] Benson | Silent error detection in numerical time-stepping schemes[END_REF]. These methods can only detect an error but not correct it. Self-stabilizing corrections after error detection in the conjugate gradient method are investigated by Sao and Vuduc [START_REF] Sao | Self-stabilizing iterative solvers[END_REF]. Also, Heroux and Hoemmen [START_REF] Heroux | Fault-tolerant iterative methods via selective reliability[END_REF] design a fault-tolerant GMRES algorithm capable of converging despite silent errors, and Bronevetsky and de Supinski [START_REF] Bronevetsky | Soft error vulnerability of iterative linear algebra methods[END_REF] provide a comparative study of detection cost for iterative methods. 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 nice instantiation of the checkpoint and verification mechanism that we study in this chapter 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. Consider a simple method such as the Preconditioned Conjugate Gradient (PCG) method: Chen's approach performs a periodic verification every d iterations, and a periodic checkpoint every d × c iterations, which is a particular case, with equispaced validations, of the approach presented later in Section 3.2. For PCG, the verification amounts to checking the orthogonality of two vectors and to recomputing and checking the residual. The cost of the verification is small if compared to the cost of an iteration, especially when the preconditioner requires many more flops than a SpMxV. As already mentioned, the approach presented in Section 3 is agnostic of the underlying error-detection technique and takes the cost of verification as an input parameter to the model.

Patterns for divisible load applications

In this section we explain how to derive the optimal pattern of interleaving checkpoints and verifications. An extended presentation of the results is available in [START_REF] Bautista-Gomez | Which verification for soft error detection[END_REF][START_REF] Benoit | Optimal resilience patterns to cope with fail-stop and silent errors[END_REF][START_REF] Cavelan | Assessing the impact of partial verifications against silent data corruptions[END_REF].

Revisiting Young and Daly's formula

Consider a divisible-load application, i.e., a (parallel) job that can be interrupted at any time for checkpointing, for a nominal cost C. To deal with fail-stop failures, the execution is partitioned into same-size chunks followed by a checkpoint, and there exist well-known formulae 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] to determine the optimal checkpointing period.

To deal with silent errors, the simplest protocol (see Figure 2) would be to perform a verification (at a cost V ) just before taking each checkpoint. If the verification succeeds, then one can safely store the checkpoint and mark it as valid. If the verification fails, then an error has struck since the last checkpoint, which is correct having been verified, and one can safely recover (which takes a time R) from that checkpoint to resume the execution of the application. This protocol with verifications zeroes out the risk of fatal errors that would force to restart the execution from scratch.

Time V C W V C W V C (Without error) Time V C W V R W V C W V C
Error Detection (With error)

Fig. 2 The simplest pattern: a work chunk W followed by a verification V and a checkpoint C.

To compute the optimal length of the work chunk W * , we first have to define the objective function. The aim is to find a pattern P (with a work chunk of length W followed by a verification of length V and a checkpoint of length C) that minimizes the expected execution time of the application. Let W base denote the base execution time of an application without any overhead due to resilience techniques (without loss of generality, we assume unit-speed execution). The execution is divided into periodic patterns, as shown in Figure 2. Let E(P) be the expected execution time of the pattern. For large jobs, the expected makespan W final of the application when taking failures into account can then be approximated by

W final ≈ E(P) W ×W base = W base + H(P) •W base
where

H(P) = E(P) W -1
is the expected overhead of the pattern. Thus, minimizing the expected makespan is equivalent to minimizing the pattern overhead H(P). Hence, we focus on minimizing the pattern overhead. We assume that silent errors are independent and follow a Poisson process with arrival rate λ . The probability of having at least a silent error during a computation of length w is given by p = 1e -λ w . We assume that errors cannot strike during recovery and verification. The following proposition shows the expected execution time of a pattern with a fixed work length W .

Proposition 1. The expected execution time of a pattern P with work length W is

E(P) = W +V +C + λW 2 + λW (V + R) + O(λ 2 W 3 ) . (1)
Proof. Let p = 1e -λW denote the probability of having at least one silent error in the pattern. The expected execution time obeys the recursive formula

E(P) = W +V + p(R + E(P)) + (1 -p)C . (2)
Equation ( 2) can be interpreted as follows: we always execute the work chunk and run the verification to detect silent errors, whose occurrence requires not only a recovery but also a re-execution of the whole pattern. Otherwise, if no silent error strikes, we can proceed with the checkpoint. Solving the recursion in Equation ( 2), we obtain

E(P) = e λW (W +V ) + e λW -1 R +C .
By approximating e λ x = 1 + λ x + λ 2 x 2 2 up to the second-order term, we can further simplify the expected execution time and obtain Equation (1).

The following theorem gives a first-order approximation to the optimal work length of a pattern. Theorem 1. A first-order approximation to the optimal work length W * is given by

W * = V +C λ . (3)
The optimal expected overhead is

H * (P) = 2 λ (V +C) + O(λ ) . (4)
Proof. From the result of Proposition 1, the expected overhead of the pattern can be computed as

H(P) = V +C W + λW + λ (V + R) + O(λ 2 W 2 ) . (5)
Assume that the MTBF of the platform µ = 1/λ is large if compared to the resilience parameters. Then consider the first two terms of H(P) in Equation ( 5): the overhead is minimal when the pattern has length W = Θ (λ -1/2 ), and in that case both terms are in the order of λ 1/2 , so that we have

H(P) = Θ (λ 1/2 ) + O(λ ).
Indeed, the last term O(λ ) becomes also negligible when compared to Θ (λ 1/2 ). Hence, the optimal pattern length W * can be obtained by balancing the first two terms in Equation ( 5), which gives Equation (3). Then, by substituting W * back into H(P), we get the optimal expected overhead in Equation ( 4).

We observe from Theorem 1 that the optimal work length W * of a pattern is in Θ λ -1/2 , and the optimal overhead H * (P) is in Θ (λ 1/2 ). This allows us to express the expected execution overhead of a pattern as

H(P) = o ef W + o rw W + O(λ )
, where o ef and o rw are two key parameters that characterize two different types of overheads in the execution, and they are defined below. Definition 1. For a given pattern, o ef denotes the error-free overhead due to the resilience operations (e.g., verification, checkpointing), and o rw denotes the reexecuted work overhead, in terms of the fraction of re-executed work due to errors.

In the simple pattern we analyze above, these two overheads are given by o ef = V +C and o rw = λ , respectively. The optimal pattern length and the optimal expected overhead can thus be expressed as

W * = o ef o rw , H * (P) = 2 √ o ef • o rw + O(λ ) .
We see that minimizing the expected execution overhead H(P) of a pattern becomes equivalent to minimizing the product o ef × o rw up to the dominating term. Intuitively, including more resilient operations reduces the re-executed work overhead but adversely increases the error-free overhead, and vice versa. This requires a resilience protocol that finds the optimal tradeoff between o ef and o rw . We make use of this observation in the next section to derive the optimal pattern in a more complicated protocol where patterns are allowed to include several chunks.

Optimal pattern

If the verification cost is small when compared to the checkpoint cost, there is room for optimization. Consider the pattern illustrated in Figure 3 with three verifications per checkpoint. There are three chunks of size w 1 , w 2 , and w 3 , each followed by a verification. Every third verification is followed by a checkpoint.

Time V C w 1 V w 2 V w 3 V C (Without error) Time V C w 1 V w 2 V R w 1 V w 2 V w 3 V C Error Detection
(With error) Fig. 3 Pattern with three verifications per checkpoint.

To understand the advantages of such a pattern, assume w 1 = w 2 = w 3 = W /3 for now, so that the total amount of work is the same as in the simplest pattern. As before, a single checkpoint needs to be kept in memory, and each error leads to re-executing the work since the last checkpoint. But detection occurs much more rapidly in the new pattern, because of the intermediate verifications. If the error strikes during the first of the three chunks, it is detected by the first verification, and only the first chunk is re-executed. Similarly, if the error strikes the execution of the second chunk (as illustrated in the figure), it is detected by the second verification, and the first two chunks are re-executed. The entire frame of work needs to be re-executed only if the error strikes during the third chunk. Under the first-order approximation as in the analysis of Theorem 1, the average amount of work to reexecute is (1 + 2 + 3)w/3 = 2w = 2W /3, that is, the re-executed work overhead becomes o rw = 2λ /3. On the contrary, in the first pattern of Figure 2, the amount of work to re-execute is always W , because the error is never detected before the end of the pattern. Hence, the second pattern leads to a 33% gain in the re-execution time. However, this comes at the price of three times as many verifications, that is, the error-free overhead becomes o ef = 3V +C. This overhead is paid in every error-free execution, and may be an overkill if the verification mechanism is too costly.

This example shows that finding the best trade-off between error-free overhead (what is paid due to the resilience method, when there is no failure during execution) and execution time (when errors strike) is not a trivial task. The optimization problem can be stated as follows: given the cost of checkpointing C, recovery R, and verification V , what is the optimal pattern to minimize the (expectation of the) execution time? A pattern is composed of several work chunks, each followed by a verification, and the last chunk is always followed by both a verification and a checkpoint. Let m denote the number of chunks in the pattern, and let w j denote the length of the j-th chunk for 1 ≤ j ≤ m. Let W = ∑ m j=1 w j . We define β j = w j /W be the relative length of the j-th chunk so that β j ≥ 0 and ∑ m j=1 β j = 1. We let

β β β = [β 1 , β 2 , . . . , β m ].
The goal is to determine the pattern work length W , the number of chunks m as well as the relative length vector β β β . Proposition 2. The expected execution time of the above pattern is

E(P) = W + mV +C + λβ β β T Aβ β β W 2 + O( √ λ ) , (6)
where A is an m × m matrix whose diagonal coefficients are equal to 1 and whose other coefficients are all equal to 1 2 .

Proof. Let p j = 1e -λ w j denote the probability of having at least one silent error in chunk j. To derive the expected execution time of the pattern, we need to know the probability q j that the chunk j actually gets executed in the current attempt. The first chunk is always executed, so we have q 1 = 1. Consider the second chunk, which is executed if no silent error strikes the first chunk, hence q 2 = 1p 1 . In general, the probability that the j-th chunk gets executed is

q j = j-1 ∏ k=1 (1 -p k ) .
Now, we are ready to compute the expected execution time of the pattern. The following gives the recursive expression:

E(P) = m ∏ k=1 (1 -p k ) C + 1 - m ∏ k=1 (1 -p k ) (R + E(P)) + m ∑ j=1 q j (w j +V ) . (7)
Specifically, line 1 of Equation [START_REF] Bosilca | Algorithm-based fault tolerance applied to high performance computing[END_REF] shows that the checkpoint at the end of the pattern is performed only when there has been no silent error in any of the chunks. Otherwise, we need to re-execute the pattern, after a recovery, as shown in line 2.

Finally, line 3 shows the condition for each chunk j to be executed. By simplifying Equation ( 7) and approximating the expression up to the second-order term, as in the proof of Proposition 1, we obtain

E(P) = W + mV +C + λ fW 2 + O( √ λ ) ,
where f = ∑ m j=1 β j ∑ m k= j β k , and it can be concisely written as f = β β β T Mβ β β , where M is the m × m matrix given by m i, j = 1 for i ≤ j 0 for i > j .

By replacing M by its symmetric part A = M+M T 2 , which does not affect the value of f , we obtain the matrix A whose diagonal coefficients are equal to 1 and whose other coefficients are all equal to 1 2 , and the expected execution time in Equation ( 6).

Theorem 2. The optimal pattern has m * equal-length chunks, total length W * and is such that:

W * = m * V +C 1 2 1 + 1 m * λ , (8)
β * j = 1 m * for 1 ≤ j ≤ m * , (9)
where m * is either max(1, m * ) or m * with

m * = C V . ( 10
)
The optimal expected overhead is

H * (P) = √ 2λC + √ 2λV + O(λ ) . (11)
Proof. Given the number of chunks m with ∑ m j=1 β j = 1, the function f = β β β T Aβ β β is shown to be minimized [10, Theorem 1 with r = 1] when β β β follows Equation ( 9), and its minimum value is given by f * = 1 2 1 + 1 m . We derive the two types of overheads as follows:

o ef = mV +C , o rw = 1 2 1 + 1 m λ .
The optimal work length W * = o ef o rw for any fixed m is thus given by Equation [START_REF] Bougeret | Checkpointing strategies for parallel jobs[END_REF]. The optimal number of chunks m * shown in Equation ( 10

Linear workflows

For an application composed of a chain of tasks, the problem of finding the optimal checkpoint strategy, i.e., of determining which tasks to checkpoint, in order to minimize the expected execution time when subject to fail-stop failures, has been solved by Toueg and Babaoglu [START_REF] Toueg | On the optimum checkpoint selection problem[END_REF], using a dynamic programming algorithm. We revisit the problem for silent errors by exploiting verification in addition to checkpoints. An extended presentation of the results is available in [START_REF] Benoit | Assessing general-purpose algorithms to cope with fail-stop and silent errors[END_REF][START_REF] Benoit | Two-level checkpointing and partial verifications for linear task graphs[END_REF].

Setup

To deal with silent errors, resilience is provided through the use of checkpointing coupled with an error detection (or verification) mechanism. When a silent error is detected, we roll back to the nearest checkpoint and recover from there. As in Section 3.1, let C denote the cost of checkpointing, R the cost of recovery, and V the cost of a verification.

We consider a chain of tasks T 1 , T 2 , . . . , T n , where each task T i has a weight w i corresponding to the computational load. For notational convenience, we also define W i, j = ∑ j k=i+1 w k to be the time to execute tasks T i+1 to T j for any i ≤ j. Once again we assume that silent errors occur following a Poisson process with arrival rate λ and that the probability of having at least one error during the execution of W i, j is given by p i, j = 1e -λW i, j .

We enforce that a verification is always taken immediately before each checkpoint, so that all checkpoints are valid, and hence only one checkpoint needs to be maintained at any time during the execution of the application. Furthermore, we assume that errors only strike the computations, while verifications, checkpoints, and recoveries are failure-free.

The goal is to find which task to verify and which task to checkpoint in order to minimize the expected execution time of the task chain. To solve this problem, we derive a two-level dynamic programming algorithm. For convenience, we add a virtual task T 0 , which is always checkpointed, and whose recovery cost is zero. This accounts for the fact that it is always possible to restart the application from scratch at no extra cost. In the following, we describe the general scheme when considering both verifications and checkpoints.

Dynamic programming

Figures 4 and5 illustrate the idea of the algorithm, which contains two dynamic programming levels, responsible for placing checkpoints and verifications, respectively, as well as an additional step to compute the expected execution time between two verifications. The following describes each step of the algorithm in detail.

T 0 V C T 1 . . . T c 1 V C T c 1 +1 . . . T c 2 V C . . . E ckpt (c 1 ) E veri f (c 1 , c 2 ) E ckpt (c 2 )
Fig. 4 First level of dynamic programming (E ckpt ).

Placing checkpoints. The first level focuses on the placement of verified checkpoints, i.e., checkpoints preceded immediately by a verification. Let E ckpt (c 2 ) denote the expected time to successfully execute all the tasks from T 1 to T c 2 , where T c 2 is verified and checkpointed. Now, to find the last verified checkpoint before T c 2 , we try all possible locations from T 0 to T c 2 -1 . For each location, say c 1 , we call the function recursively with E ckpt (c 1 ) (for placing checkpoints before T c 1 ), and compute the expected time to execute the tasks from T c 1 +1 to T c 2 . The latter is done through E veri f (c 1 , c 2 ), which also decides where to place additional verifications between T c 1 +1 and T c 2 . Finally, we add the checkpointing cost C (after T c 2 ) to E ckpt (c 2 ).

Overall, we can express E ckpt (c 2 ) as follows:

E ckpt (c 2 ) = min 0≤c 1 <c 2 {E ckpt (c 1 ) + E veri f (c 1 , c 2 ) +C} .
Note that a location c 1 = 0 means that no further checkpoints are added. In this case, we simply set E ckpt (0) = 0, which initializes the dynamic program. The total expected time to execute all the tasks from T 1 to T n is thus given by E ckpt (n).

. . .

T c 1 V C T c 1 +1 . . . T v 1 V T v 1 +1 . . . T v 2 V . . . E veri f (c 1 , v 1 ) E(c 1 , v 1 , v 2 ) E veri f (c 1 , v 2 )
Fig. 5 Second level of dynamic programming (E veri f ) and computation of expected execution time between two verifications (E).

Placing additional verifications. The second level decides where to insert additional verifications between two tasks with verified checkpoints. The function is initially called from the first level between two checkpointed tasks T c 1 and T c 2 , each of which also comes with a verification. Therefore, we define E veri f (c 1 , v 2 ) as the expected time to successfully execute all the tasks from T c 1 +1 to T v 2 , knowing that the last checkpoint is right after task T c 1 , and there is no additional checkpoint between T c 1 +1 and T v 2 . Note that E veri f (c 1 , v 2 ) accounts only for the time required to execute and verify these tasks. As before, we try all possible locations for the last verification between T c 1 and T v 2 and, for each location v 1 , we call the function recursively with E veri f (c 1 , v 1 ). Furthermore, we add the expected time needed to successfully execute the tasks T v 1 +1 to T v 2 , denoted by E(c 1 , v 1 , v 2 ), given the position c 1 of the last checkpoint. Overall, we can express E veri f (c 1 , v 2 ) as follows:

E veri f (c 1 , v 2 ) = min c 1 ≤v 1 <v 2 {E veri f (c 1 , v 1 ) + E(c 1 , v 1 , v 2 )} . (12)
Again, the case v 1 = c 1 means that no further verification is added, so we initialize the dynamic program with E veri f (c 1 , c 1 ) = 0. Note that the verification cost V at the end of task T v 2 will be accounted for in the function

E(c 1 , v 1 , v 2 ).
Computing expected execution time between two verifications. Finally, to compute the expected time to successfully execute several tasks between two verifications, we need the position of the last checkpoint c 1 , as well as the positions of the two verifications v 1 and v 2 . First, we pay W v 1 ,v 2 by executing all the tasks from T v 1 +1 to T v 2 , followed by the cost of verification V after T v 2 . During the execution, there is a probability p v 1 ,v 2 = 1e -λW v 1 ,v 2 of having a silent error, which will be detected by the verification after T v 2 . In this case, we need to perform a recovery from the last checkpoint after T c 1 with a cost R (set to 0 if c 1 = 0), and re-execute the tasks from there by calling the function E veri f (c 1 , v 1 ) followed by E(c 1 , v 1 , v 2 ). Therefore, we can express E(c 1 , v 1 , v 2 ) as follows:

E(c 1 , v 1 , v 2 ) = W v 1 ,v 2 +V + p v 1 ,v 2 R + E veri f (c 1 , v 1 ) + E(c 1 , v 1 , v 2 ) . ( 13
)
Simplifying Equation ( 13), we get

E(c 1 , v 1 , v 2 ) = e λW v 1 ,v 2 (W v 1 ,v 2 +V ) + e λW v 1 ,v 2 -1 R + E veri f (c 1 , v 1 ) .
Complexity. The complexity is dominated by the computation of the expected completion time table E veri f (c 1 , v 2 ), which contains O(n 2 ) entries, and each entry depends on at most n other entries that are already computed. All tables are computed in a bottom-up fashion, from the left to the right of the task chain. Hence, the overall complexity of the algorithm is O(n 3 ).

ABFT and checkpointing for linear algebra kernels

In this section we introduce ABFT (Algorithm Based Fault Tolerance) as an application-specific technique which allows for both error detection and correction. We streamline our discussion on the CG method, however, the techniques that we describe are applicable to any iterative solver that uses sparse matrix vector multiplies and vector operations. This list includes many of the non-stationary iterative solvers such as CGNE (Conjugate Gradient on Normal Equations), BiCG (Bi-Conjugate Gradient), BiCGstab (Bi-Conjugate Gradient Stabilized), where sparse matrix transpose vector multiply operations also take place. Preconditioned variants of these solvers with an approximate inverse preconditioner (applied as an SpMxV, or two SpMxVs) can also be made fault-tolerant with the proposed scheme. The extension to PCG is described in [START_REF] Fasi | A backward/forward recovery approach for the preconditioned conjugate gradient method[END_REF]. In Section 5.1, we first provide a background on the CG method and give an overview of 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 ABFT protection schemes. Then we detail ABFT techniques for the SpMxV kernel.

CG and fault tolerance mechanisms

The code for the CG method is shown in Algorithm 1. The main loop features a sparse matrix-vector multiply, two inner products (for p i q and r i+1

2 ), and three vector operations of the form axpy.

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.

We investigate three fault tolerance mechanisms. The first one is ONLINE-DETECTION; this is Chen's original approach modified to save the matrix A in addition to the current iteration vectors. This is needed when a silent error is de-Algorithm 1 The Conjugate Gradient algorithm for a positive definite matrix A.

Input: A ∈ R n×n , b, x 0 ∈ R n , ε ∈ R Output: x ∈ R n : Ax -b ≤ ε 1: r 0 ← b -Ax 0 ; 2: p 0 ← r 0 ; 3: i ← 0; 4: while r i > ε ( A • r 0 + b ) do 5: q ← Ap i ; 6:
α i ← r i 2 /p i q; 7:

x i+1 ← x i + α p i ; 8:

r i+1 ← r i -α q; 9:
β ← r i+1 2 / r i 2 ; 10:

p i+1 ← r i+1 + β p i ; 11:
i ← i + 1; 12: end while 13: return x i ; tected: if this error comes for a corruption in data memory, we need to recover with a valid copy of the data matrix A. The second one is ABFT-DETECTION, which detects errors and restarts from the most recent checkpoint. The thirds one is ABFT-CORRECTION, which detects errors and corrects if there was only one, otherwise restarts from the last checkpoint. The three methods under the study keep a valid copy of A and have exactly the same checkpoint cost.

We now introduce the ingredients of our own protection and verification mechanisms ABFT-DETECTION and ABFT-CORRECTION. We use ABFT techniques to protect the SpMxV, its result (hence the vector q), the matrix A and the input vector p i . As ABFT methods for vector operations is as costly as a repeated computation, we use TMR for them for simplicity. That is we do not protect p i , q, r i , and x i of the ith loop beyond the SpMxV at line 5 with ABFT, but we compute the dots, norms and axpy operations in resilient mode.

Although theoretically possible, constructing 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 errors. Therefore, we focus on detecting up to two errors and correcting single errors. 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 results out of three are correct (we assume errors do not strike the vector data here). Our fault-tolerant CG versions thus have the following ingredients: ABFT to detect up to two errors in the SpMxV and correct up to one; TMR for vector operations; and checkpoint and roll-back in case errors are not corrected. In the rest of this section, we discuss the proposed ABFT method for the SpMxV (combining ABFT with checkpointing is later in Section 5.3).

ABFT-SpMxV

The overhead of the standard single error correcting ABFT technique is too high for the sparse matrix-vector product case. Shantaram 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 detection of single errors striking either the computation or the memory representation of the two input operands (matrix and vector). As their results depend on the sparse storage format adopted, throughout this section we assume that sparse matrices are stored in the compressed storage format by rows (CSR) format [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]Sec. 3.4], that is by means of three distinct arrays, namely Colid ∈ N nnz(A) , Val ∈ R nnz(A) and Rowidx ∈ N n+1 .

Shantaram 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=1 a i, j (14)
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 the 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, that seems to restrict too much the practical applicability of their ABFT scheme. Shantaram et al. need this condition to ensure the 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. Let us 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 [START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF], 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 one 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, (15)
Algorithm 2 ABFT-protected SpMxV, detection of 2 errors, correction of 1 error

Input: A ∈ R n×n (asVal ∈ R nnz(A) , Colid ∈ N nnz(A) , Rowidx ∈ R n ), x ∈ R n Output: y = Ax,
d x = W y -C x; 22: d x = W (x -y) -M x; 23: if d r = 0 ∧ d x = 0 ∧ d x = 0 then 24:
return y; 25: for any i, if and only if it lies on h i . As the Lebesgue measure in R n of an hyperplane of R n itself is zero, the union of these hyperplanes is measurable with m n ( n i=1 h i ) = 0, where m n denotes the Lebesgue measure of R n . Therefore, the probability that a vector w randomly picked in R n does not satisfy condition [START_REF] Dongarra | The international exascale software project: a call to cooperative action by the global high-performance community[END_REF] 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 bigger 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 execution time and rounding errors. 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]Chapter 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 of its entries a constant value chosen so that all of the elements of the new vector are different from zero. We give the result in Theorem 3 for the simpler case of single error detection without correction, in which case Algorithm 2 has W = (1, . . . , 1) at line 1 and raises an error at line 26 (instead of correcting the error) if the tests at line 23 are not passed. The cases of multiple error detection and single error correction are proved in a technical report [START_REF] Fasi | Combining Algorithm-based Fault Tolerance and Checkpointing for Iterative Solvers[END_REF]Section 3.2].

Theorem 3 (Correctness of Algorithm 2 for error detection). 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, (16)
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

The proof of this theorem is technical and is available elsewhere [START_REF] Fasi | Combining Algorithm-based Fault Tolerance and Checkpointing for Iterative Solvers[END_REF]Theorem 1]. The function COMPUTECHECKSUM in Algorithm 2 requires just the knowledge of the matrix. Hence in the common scenario of many SpMxVs 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 the checksumming techniques.

Extensions to k ≥ 2 errors are discussed elsewhere [START_REF] Fasi | Combining Algorithm-based Fault Tolerance and Checkpointing for Iterative Solvers[END_REF]Section 3.2], where the following are detailed. The method just described can be extended to detect up to a total of k errors anywhere in the computation, in the representation of A, or in the vector x. Building up the necessary structures requires O (k nnz(A)) time, and the overhead per SpMxV is O(kn). For the particular case of k = 2 a result similar to that in Theorem 3 is also shown.

We now discuss error correction. If at least one of the tests at line 23 of Algorithm 2 fails, the algorithm invokes CORRECTERRORS in order to determine whether just one error struck either the computation or the memory and, in that case, to correct it. Indeed, whenever a single error is detected, disregarding its location (i.e., computation or memory), it can be corrected by means of a succession of various steps, as explained below; if need be, partial recomputations of the result are performed.

Specifically, we proceed as follows. To detect errors striking Rowidx, we compute the ratio d 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 it is so, the algorithm concludes that the d-th element of Rowidx is faulty 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 d has to be recomputed. • if z C > 2, then errors can be detected but not corrected, and an error is emitted.

To correct errors striking x, the algorithm computes d, 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 icx 1 and corrects x d = x dτ. The result is updated by subtracting from y the vector y τ = Ax τ , where x τ ∈ R n×n is such that x τ d = τ and x τ i = 0 otherwise. Finally, note that double errors could be shadowed when using Algorithm 2, but the probability of such an event is negligible. Still, there exists an improved version which avoids this issue by adding a third checksum [20, Section 3.2].

Performance model

The performance model is a simplified instance of the one discussed in Section 4, and we instantiate it for the three methods that we are considering, namely ONLINE-DETECTION, ABFT-DETECTION and ABFT-CORRECTION. We have a linear chain of identical tasks, where each task corresponds to one or several CG iterations. We execute T 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 the checkpoint, recovery and verification operations are error-free. For each method below, we let C, R and V be the respective cost of these operations. Finally, and as before, assume a Poisson process for errors and let q be the probability of successful execution for each chunk: q = e -λ T , where λ is the fault rate.

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 the following parameters:

• We have chunks of d iterations, hence T = dT iter , where T iter is the raw cost of a CG iteration without any resilience method.

• The verification time V is the cost of the operations described in Section 5.1.

• 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 cumulated rate of λ = λ a + λ m , and the success probability for a chunk is q = e -λ T . The optimal values of d and s can be computed by the same method as in Section 4.

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, V 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 iteration. Just as for ABFT-DETECTION, a chunk is a single iteration, and T = T iter , but V 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-DETEC-TION. The probability of exactly k errors in time T is (λ T ) k k! e -λ T [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 λ Te -λ T , so that q = e -λ T + λ Te -λ T .

Experiments

Comprehensive tests were performed and reported in the technical report [START_REF] Fasi | Combining Algorithm-based Fault Tolerance and Checkpointing for Iterative Solvers[END_REF]. The main observation is that ABFT-CORRECTION outperforms both ONLINE-DETEC-TION 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

Both fail-stop errors and silent data corruptions are major threats to executing HPC applications at scale. While many techniques have been advocated to deal with failstop errors, the lack of an efficient solution to handle silent errors is a real issue.

We have presented both a general-purpose solution and application-specific techniques to deal with silent data corruptions, with a focus on minimizing the overhead. For a divisible load application, we have extended the classical bound of Young/Daly to handle silent errors by combining checkpointing and verification mechanisms. For linear workflows, we have devised a polynomial-time dynamic programming algorithm that decides the optimal checkpointing and verification positions. Then, we have introduced ABFT as an application-specific technique to both detect and correct silent errors in iterative solvers that use sparse matrix vector multiplies and vector operations.

Our approach only addresses silent data corruptions. While several techniques have been developed to cope with either type of errors, few approaches are devoted to addressing both of them simultaneously. Hence, the next step is to extend our study to encompass both fail-stop and silent data corruptions in order to propose a comprehensive solution for executing applications on large scale platforms.

Fig. 1

 1 Fig. 1 Error and detection latency.

  ) is obtained by minimizing F(m) = o ef × o rw . The number of chunks in a pattern can only be a positive integer, so m * is either max(1, m * ) or m * , since F(m) is a convex function of m. Finally, substituting Equation (10) back into H * (P) = 2 √ o ef × o rw + O(λ ) gives rise to the optimal expected overhead as shown in Equation (11).

else 26 :

 26 CORRECTERRORS(Val, Colid, Rowidx, x, x , C, M, d r , d x , d x , c r , c x );

  , performs the correction by subtracting the first component of d r to Rowidx d , and recomputes y d and y d-1 , if the error in Rowidx d is a decrement; or y d+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 corrected together. Let now d be the ratio of the second component of d x to the first one. If d is near enough to an integer, the algorithm computes the checksum matrix 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 d , and can be corrected by simply recomputing this value. • If z C = 1, then the error concerns 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 d 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 d. 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 .

  M, c r , c x ] = COMPUTECHECKSUMS(Val, Colid, Rowidx); 5: return SPMXV(Val, Colid, Rowidx, x, x , M, c r , c x );

	correction of single error or detection of double error 1: global W ← [ 1 1 1 2 ... ... 1 n ] ∈ R 2×n ; 2: global W ← W 1 n+1 ∈ R 2×n+1 ;
	3: x ← x;
	4: [C, 6: function COMPUTECHECKSUMS(Val, Colid, Rowidx)
	7:	C ← W A;
	8:	M ← W -C;
	9:	c r ← W Rowidx;
	10:	c x ← W x;
	11:	return C, M, c r , c x ;
	12: function SPMXV(Val, Colid, Rowidx, x, x , C, M, c r , c x )
	13:	s r ← 0 ∈ R 2×1 ;
	14:	for i ← 1 to n do
	15: 16:	y i ← 0; s r ← s r +	w 1,i w 2,i Rowidx i ;
	17:	for j ← Rowidx i to Rowidx i+1 -1 do
	18:	ind ← Colid j ;
	19:	y i ← y i + Val j • x ind ;
	20:	d r = c r -s r ;
	21:

an auxiliary checksum vector c = [∑ n i=1 a i,1 + k, . . . , ∑ n i=1 a i,n + k],where k is such that ∑ n i=1 a i, j + k = 0 for 1 ≤ j ≤ n,

an auxiliary checksum y n+1 = k ∑ n i=i x i ,

an auxiliary counter s r initialized to 0 and updated at runtime by adding the value of the hit element each time the Rowidx array is accessed,

an auxiliary checksum c r = ∑ n i=1 Rowidx 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 , difference is in d x at line 21, ii. c x = ∑ n+1 i=1 y i , difference is in d x at line 22; iii. s r = c r , difference is in d r at line 20.