<figure><img src="image_1.png" /></figure>

HAL
open science

# Fault-tolerant numerical iterative algorithms at scale

Alix Tremodeux, Anne Benoit, Emmanuel Agullo, Thomas Herault, Luc Giraud, Yves Robert

<figure><img src="image_2.png" /></figure>

## To cite this version:

Alix Tremodeux, Anne Benoit, Emmanuel Agullo, Thomas Herault, Luc Giraud, et al., Fault-tolerant numerical iterative algorithms at scale. International Journal of High Performance Computing Applications, 2025. hal-05234063

HAL Id: hal-05234063

https://inria.hal.science/hal-05234063v1

Submitted on 1 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fault-tolerant numerical iterative algorithms at scale

Alix Tremodeux<fnref n="1" />, Anne Benoit<fnref n="1" />,<fnref n="2" />, Emmanuel Agullo<fnref n="3" />, Thomas Herault<fnref n="3" />,<fnref n="4" />, Luc Giraud<fnref n="3" /> and Yves Robert<fnref n="1" />

## Abstract

This work investigates how to protect numerical iterative algorithms from all types of errors that can strike at scale: fail-stop errors (a.k.a. failures) and silent errors, striking both as computation errors and memory bit-flips. We combine various techniques: detectors for computation errors, checksums for memory errors, and checkpoint/restart for failures. The objective is to minimize the expected time per iteration of the algorithm. We design a hierarchical pattern that combines and interleaves all these fault-tolerance mechanisms, and we determine the optimal periodic pattern that achieves this objective. We instantiate these results for the performance analysis of the Preconditioned Conjugate Gradient (PCG) algorithm: we report several scenarios where the optimal pattern dramatically decreases the overhead due to error mitigation.

## Keywords

Numerical iterative algorithms, computation errors, memory bit-flips, fail-stop errors, Preconditioned Conjugate Gradient.

## Introduction

Large-scale simulation is one of the main tools of modern science. To achieve new scientific challenges, simulations are processed at larger and larger scale, on platforms that comprise a multitude of computational units, and execute during long periods of time, up to several days. As a consequence, High Performance Computing (HPC) experiments are increasingly exposed to errors¹⁵.

Errors are usually classified in two main categories: fail-stop errors, that are immediately detected, and silent data corruptions, that can be detected through some verification mechanisms²⁶. Fail-stop errors correspond to permanent failures, e.g., processor crashes. Silent errors are disruptions that strike and stay undetected until they manifest eventually through strange application behavior. Silent errors arise from two main sources: computation errors and memory bit-flips. Protecting algorithms and software libraries from all these errors is a major concern within the HPC community. Although there have been many contributions in that direction (e.g. see the survey<fnref n="1" />), the design of a resilient holistic methodology to protect HPC applications from all error types is far from having been achieved. In particular, most studies are specifically restricted to only one type of errors, either fail-stop or silent.

This work aims at studying how to protect numerical iterative algorithms from all types of errors that can strike during execution, may they be either fail-stop or silent errors. We rely on system-level and numerical state-of-the-art techniques designed for each error type, and we show how they can be combined efficiently at scale. Our comprehensive approach mitigates the impact of all these errors and aims at guaranteeing a correct execution, regardless of the number and nature of errors that will strike. We combine various techniques: detectors for computation errors, checksums for

memory errors, and checkpoint/restart for fail-stop errors. We stress that our approach is agnostic of the type, nature and cost of all these techniques, as well as of the details of the target iterative algorithm. In essence, after each iteration of the algorithm, we decide whether or not it is worth to insert some fault-tolerant mechanisms, and if yes, which ones.

More precisely, we design and evaluate an execution pattern that repeats periodically, and within which each technique is applied with a specific frequency. In the scenarios that we explore for evaluation (see Section 13), computation errors strike more frequently than memory errors, which in turn strike more frequently than fail-stop errors. This motivates the following hierarchical shape for the periodic pattern:

* **Chunk**: after executing a certain number of iterations, we perform a verification for computation errors;

* **Segment**: after several chunks, we take a memory checkpoint to be able to recover from a memory error;

* **Pattern**: after several segments, we take a (costly) checkpoint on stable storage in order to be able to recover from fail-stop errors.

The size of the chunks, the number of chunks per segment, and the number of segments in the pattern are key parameters of our approach. All these parameters are to be determined according to the duration of each iteration of the algorithm,

## Corresponding author

Anne Benoit, Laboratoire LIP, ENS Lyon, 69364 Lyon Cedex 07, France.

Email: Anne.Benoit@ens-lyon.fr

<footnote n="1">¹ ENS Lyon, Laboratoire LIP & Inria, Lyon, France</footnote>

<footnote n="2">² Institut Universitaire de France</footnote>

<footnote n="3">³ Inria Centre at the University of Bordeaux, France</footnote>

<footnote n="4">⁴ University Tennessee Knoxville, TN, USA</footnote>
the frequency of each error type, and the cost of each fault-
tolerance mechanism.

The first main contribution of this work is to introduce
such a periodic execution pattern capable of coping with all
error sources that can strike the application simultaneously;
the second main contribution is to provide a performance
analysis that enables to characterize the optimal pattern as a
function of platform, application and resilience parameters.
The optimal pattern is the one that minimizes the expected
time per iteration, when both failure-free overhead (cost of
resilience mechanisms) and failure-induced overhead (time
lost by downtime, recovery and re-execution) are taken into
account. The interplay between fail-stop and silent errors
dramatically complicates the derivation of the expected
time per iteration, which characterizes the progress of the
application.

While the design and analysis methodology of the
execution pattern is agnostic of the target iterative algorithm,
we focus on the Preconditioned Conjugate Gradient (PCG)
algorithm for the practical evaluation. PCG is one of
the most widely used algorithms for the solution of
large problems that appear in many numerical simulations
(ocean, atmosphere³⁶, incompressible fluid flows¹⁹, plasma
physics¹², to name a few). In addition, PCG is used for
the TOP500 HPCG benchmark⁴⁹. We evaluate execution
patterns for several realistic scenarios of the deployment of
PCG at scale, in order to assess the impact of each parameter.
Our in-house toolkit is publicly available for the sake of
reproducibility of the results, and to enable further tests and
assessments with different execution scenarios.

We stress that this work provides neither actual
experiments on a large-scale platform, nor in-house Monte
Carlo simulations with error traces (or artificial error
injection methods). Instead, we make a toolkit available
for the community to study the impact of errors on the
performance of numerical iterative algorithms at scale. Our
study lays the foundations for deciding which resilience
mechanism to deploy, given its cost, the frequency of the
corresponding error type, and the detailed characteristics of
the target algorithm.

This paper is organized as follows. Section surveys
related work. Section outlines the generic framework, with
a short presentation of each error source, together with their
mitigation techniques. Section is devoted to the design and
theoretical analysis of the optimal pattern, which constitutes
the heart of our algorithmic contribution. Section details the
nature and cost of each fault-tolerant mechanism for PCG.
Section 13 presents the evaluation framework and reports the
results obtained from an extensive set of scenarios. Finally,
Section 13 provides concluding remarks and hints for future
work.

Related work

We survey related work on fault-tolerance in Section .
Section , discusses several papers that focus on our case
study, namely the PCG algorithm. We conclude with a brief
summary in Section .

Fault-tolerance techniques

Fault-tolerant approaches cover a wide spectrum of
techniques, depending upon the type of errors that are
addressed and mitigated. For fail-stop errors, checkpoint-
restart is de-facto standard strategy. After each failure,
the application executing on the faulty processor (and
likely on many other processors for a large parallel
application) is interrupted and must be restarted. Without
checkpointing, all the work executed for the application
is lost. With checkpointing, the execution can resume
from the last checkpoint, after some downtime (enroll
a spare to replace the faulty processor) and a recovery
(read the checkpoint). Several variants of this policy have
been studied; coordinated, hierarchical, multi-level, and in-
memory, to name a few, as described in an overview<sup>26</sup>.
A natural strategy is to checkpoint periodically, and one
must decide how often to checkpoint, i.e., derive the optimal
checkpointing period. An optimal strategy is defined as a
strategy that minimizes the expectation of the execution time
of the application. For a preemptible application, where one
can checkpoint at any time, the classical formula due to
Young<sup>50</sup> and Daly<sup>20</sup> states that the optimal checkpointing
period is &radic;2&mu;C, where &mu; is the application Mean Time
Between Failures (MTBF) and C the checkpoint cost.

While fail-stop errors lead to fatal interruptions (such as
processor crashes) and cause the loss of the entire memory
of the processor, silent errors, a.k.a. silent data corruptions,
only impact a given process and lead to incorrect results.
But a silent error strikes undetected and the processor can
continue its execution; sometimes the silent error can be
detected and corrected, and some other times it degenerates
into a fatal fail-stop error. We classify silent errors into
two categories: computation errors, which are transient and
caused by arithmetic errors in the Arithmetic and Logic
Unit (ALU), and memory errors, which are persistent and
correspond to bit flips in memory. These bit flips can strike in
the dynamic random-access memory (DRAM) due to cosmic
radiation, overheat and other sources, but also in the L1 cache
which is usually not well protected, or in the L2 cache which
might be protected by one parity bit<sup>40,41,51,52</sup>.

There are several software and hardware mechanisms to
detect and correct silent errors, such as parity bits, error
correcting codes (ECCs), and Chip-kill technology. They
have been implemented to protect the DRAM and different
cache layers to some extent. However, the closer the data
is to the processing unit, the more frequent the access to
that data and therefore the higher the overhead of these
methods. Thus, processor caches are not protected by ECC
in general, but by weaker mechanisms, like simple parity,
exposing a higher risk of undetectable error in case of
multiple simultaneous bit flips. Buses are also often a weak
link in the chain of protection, making all data transfers
at higher risk. In addition, the constant effort to reduce
operational voltage to save energy increases the likelihood
of silent errors. Although many silent errors caused by one
or multiple bits that spontaneously flip to the opposite state
are caught by the above-mentioned hardware mechanisms,
in reality, some bit flips still manage to pass undetected<sup>7,48</sup>.
In a nutshell, silent errors have become a major threat due
to the increase in problem size<sup>46</sup>: the larger the problem, the
more memory to be used to store the data, the more frequent
the errors, and the higher the probability of overriding ECC protection, generating multiple errors.

A major problem with silent errors is detection latency: contrarily to a fail-stop error whose detection is immediate, a silent error is identified only when the corrupted data is activated and/or leads to an unusual application behavior. Because checkpoint and rollback recovery assumes instantaneous error detection, this raises a new difficulty: if the error stroke before the last checkpoint, and is detected after that checkpoint, then the checkpoint is corrupted and cannot be used to restore the application. Some verification mechanisms, or detectors, must hence be enforced to guarantee that chekpoints are not corrupted. As stated above, hardware mechanisms, such as ECC memory, can detect and even correct a fraction of errors, but in practice, they are complemented with software techniques. The only general-purpose method is to replicate the execution of the target computational kernel on two sets of processors (i.e., duplication) and to compare the results of both executions. If they do not coincide, an error has been detected, and the application must be executed a third time. To avoid a-posteriori re-execution, triplication (originally known as triple modular redundancy and voting<sup>35</sup>) can be enforced, which allows for error correction in addition to error detection, using a simple majority vote.

Application-specific information can be useful to enable ad-hoc solutions, which dramatically decreases the cost of detection. Many techniques have been advocated. They include memory scrubbing<sup>30</sup> and Algorithm-Based Fault Tolerance (ABFT) techniques<sup>13,29,44</sup>, such as coding for the sparse-matrix vector multiplication kernel<sup>44</sup>, blockwise checksum calculation for error-bounded lossy compressor<sup>33</sup>, and coupling a higher-order with a lower-order scheme for Partial Differential Equations (PDE)<sup>11</sup>. These methods can only detect an error but do not correct it.

## Protecting the PCG algorithm

In this section, we focus on related work on the protection of the Preconditioned Conjugate Gradient (PCG) algorithm from silent errors. We refer to¹ for a more extensive survey on the resilience of numerical algorithms. The main idea from the literature is to design detection mechanisms based upon intrinsic properties of the PCG algorithm<sup>14,16,27</sup>. A first criterion is to check the orthogonality of some vectors; in finite precision arithmetic, the scalar product may not be computed exactly, hence a threshold has been introduced<sup>17</sup>. Another widely used criterion is to recompute the residual<sup>2,17</sup>. Other properties from the matrix² or the preconditioner³⁷ have been introduced to reduce the cost of the detection and improve its recall. Many of these works do not discuss correction techniques after having detected an error; checkpoint and rollback is used in<sup>2,17</sup>. The work²³ goes beyond and optimizes the ratio of the detection frequency and the checkpoint frequency in order to minimize execution time. All these techniques rely completely on the detection of error by the program. A completely different approach⁴³ advocates no rollbacks, but instead it aims at self-stabilizing the algorithm by recomputing the residual at regular intervals; however, this could impact the convergence rate in a non predictable way.

## Synthesis

To the best of our knowledge, this work is the first attempt to protect numerical iterative applications from all types of errors: fail-stop, computation and memory. Previous works focused either on an optimal pattern independent from a particular algorithm<sup>9,10</sup>, or they were limited to the detection and/or correction of transient errors in PCG<sup>2,17,23,37,43</sup>.

## Framework

This section details the fault-tolerance framework that we envision for numerical iterative algorithms. We start with error sources and mitigation techniques in Section. Error probability distributions are specified in Section. Then in Section, we describe the periodic hierarchical pattern that is the key building block of our holistic approach. All the notations introduced in this section are summarized in Table 1.

## Error sources and mitigation techniques

As discussed earlier, we aim at protecting numerical iterative algorithms from three error categories: fail-stop errors, and two types of silent errors, namely computation and memory. The techniques described in this general framework are agnostic to the specific target algorithm. We refer to Section for a case study with the Preconditioned Conjugate Gradient (PCG) algorithm.

For fail-stop errors, we use checkpoint-restart, with a global checkpoint that saves the entire footprint of the algorithm on stable storage. This checkpoint includes both static and dynamic data. Typically, the matrix and the preconditioner form the static data, while all the vectors updated at each iteration correspond to the dynamic data. We let $C_{fs}$ denote this global checkpoint, and we also use $C_{fs}$ for its associated cost. Similarly, $R_{fs}$ denote the recovery (cost) from a $C_{fs}$ checkpoint.

For memory errors (e.g., bit flips), we recompute checksums to verify that the memory has not been corrupted. The cost of this memory verification is $V_m$. If the verification is successful, we take a local in-memory checkpoint of cost $C_{cm}$. The $V_m$ verification detects any memory error, both for static or dynamic data, while the $C_{cm}$ checkpoint stores only the dynamic data (updated vectors). This is in contrast to the global checkpoint $C_{fs}$, whose cost is expected much higher than that of $C_{cm}$. We let $R_{cm}$ denote the recovery from a $C_{cm}$ checkpoint.

For computation errors (e.g., incorrect arithmetic operations), we use an application-specific computation verification $V_c$. In addition to detecting computation errors, the $V_c$ verification is likely to detect some (but not all) memory errors. We expect the cost of $V_c$ to be much smaller than the cost of $V_m$.

In summary, the $V_c$ verification detects computation errors, the $V_m$ verification detects memory errors, the $C_{cm}$ checkpoint stores only the dynamic data (vectors) in memory, and the $C_{fs}$ checkpoint stores all data (vectors, matrix and preconditioner) on stable storage. For simplicity, we make the following assumptions:

1. fail-stop errors can strike anytime except during $C_{fs}$ and $R_{fs}$;
2. silent errors can strike neither during verifications $V_c$ and $V_m$, nor during checkpoints $C_{cm}$ and $C_{fs}$, nor during recoveries $R_{cm}$ and $R_{fs}$;

3. verifications $V_c$ and $V_m$ have perfect recall and precision 1, which means that we have neither false negative (undetected errors) nor false positive (false alarms).

These three assumptions are commonly used in the literature²⁶. The third one is critical to the correctness of the approach: we always perform $V_c$ and $V_m$ verifications before taking a local $C_{cm}$ checkpoint, to ensure its validity. Similarly, we only take a global $C_{fs}$ checkpoint after a $V_c$, $V_m$, $C_{cm}$ sequence for the same reason.

We acknowledge that assuming perfect verifications may be unrealistic in some practical settings. However, in this work, we aim at avoiding a prohibitively costly approach based on full redundancy via replication, or even triplication. Using perfect application-specific detectors is the state-of-the-art technique to guarantee that checkpoints are not corrupted⁶,²².

## Probability distributions

The probability distribution of inter-arrival times for memory and fail-stop errors both obey an exponential (memoryless) distribution²⁶. Let $\lambda_{mem}$ and $\lambda_{fs}$ be the parameter for respectively memory errors and fail-stop errors. Hence, $e^{-\lambda_{mem}t}$ and $e^{-\lambda_{fs}t}$ are, respectively, the memory and fail-stop probabilities of successful execution during a duration $t$ (i.e., no error of the specified type strikes during time $t$). The computation errors follow a geometric distribution of parameter $f_{calc}$, which is the success probability per iteration (i.e., no computation error strikes during the iteration). We need a discrete probability law for computation errors because they are counted per iteration, but it remains very close to an exponential distribution of parameter $\lambda_{calc}$, with $f_{calc} = e^{-\lambda_{calc}I}$, and where $I$ is the duration of an iteration. All three error probability distributions are considered as mutually independent.

## Hierarchical pattern

We aim at creating a flexible and holistic pattern where verifications and checkpoints are interleaved after some iterations, and which is repeated periodically during execution. The frequency of each resilience mechanism within the pattern will depend upon many parameters: the probability of each error source, the cost of verifications and checkpoints, and the time spent per iteration.

The key hypothesis that drives the design of the pattern is that computation errors are more frequent than memory errors, which in turn are more frequent than fail-stop errors. This calls for a hierarchical pattern where computation verifications $V_c$ are inserted more frequently than memory verifications $V_m$ and local checkpoints $C_{cm}$, which in turn are inserted more frequently than global checkpoints $C_{fs}$.

The periodic pattern begins with $n_{vc} \ge 1$ iterations, each taking $I$ time units. After these $n_{vc}$ iterations, we take a $V_c$ verification. This constitutes a chunk, of duration $T_{calc} = n_{vc}I + V_c$. Then, we repeat $n_{cm} \ge 1$ chunks before executing a $V_m$ verification followed by a $C_{cm}$ checkpoint. This constitutes a segment, of duration

$T_{mem} + C_{cm}$, where $T_{mem} = n_{cm}(n_{vc}I + V_c) + V_m$. We repeat $n_{fs} \ge 1$ segments before terminating the pattern with a $C_{fs}$ checkpoint. As illustrated in Figure 1, this describes the full pattern, of total length $T_{fs} + C_{fs}$, where $T_{fs} = n_{fs}(n_{cm}(n_{vc}I + V_c) + V_m + C_{cm})$. There are $n_{vc} \times n_{cm} \times n_{fs}$ iterations executed within the pattern.

## Objective function

The objective is to minimize the slowdown, which is defined as the ratio of the expected time needed per iteration over the time per iteration in a failure-free execution, with no fault-tolerance mechanism. Given a pattern with parameters $(n_{vc}, n_{cm}, n_{fs})$, let $E(n_{vc}, n_{cm}, n_{fs})$ be the expected time to execute the pattern successfully, from its beginning up to the last $C_{fs}$ checkpoint. We aim at minimizing the slowdown $S(n_{vc}, n_{cm}, n_{fs})$, where

$$S(n_{vc}, n_{cm}, n_{fs}) = \frac{E(n_{vc}, n_{cm}, n_{fs})}{n_{vc}n_{cm}n_{fs}I}.$$

We show how to derive a closed-form analytical formula for $E(n_{vc}, n_{cm}, n_{fs})$ in the next section.

## Optimal pattern

In this section, we first detail the derivation of $E(n_{vc}, n_{cm}, n_{fs})$ in Section. Then, we explain how to derive the optimal pattern in Section.

### Formula for $E(n_{vc}, n_{cm}, n_{fs})$

For simplicity, we write $E$ for $E(n_{vc}, n_{cm}, n_{fs})$ in this section. The main difficulty to compute $E$ is that we cannot rely on previous approaches when combining all different error sources. For fail-stop errors, we could try and dimension the pattern so that the time to execute the useful work $W$ would correspond to the Young/Daly period²⁰,⁵⁰; however, whenever a silent error strikes within the period, we need to rollback and re-execute some work (depending upon whether it is a computational or memory error), thereby leading to an actual execution time that may well largely exceed the failure-free work $W$. Furthermore, different error types induce different overheads in their recovery process.

An important observation is that the expected time to execute a given segment depends upon the position of that segment within the pattern. Therefore, we break down the expected time $E$ for the whole pattern into

$$E = \sum_{k=1}^{n_{fs}} E_k + C_{fs},$$

where $E_k$ is the mean time to complete segment $k$ (i.e., to successfully execute its memory checkpoint). If a silent error is detected within segment $k$, then only this segment is restarted, while if a fail-stop error is detected, then all previous $k-1$ segments have to be completed again before the current segment $k$ is restarted. Now, the formula for $E_k$ is obtained by a case analysis, depending upon whether there is no error, or which error source is the first to be detected during execution. We start with some notations:

* $P_{succ}^{fs} = e^{-\lambda_{fs}(T_{mem}+C_{cm})}$ is the probability of no fail-stop error within the segment;
<figure><img src="image_4.png" /><figcaption>Figure 1. Example of a pattern with $n_{cm} = 3$ and $n_{fs} = 2$.</figcaption></figure>

<table><caption>Table 1. Table of notations.</caption><thead><tr><th>Symbol</th><th>Definition</th><th>Symbol</th><th>Definition</th></tr></thead><tbody><tr><td>I</td><td>One iteration</td><td>$n_{vc}$</td><td>Number of iterations in a chunk</td></tr><tr><td>V<sub>c</sub></td><td>Computation verification</td><td>$n_{cm}$</td><td>Number of chunks in a segment</td></tr><tr><td>V<sub>m</sub></td><td>Memory verification</td><td>$n_{fs}$</td><td>Number of segments in a pattern</td></tr><tr><td>C<sub>cm</sub></td><td>Local memory checkpoint (vectors of the algorithm)</td><td>$T_{calc}$</td><td>= $n_{vc}I + V_c$</td></tr><tr><td>C<sub>fs</sub></td><td>Checkpoint (in case of a fail-stop error)</td><td>$T_{mem}$</td><td>= $n_{cm}(n_{vc}I + V_c) + V_m$</td></tr><tr><td>R<sub>cm</sub></td><td>Local recovery process from C<sub>cm</sub> checkpoint</td><td>$T_{fs}$</td><td>= $n_{fs}(n_{cm}(n_{vc}I + V_c) + V_m + C_{cm})$</td></tr><tr><td>R<sub>fs</sub></td><td>Recovery process from C<sub>fs</sub> checkpoint</td><td>$W$</td><td>= $n_{vc}n_{cm}n_{fs}I$</td></tr><tr><td>P<sup>z</sup><sub>succ</sub></td><td>Probability of no error of type z</td><td>$f_{calc}$</td><td>Parameter of geometric distrib. for computation errors</td></tr><tr><td>P<sup>z</sup><sub>fail</sub></td><td>= $1 - P<sup>z</sup><sub>succ</sub></td><td>$\lambda_{mem}$</td><td>Parameter of exponential distrib. for memory errors</td></tr><tr><td>P<sup>calc</sup><sub>fail</sub>(i)</td><td>Proba. that first computation error is in chunk i of segment</td><td>$\lambda_{fs}$</td><td>Parameter of exponential distrib. for fail-stop errors</td></tr><tr><td>P<sub>no.fs</sub></td><td>1 - P<sub>no.fs</sub>: the first error to strike is a fail-stop error</td><td></td><td></td></tr></tbody></table>

* $P_{succ}^{mem} = e^{-\lambda_{fs}T_{mem}}$ is the probability of no memory error within the segment;

* $P_{succ}^{calc}(i) = f_{calc}^{n_{vc}i}$ is the probability of no computation error during the first i chunks (hence during a total of $n_{vc}i$ iterations), and $P_{fail}^{calc}(i) = P_{succ}^{calc}(i-1)(1-f_{calc}^{n_{vc}})$ is the probability that the first computation error strikes during chunk i;

* Finally, $1 - P_{no.fs}$ is the probability that the first error in the segment is a fail-stop error. Its complicated value is determined below.

We can now write down the overall formula for $E_k$:

$$ 
\begin{aligned}
E_k = & P_{succ}^{fs} P_{succ}^{mem} P_{succ}^{calc}(n_{cm})[T_{mem} + C_{cm}] \quad (\text{Case 1}) \\
& + (1 - P_{succ}^{mem}) P_{succ}^{fs} P_{succ}^{calc}(n_{cm})[T_{mem} + R_{cm} + E_k] \quad (\text{Case 2})
\end{aligned}
$$

$$ + \sum_{i=1}^{n_{cm}} e^{-\lambda_{fs}iT_{calc}} P_{fail}^{calc}(i)[iT_{calc} + R_{cm} + E_k] \quad (\text{Case 3}) $$

$$ + (1 - P_{no.fs})[E_{lost} + R_{fs} + \sum_{i=0}^{k} E_i] \quad (\text{Case 4}) $$

Case 1 is when no error struck during the computation of the whole segment. It is weighted by the probability of success for all error types. Note that $P_{succ}^{calc}(n_{cm})$ is the probability of having no computation error throughout the segment.

Case 2 is when the first detected error in the segment is a memory error. This implies that no fail-stop nor computation error took place before the end of the memory verification at the end of the segment before the local in-memory checkpoint, hence the weight by the corresponding probabilities.

Case 3 is when the first detected error is a computation error. This error type can strike only during iterations. The

lost time induced by this error depends upon the chunk where it struck. If it is the i-th chunk of the segment, the lost time includes all the first i chunks, hence a time $iT_{calc}$. Note that in this situation, no fail-stop error has struck before the computation error, while a memory error may happen anytime before the computation error and it will not be detected. Again, we weight this case with its probability to happen.

Case 4 is when the first error to strike is a fail-stop error, which is immediately detected. When such an error strikes, the re-execution starts at the very beginning of the pattern, so all the previously completed segments have to be recomputed. Here, $E_{lost}$ is the expected lost time during the segment, whose value is

$$ 
\begin{aligned}
E_{lost} &= \int_{0}^{\infty} xP(X = x|X < T_{mem} + C_{cm})dx \\
&= \frac{1}{1 - e^{-\lambda_{fs}(T_{mem} + C_{cm})}} \int_{0}^{T_{mem} + C_{cm}} x\lambda_{fs}e^{-\lambda_{fs}x} dx \\
&= \frac{1}{\lambda_{fs}} - \frac{T_{mem} + C_{cm}}{e^{\lambda_{fs}(T_{mem} + C_{cm})} - 1},
\end{aligned}
$$

where X is an exponential random variable of parameter $\lambda_{fs}$ (see$^{26}$ for details).

There remains to evaluate $P_{no.fs}$, the probability that the first error to be spotted is not a fail-stop error. But by definition, the value of $P_{no.fs}$ is the sum of the probability
weights of the first three cases, hence

$$
\begin{align*}
P_{no,fs} &= P_{succ}^{fs} P_{succ}^{mem} P_{succ}^{calc}(n_{cm}) \\
&\quad + (1 - P_{succ}^{mem}) P_{succ}^{fs} P_{succ}^{calc}(n_{cm}) \\
&\quad + \sum_{i=1}^{n_{cm}} e^{-\lambda_{fs} i T_{calc}} P_{fail}^{calc}(i) \\
&= e^{-T_{mem}(\lambda_{fs} + \lambda_{mem})} - C_{cm} \lambda_{fs} (f_{calc})^{n_{vc} n_{cm}} \\
&\quad + (1 - e^{-\lambda_{mem} T_{mem}}) e^{-\lambda_{fs} T_{mem}} (f_{calc})^{n_{vc} n_{cm}} \\
&\quad + (1 - f_{calc}^{n_{vc}}) e^{-\lambda_{fs} T_{calc}} \left( \frac{(f_{calc}^{n_{vc}} e^{-\lambda_{fs} T_{calc}})^{n_{cm}} - 1}{f_{calc}^{n_{vc}} e^{-\lambda_{fs} T_{calc}} - 1} \right).
\end{align*}
$$

Altogether, after some re-ordering, and double-checking
the derivation with a computer algebra system, we are led to
the following recurrence equation:

$$
E_k = M + (1 - P_{no,fs}) \sum_{i=1}^{k-1} E_i \\
\qquad + (1 - P_{succ}^{fs} P_{succ}^{mem} P_{succ}^{calc}(n_{cm})) E_k,
$$

whose complexity is hidden into the value of M, which
is the part of the formula that does not depend upon k:

$$
\begin{align*}
M &= e^{-T_{mem}(\lambda_{fs} + \lambda_{mem})} - C_{cm}\lambda_{fs}(f_{calc})^{n_{vc}n_{cm}}[T_{mem} + C_{cm}] \\
  &\quad + (1 - e^{-\lambda_{mem}T_{mem}})e^{-\lambda_{fs}T_{mem}}(f_{calc})^{n_{vc}n_{cm}}[n_{cm}(n_{vc}I + V_c) + V_m + R_{cm}] \\
  &\quad + (1 - f_{calc}^{n_{vc}})e^{-\lambda_{fs}T_{calc}}\left(\frac{(f_{calc}^{n_{vc}}e^{-\lambda_{fs}T_{calc}})^{n_{cm}} - 1}{f_{calc}^{n_{vc}}e^{-\lambda_{fs}T_{calc}} - 1}\right)[R_{cm}] \\
  &\quad + (1 - f_{calc}^{n_{vc}})e^{-\lambda_{fs}T_{calc}}\left(\frac{(f_{calc}^{n_{vc}}e^{-\lambda_{fs}T_{calc}})^{n_{cm}}(n_{cm}(f_{calc}^{n_{vc}}e^{-\lambda_{fs}T_{calc}} - 1) - 1) + 1}{(f_{calc}^{n_{vc}}e^{-\lambda_{fs}T_{calc}} - 1)^2}\right)T_{calc} \\
  &\quad + \left[1 - \left[e^{-(n_{cm}(n_{vc}I + V_c) + V_m)(\lambda_{fs} + \lambda_{mem})} - C_{cm}\lambda_{fs}(f_{calc})^{n_{vc}n_{cm}}\right.\right. \\
  &\qquad \left. + (1 - e^{-\lambda_{mem}T_{mem}})e^{-\lambda_{fs}T_{mem}}(f_{calc})^{n_{vc}n_{cm}}\right. \\
  &\qquad \left. + (1 - f_{calc}^{n_{vc}})e^{-\lambda_{fs}T_{calc}}\left(\frac{(f_{calc}^{n_{vc}}e^{-\lambda_{fs}T_{calc}})^{n_{cm}} - 1}{f_{calc}^{n_{vc}}e^{-\lambda_{fs}T_{calc}} - 1}\right)\right]\right] \\
  &\quad \times \left[\frac{1}{\lambda_{fs}} - \frac{n_{cm}(n_{vc}I + V_c) + V_m + C_{cm}}{e^{\lambda_{fs}(T_{mem} + C_{cm})} - 1} + R_{fs}\right].
\end{align*}
$$

By induction on *k*, we derive a compact form for $E_k$:

**Lemma 1.**

$$
E_k = \frac{M}{P_{succ}^{fs} P_{succ}^{mem} P_{succ}^{calc}(n_{cm})} \left( \frac{1 - P_{no,fs}}{P_{succ}^{fs} P_{succ}^{mem} P_{succ}^{calc}(n_{cm})} + 1 \right)^k \quad (1)
$$

**Proof.** The proof of Equation (1) is done by induction on *k*. Let us denote $P_{succ}^{fs,mem,calc} = P_{succ}^{fs} P_{succ}^{mem} P_{succ}^{calc}(n_{cm})$. Initially,

$$
\begin{align*}
E_0 &= \frac{M}{P_{succ}^{fs,mem,calc}} \\
    &= \frac{M}{P_{succ}^{fs,mem,calc}} \left( \frac{1 - P_{no,fs}}{P_{succ}^{fs,mem,calc}} + 1 \right)^0.
\end{align*}
$$

Now, assume that the formula is true for all
$i \in 0, 1, \dots, k-1$. Then,

$$
\begin{align*}
E_k &= \frac{M}{P_{succ}^{fs,mem,calc}} \\
&\quad + \frac{(1 - P_{no,fs})}{P_{succ}^{fs,mem,calc}} \sum_{i=0}^{k-1} \frac{M}{P_{succ}^{fs,mem,calc}} \left( \frac{1 - P_{no,fs}}{P_{succ}^{fs,mem,calc}} + 1 \right)^i \\
&= \frac{M}{P_{succ}^{fs,mem,calc}} \\
&\quad \times \left( 1 + \frac{(1 - P_{no,fs})}{P_{succ}^{fs,mem,calc}} \frac{\left( \frac{1 - P_{no,fs}}{P_{succ}^{fs,mem,calc}} + 1 \right)^k - 1}{\frac{1 - P_{no,fs}}{P_{succ}^{fs,mem,calc}} + 1 - 1} \right) \\
&= \frac{M}{P_{succ}^{fs,mem,calc}} \left( \frac{1 - P_{no,fs}}{P_{succ}^{fs,mem,calc}} + 1 \right)^k,
\end{align*}
$$

which corresponds to Equation (1).

<figure><img src="image_5.png" /></figure>

Equation (1) is the last step before the final derivation of
$E = \sum_{k=1}^{n_{fs}} E_k + C_{fs}$. We finally obtain:

$$
E = \frac{M}{1 - P_{no,fs}} \left[ \left( \frac{1 - P_{no,fs}}{P_{succ}^{fs} P_{succ}^{mem} P_{succ}^{calc}(n_{cm})} + 1 \right)^{n_{fs}} - 1 \right] + C_{fs}.
$$

Optimal pattern

The formula for *E* is finally determined! We aim at
minimizing the slowdown *S* = *E* / (*n*<sub>vc</sub> *n*<sub>cm</sub> *n*<sub>fs</sub> I) over all possible
pattern parameters *n*<sub>vc</sub>, *n*<sub>cm</sub> and *n*<sub>fs</sub>. Given the complexity
of the formula, it is impossible to optimize *S* analytically.
However, when given all the other platform/application
parameters (the duration of an iteration, the distribution
of the three error types and the cost of all resilience
mechanisms), it is easy to evaluate *S* numerically and to
determine its minimum via an exhaustive search over the
three pattern parameters. It takes only a few milliseconds
to sweep over 10<sup>7</sup> values, letting (conservatively) *n*<sub>vc</sub> range
from 1 to 1000, and *n*<sub>cm</sub> and *n*<sub>fs</sub> range from 1 to 100.

It is then interesting to compare the expected cost of
the optimal pattern over that of the naive approach where
both verifications and checkpoints are taken after each
iteration, and which corresponds to the simplest pattern with
$n_{vc} = n_{cm} = n_{fs} = 1$. This naive approach is often the first
approach used to design a fault-tolerant algorithm³⁷.

**Case study: PCG**

This section focuses on the Preconditioned Conjugate
Gradient (PCG) algorithm as a case-study. First, we give
a quick background on the algorithm in Section . Then,
application-specific detection and correction techniques for
PCG are presented in Section 13.

PCG algorithm

The PCG algorithm is the iterative method of choice for
solving a linear system Ax = b, with A a real symmetric
positive definite matrix, and b and x real vectors. At
iteration i, the algorithm uses the computed residual ri, the
direction pi and a coefficient αi to move xi closer to the
solution x of the system²⁴. The preconditioner matrix M
aims at making convergence faster. $z_i$ is an additional vector related to the use of the preconditioner. The PCG algorithm is depicted in Algorithm 1. $\gamma_i$, $\alpha_i$, $\beta_i$ are scalars. The initialization of the algorithm is on lines [1,4]. The initial guess $x_0$ is chosen arbitrarily. The first line computes the first residual $r_0$. Then, the iteration is on lines [6,13]. At each iteration, the algorithm adds the approximate part of the solution depending on the direction given by the previous iteration, to $x_{i+1}$. Finally, the new direction $p_{i+1}$ is estimated and conjugated to all previous directions.

## Error mitigation techniques for PCG

In exact arithmetic, the PCG algorithm is a sophisticated and elegant numerical scheme that has many properties induced by the symmetric positive definiteness of the matrix $A$. Unfortunately, most characteristic properties are no longer valid in finite precision calculation, such as orthogonality among the residuals $r_i$ (i.e., $r_i^T r_j = \delta_{i,j}$), A-conjugacy of the descent direction $p_i$ (i.e., $p_i^T A p_j = \delta_{i,j}$), or equality between the computed residual $r_i$ and the true residual $b - Ax_i$ (i.e., $r_i = b - Ax_i$). Unfortunately, it is known that these equalities are no longer valid in finite precision arithmetic. Consequently, they cannot be used to detect silent errors. It is argued and shown through extensive experiments² that a robust detection mechanism can be designed by combining a check on a (first) bound on the residual gap, together with a (second) bound on $\alpha$. The first bound, on the so-called residual gap, defined by $\|r_i - (b - Ax_i)\|$, is derived using classical rounding error analysis techniques. This computation is the core of the additional cost associated with numerical verification, requiring an additional matrix-vector product, vector addition (axpy), and norm computation. The second check, among the numerous relationships that exist between the quantities computed by PCG, is one that is possibly less prone to defection due to finite precision calculation. Indeed it is a bound and not a strict equality as for the orthogonality or A-orthogonality properties. It was already presented in the original paper on CG²⁸ Th. 5.5:

$$ \forall i \quad \frac{1}{\lambda_{\max}(A)} < \alpha_i < \frac{1}{\lambda_{\min}(A)}, \quad (2) $$

where $\lambda_{\max}(A)$ and $\lambda_{\min}(A)$ denote the largest and smallest eigenvalues of $A$, respectively. It was shown² that checking

<table><caption>**Algorithm 1: PCG algorithm**</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>$r_0 \leftarrow b - Ax_0$</td></tr><tr><td>2</td><td>$z_0 \leftarrow M^{-1}r_0$</td></tr><tr><td>3</td><td>$p_0 \leftarrow z_0$</td></tr><tr><td>4</td><td>$\gamma_0 \leftarrow r_0^T z_0$</td></tr><tr><td>5</td><td>for i = 0 ... until convergence do</td></tr><tr><td>6</td><td>$q_i \leftarrow A p_i$</td></tr><tr><td>7</td><td>$\alpha_i \leftarrow \gamma_i / p_i^T q_i$</td></tr><tr><td>8</td><td>$x_{i+1} \leftarrow x_i + \alpha_i p_i$</td></tr><tr><td>9</td><td>$r_{i+1} \leftarrow r_i - \alpha_i q_i$</td></tr><tr><td>10</td><td>$z_{i+1} \leftarrow M^{-1} r_{i+1}$</td></tr><tr><td>11</td><td>$\gamma_{i+1} \leftarrow r_{i+1}^T z_{i+1}$</td></tr><tr><td>12</td><td>$\beta_{i+1} \leftarrow \gamma_{i+1} / \gamma_i$</td></tr><tr><td>13</td><td>$p_{i+1} \leftarrow z_{i+1} + \beta_i p_i$</td></tr></tbody></table>

only the lower bound $\frac{1}{\lambda_{\max}(A)} < \alpha_i$ is sufficient to make the overall detection robust. This is important from a practical point of view as the calculation or the tight approximation of $\lambda_{\min}(A)$ is generally expensive while $\lambda_{\max}(A)$ can often be cheaply approximated, using for instance randomized techniques²⁵. Once the approximation of $\lambda_{\max}(A)$ has been obtained (in a preliminary setup phase), the check itself $\frac{1}{\lambda_{\max}(A)} < \alpha_i$ (triggered at iteration $i$) comes at cost zero as $\alpha_i$ is a by-product of PCG. In the following, we assume the numerical verification is ensured with such a check combining a verification on the residual gap bound (at the cost of a matrix-vector product, an axpy and a norm computation) and alpha verification (at cost zero). We furthermore assume that the memory checkpoint stores only the vectors, while the fail-stop checkpoint stores the matrix and the preconditioner in addition to the vectors.

## Evaluation

This section aims at evaluating the efficiency of the approach on a variety of realistic scenarios. We first recall the key input parameters in Section 13. Then, we discuss scenarios that instantiate these parameters for PCG in Section 13, with the detail on the memory complexity of the preconditioners. Section 13 deals with error sources. Finally, the results and analyses of the scenarios are provided in Section 13.

### Input parameters

The formula for the slowdown $S(n_{vc}, n_{cm}, n_{fs})$ needs the values of several parameters:

* Time $I$ spent per iteration;

* Cost of resilience parameters: $V_c$ and $V_m$ for computation and memory verifications, $C_{cm}$ and $R_{cm}$ for local in-memory checkpoint and recovery, and $C_{fs}$ and $R_{fs}$ for global checkpoint and recovery;

* Frequency of each error source: $f_{calc}$ for the geometric distribution of computation errors, and $\lambda_{mem}$ and $\lambda_{fs}$ for the exponential distributions of memory and fail-stop errors.

### Scenarios for PCG

We consider the solution of linear systems, $Ax = b$, arising from the discretization of 3D elliptic PDE, which appear in the modeling of many complex systems, such as for atmosphere or ocean simulation³⁶, incompressible fluid flows¹⁹, plasma physics¹², to name a few. For these large dimension systems, iterative solvers are the only possible solvers and PCG is the best suited one. To speed up its convergence, preconditioning techniques are used. For linear systems arising from 3D PDE discretization, multigrid approaches are probably the best option for sufficiently smooth problems because of their linear complexity of required computation and memory with respect to the problem size. Multigrid preconditioning is used in the TOP500 HPCG benchmark⁴⁹. In the case of stiff problems, alternative preconditioning techniques that may be employed include domain decomposition methods (DDM). It should be noted that these preconditioning techniques may require greater computational and memory resources to ensure
further numerical robustness, especially when using a coarse grid correction as we consider here$^{3,47}$.

We consider both these multigrid and DDM classes of preconditioners. We first derive their memory costs on an analytical problem. We then consider large-scale real-life experiments from the literature to derive realistic timing parameters. Finally, we conduct the evaluation based on these parameters.

Assuming a 7-point stencil discretization of the PDE using a uniform Cartesian grid in the unit cube of side $m$, the square matrix $A$ is of order $N = m^3$, sparse, and it contains $7N$ non-zero elements. Assuming a COOordinates (COO) storage, it is composed of three vectors: row indices, column indices and values, inducing a storage of $7N$ real numbers plus $14N$ integers, i.e., $21 \times 8 \times N$ bytes in 64 bits integer and real arithmetic. In addition, there are six vectors $(x, b, r, p, q, z)$ in Algorithm 1, each of size $N$, requiring $6 \times 8 \times N$ bytes. Hence, PCG uses at least $(21+6) \times 8 \times N$ bytes of memory. Some additional memory is required to store the preconditioner, as detailed below.

## Memory complexity of the preconditioners

**Memory complexity of the multigrid preconditioner.** The geometric multigrid algorithm is based on a series of nested 3D grids where only every second grid point is retained in each direction from one grid level ($\ell$) to the next ($\ell+1$). This means that the number of unknowns from one grid to the next is divided by eight, so that asymptotically the number of unknowns on all grids is $\sum_{\ell=0}^{\infty} N/8^\ell = 8/7N$. This means that the three vectors (current iterate, residual, correction) that need to be stored on each grid level have a total memory requirement of $24/7 \times 8 \times N$ bytes (assuming double precision real coefficients). In addition, the restriction matrices that allow the vectors defined on one grid to be moved to the next have to be stored. If we consider a fully weighted restriction, the corresponding matrix is rectangular with 8 times less rows than columns, with at most 27 non-zero coefficients per row. On the grid hierarchy this is $\sum_{\ell=1}^{\infty} 27 \times N/8^\ell = 27/7 \times N \approx 4 \times N$ coefficients, which translates to $\approx 3 \times 4 \times 8$ bytes if these matrices are stored in COO format (assuming double precision real coefficients and 64 bits integers). The prolongation is the transpose of the constraint and does not need to be stored. Finally, the discretization matrices are computed using a Galerkin approach at each grid level, i.e., $A_\ell = R_\ell A_{\ell-1} R_\ell^T$ where $R_\ell$ is the restriction from level $\ell-1$ to level $\ell$. The matrix $A_\ell$ essentially corresponds to a 27-point grid, so the number of coefficients on the sequence of grid levels is $\sum_{\ell=1}^{\infty} 27 \times N/8^\ell = 27/7N \approx 4 \times N$. Stored in COO format, the memory footprint is $3 \times 4 \times 8 \times N$ bytes. The total memory footprint is $(24/7 + 2 \times 3 \times 4) \times 8 \times N \approx 224 \times N$.

**Memory complexity of the domain decomposition preconditioner.** For the domain decomposition preconditioner, we consider a two-stage Schwarz method based on the GENEIO approach$^{3,47}$. Formally, the preconditioner writes:

$$M = \sum_{i=1}^{P} R_{i}^{T} (R_{i} A R_{i}^{T})^{-1} R_{i} + U_{0} (U_{0}^{T} A U_{0})^{-1} U_{0}^{T},$$

where $P$ denotes the number of subdomains, $R_i$ is the canonical restriction that maps the entries of a vector defined on the full domain to its values in the subdomain number $i$, and the columns of $U_0$ span the coarse space defining the second level of the preconditioner. Note that $R_i$ is required for the mathematical formalism, but is only used for the data mapping in a parallel implementation, so it has no additional memory cost within the iteration. In the GENEIO approach, $U_0$ has $n_{ev} \times P$ columns; each subdomain contributes $n_{ev}$ local vectors resulting from the solution of a local generalised eigenproblem. For the sake of complexity, we neglect the overlap between the subdomains, we consider a 3D decomposition of the cube and one process per subdomain, so that $P$ also corresponds to the number of computing nodes. Each of the $P$ subdomains lies on a subcube of side $m_{dom} = m \times P^{-1/3}$. A sparse direct solver is used to solve both the local problem defined on the subcube and the coarse space problem $A_0 = (U_0^T A U_0)$. In both cases the matrices correspond to the solution of a linear system defined by a 3D-stencil on a Cartesian grid inside a cube. The computational complexity of a sparse direct solver on a cube of side $d$ is$^{21}$: $d^4 \times 8$ bytes in memory. We can now calculate the memory footprint of the two-level preconditioner. For each subdomain we have to store the factors of $(R_i A R_i^T)$, which is $m_{dom}^4 \times 8$ bytes, the local column vectors contributing to $U_0$, which is $m_{dom}^3 \times n_{ev}$. The sum over the subdomain is $(P \times m_{dom}^4 + P \times m_{dom}^3 \times n_{ev}) \times 8 = (P \times (N/P)^{4/3} + N \times n_{ev}) \times 8$ bytes. We also need to store the Cholesky factor of $A_0$, which is a sparse matrix of dimension $(P \times n_{ev}) \times (P \times n_{ev})$, resulting from the discretisation of a template operator on a coarse cube of side $(n_{ev} \times P)^{1/3}$. Its memory footprint is $(n_{ev} \times P)^{4/3} \times 8$ bytes. In total, the domain decomposition preconditioner requires $(P \times (N/P)^{4/3} + N \times n_{ev} + (n_{ev} \times P)^{4/3}) \times 8$ bytes. Usually $n_{ev} \approx 10$ is sufficient to ensure good numerical scalability of the DDM preconditioner. We therefore consider $n_{ev} = 10$ in the evaluation.

## Scenarios

We now proceed to consider four real-life scenarios based on large-scale experiments from the literature$^{31,32,39}$. It should be noted that there is no exact correspondence between the parameters reported in these articles and the input parameters of Section 13. It is therefore necessary to extrapolate some values. While we have endeavoured to provide a rationale for our estimations, alternative extrapolations could yield equally valid results. In the first two scenarios, the use of a multigrid preconditioner is assumed; in contrast, the latter two scenarios assume the use of a DDM preconditioner.

**Scenarios 1 and 2 using a multigrid preconditioner.** Scenario 1 is motivated by$^{32}$, where a Stokes problem of order $N = 3.6 \times 10^{12}$ unknowns is being processed on $P = 3 \times 10^3$ computing nodes. Although the original article employs multigrid as a solver (as opposed to a preconditioner for PCG), the reported timings indicate that a PCG iteration can be estimated as follows: The application of the preconditioner requires 10 seconds, the matrix-vector product 2 seconds, and the remaining operations 1 second. This yields an estimation of $I = 13$ seconds for the duration
of a complete PCG iteration. As the cost of the numerical verification is dominated by the cost of the matrix-vector product, we estimate $V_c = 2$ seconds. For the other values, we deduce them using the size $N$ of the problem, the time of the computation and the balance that appears in Scenario 1. Hence, the cost to copy vectors in the local memory means that $C_{cm} = 0.5$ seconds. The fail-stop checkpoint is global; the filesystem bandwidth of the target machine is typically within the range $[4.6 \cdot 10^{12}, 3.7 \cdot 10^{13}]$ byte per second, depending on the technology used to take and store the checkpoint⁵, so for the overall amount of memory to save, we take $C_{fs} = 180$ seconds. Finally, the checksum cost implies that $V_m = 6$ seconds. In all scenarios, we use $R_{cm} = C_{cm}$ and $R_{fs} = C_{fs}$.

Scenario 2 is motivated by³⁹, where a mantle convection simulation problem of order $N = 1.1 \times 10^{13}$ is processed on $P = 2 \times 10^4$ nodes. Similar extrapolations lead us to estimate a cost of $I = 110$ seconds for a PCG iteration, which can be decomposed into 85 seconds for the preconditioner, 17 seconds for the matrix-vector product, and 8 seconds for the other operations. The numerical verification being roughly the cost of one matrix-vector product, we consider $V_c = 17$ seconds. Then, although the size of the problem is 3.3 times larger than in Scenario 2, there are 6.7 times more nodes. Hence, the local checkpoint is quick, $C_{cm} = 0.25$ seconds. The memory verification takes $V_m = 3$ seconds. Finally, the fail-stop checkpoint is more expensive because of the number of nodes and the unknowns of the problem: $C_{fs} = 540$ seconds.

### Scenarios 3 and 4 using a DDM preconditioner.

Scenarios 3 and 4 are the domain decomposition counterpart of the previous two scenarios that involved a multigrid preconditioner. The memory footprint changes as described earlier. The computation times are estimated using parallel experiments performed at scale and reported in³¹ page 99. These experiments essentially show that for stiff problems, the domain decomposition preconditioner executes four times fewer iterations, but is only twice as fast, which means that an iteration of PCG with domain decomposition takes twice as long as an iteration of PCG with a multigrid preconditioner. A summary of all the parameter values for each scenario is given in Table 2.

## Error probability distributions

There is limited information about potential values for the probability distributions of the three error sources. Let MTBF($fs$), MTBF($mem$) and MTBF($calc$) denote the MTBF for each error type; it is the inverse of the parameters $\lambda_{fs}$, $\lambda_{mem}$ and $f_{calc}$ of their probability distributions. On the

one hand, there are good sources of data for fail-stop errors in the literature. The MTBF for fail-stop errors of recent leadership HPC machines is usually of a few hours⁴,²⁶,³⁴,⁴²,⁴⁵. We focus on a range of 1h to 8h for MTBF($fs$); experiments with a range up to 24h were also done, leading to flat plots and just slightly more iterations between checkpoints. Hence, we focus on the most interesting part of the plots. On the other hand, there is no extensive study on memory and corruption errors in HPC to the best of our knowledge. A few studies have artificially modified the environment of a few nodes of HPC systems to increase the probability of memory and CPU failures due to particles of high energy interacting with the electronics (see¹⁸,³⁸). They then scale down the probability of failure linearly with the chance of normal encounter of such events, and scale it up linearly with the number of nodes. Other studies have disabled ECC mechanisms at scale on significant platforms and studied how many real failures were observed (see⁸). These studies show that the MTBF of unrecoverable and/or undetected errors vary significantly depending on the hardware used, and even on the application running on the machine.

We did not find studies on recent hardware used in exascale supercomputers, and resort to scaling the error rates observed in previous machines to provide a rough estimate of reasonable values for the MTBF($calc$) and MTBF($mem$). For MTBF($mem$), we could assume that silent memory corruptions scale with the amount of memory available. Using this heuristic to scale the observation of¹⁸ on the Sequoia supercomputer to the Frontier supercomputer would lead to an MTBF($mem$) of about 2.25h, and about 5h for the observations of³⁸ on the Roadrunner supercomputer. For MTBF($calc$), we assume that the risk of failure is proportional to the amount of computation, which leads to much lower values (as computation capability has increased much faster than the amount of memory available).

It is to be noted, however, that we lack precise values to estimate MTBF($calc$) and MTBF($mem$), as the technology has evolved. Authors of¹⁸ also note that the computational intensity and type of operation of the application strongly impacts these values (with a factor 50 observed between the most reliable application and the least reliable one). Thus, we use these studies to give us an idea of the order of magnitude of risk for the different MTBFs, and study a large range of possible values for these different parameters.

In order to provide a unified scaling mechanism, we introduce a new parameter $x$, that links the different failure probabilities together. We say that:

$$MTBF(fs) = x, MTBF(mem) = \frac{x}{2}, \text{ and } MTBF(calc) = \frac{x}{20}.$$

$$\text{It means that } \lambda_{fs} = \frac{1}{x}, \lambda_{mem} = \frac{2}{x}, \text{ and } f_{calc} = e^{-\frac{20}{x}}.$$

We made this choice because if the majority of errors are fail-stop errors, then the system would either be really secure or not reliable at all, and our approach would not help in these cases. Moreover, even if we say that there are computation errors and memory errors, it can be difficult in the execution to differentiate between the two of them: in practice, there is a significant chance that what is originally a memory error is caught during the computation verification as it is the first verification to be encountered. In this work, we consider the case when the MTBF of computation errors is shorter than

<table><caption>**Table 2.** Parameter values (in seconds) for the different scenarios with PCG.</caption><thead><tr><th>Scenario</th><th>Scenario 1</th><th>Scenario 2</th><th>Scenario 3</th><th>Scenario 4</th></tr></thead><tbody><tr><td>Pb size (N)</td><td>3.6 × 10<sup>12</sup></td><td>1.1 × 10<sup>13</sup></td><td>3.6 × 10<sup>12</sup></td><td>1.1 × 10<sup>13</sup></td></tr><tr><td>#nodes (P)</td><td>3 × 10<sup>3</sup></td><td>2 × 10<sup>4</sup></td><td>3 × 10<sup>3</sup></td><td>2 × 10<sup>4</sup></td></tr><tr><td>I</td><td>13</td><td>110</td><td>26</td><td>220</td></tr><tr><td>V<sub>c</sub></td><td>2</td><td>17</td><td>2</td><td>17</td></tr><tr><td>V<sub>m</sub></td><td>6</td><td>3</td><td>22.8</td><td>11.4</td></tr><tr><td>C<sub>cm</sub> = R<sub>cm</sub></td><td>0.5</td><td>0.25</td><td>1.9</td><td>0.95</td></tr><tr><td>C<sub>fs</sub> = R<sub>fs</sub></td><td>180</td><td>540</td><td>684</td><td>2052</td></tr></tbody></table>
the MBTF of memory errors. This was the main rationale to
design the hierarchical pattern shown in Figure 1.

Results

The code for the evaluation is publicly available at https://doi.org/10.5281/zenodo.14572919. This code computes the optimal values for $n_{vc}$, $n_{cm}$, $n_{fs}$ that minimize the slowdown *S* using the exhaustive search described in Section, given the parameters of the different scenarios of Table 2.

Figures 2 to 5 show the slowdown ratio *S*, with the optimal values for *n*<sub>vc</sub>, *n*<sub>cm</sub>, *n*<sub>fs</sub>, and with the naive values (*n*<sub>vc</sub> = 1, *n*<sub>cm</sub> = 1, *n*<sub>fs</sub> = 1) in some cases, as a function of the error parameter *x*, which represents the MTBF of fail-stop errors (MTBF(*fs*)), the other MTBFs (MTBF(*calc*) and MTBF(*mem*)) being scaled accordingly, as described in Section 13.

When *x* increases, all error types are less likely to strike. The y-axis is the slowdown, which is a ratio of the expected completion time of the strategy to the execution time without faults and without fault-tolerance mechanism. We measure the average of added costs consisting on the impact of errors on the algorithm with the recoveries and the added cost of the checkpoints and verifications that protect the algorithm. A value of one means that there is no added cost. It is safe to say that the slowdown decreases when the MTBF increases as less errors means less recoveries and there is less need to take verifications and checkpoints as often.

The triplets of numbers that appears on the plots are ($n_{vc}$, $n_{cm}$, $n_{fs}$) and correspond respectively to the number of iterations in a chunk, the number of chunks and the number of segments in the pattern. The product of these three variables is the number of iterations for the optimal pattern as a function of the error probabilities.

The result of Scenario 1 (Figure 2) shows a slowdown less than 2 for all values above 1h of fail-stop MTBF, and quickly drops to below 1.5 for values above 3h (so the application is at most 50% slower with faults and the optimal fault tolerance pattern than in a fully reliable machine). However, in the context of the PCG algorithm, that scenario exposes a possible limit of the approach via the number of iterations in a pattern. For *x* = 4 hours, a full pattern requires 3 × 2 × 22 = 132 iterations. A simple CG algorithm in exact arithmetic takes at most the size of the problem as the number of iterations to converge. It converges even faster with the help of a preconditioner. In the case we consider, this means that in practice, only a few full patterns are expected to execute.

The naive approach consists of checkpointing at every
iteration. We do not show the performance of the naive
approach for this scenario, as it exhibits a slowdown
systematically above 16 for the entire set of x values. This
reflects the efficiency of the optimized approach.

Scenario 2 (Figure 3) shows a similar behavior, with a
higher slowdown. As the number of nodes and the problem
size increase, so does the duration of all other parameters,
which directly impacts the overheads, resulting in a higher
slowdown. The problem size, however, may require to
distribute the application over more nodes. On a machine
with the same levels of reliability, this leads to a minimum

slowdown about twice larger (but still under a factor two for
fairly reliable machines).

In this case, the optimal pattern falls back to taking
verifications and in-memory checkpoints at each iteration
($n_{vc}$ = 1 and $n_{cm}$ = 1 for all optimal patterns), and only
separates full fail-stop filesystem based checkpoints by an
increasing number of iterations as failures become more
rare. This is easily explained by the relative cost of the
different fault tolerance techniques: compared to the high
filesystem checkpoint cost of 540s, the other overheads are
small and worth paying for the added benefits of avoiding
full rollbacks.

This figure is also the first to show the slowdown of the
naive approach. Compared to the optimal pattern, this would
lead to a slowdown about three times higher than what is
achievable, and above 6 times slower than the failure-free
execution.

Scenarios 3 and 4 (Figures 4 and 5) consider DDM
preconditioner, which exhibit higher overheads for the
fail-stop checkpoints and memory protection. The cost of
each iteration is also increased compared to the previous
approach. As a consequence, the slowdowns are also
increased, up to a factor 3.5 for Scenario 3, and above
10 for Scenario 4, if the system is very unreliable (1h of
fail-stop MTBF). However, in more reliable setups (8h of
fail-stop MTBF, which corresponds to the type of MTBF
expected today on large supercomputers), these slowdown
factors become more reasonable (about 1.75 for Scenario 3,
and 2.5 for Scenario 4). As for Figure 2, the naive approach
is too slow to appear on the figure; in this case, for all
considered MTBF(*fs*) values, the naive approach slowdown
is above 29.

At a given problem size, a lower number of nodes incurs
a lower slowdown. This is mostly driven by the increase
in iteration cost and in the increase in checkpointing time
for the filesystem checkpoint. If a scalable approach for
checkpointing can be implemented, most of these overheads
could be absorbed providing better efficiency for this
hierarchical approach.

Finally, note that while the slowdown of the naive
approach turns out to be very high even with large values of
MTBF (i.e., few failures), the slowdowns using the optimal
pattern are reasonable, in particular with large values of
MTBF. For an MTBF of 12h, we observe a slowdown of no
more than 1.2 in scenarios 1 and 2, and 1.3 in scenarios 3
and 4.

All previous experiments set a fixed relation between
MTBF(*fs*), MTBF(*calc*), and MTBF(*mem*). This was done
to allow us representing curves in a simple space. In Figure 6,
we consider a single point of Scenario 2 (Figure 3), when
MTBF(*fs*) = *x* = 4h, and we vary MTBF(*mem*) relative
to MTBF(*fs*) on the x-axis and MTBF(*calc*) relative to
MTBF(*fs*) on the y-axis. The color shade represents the
minimum Slowdown *S* for this set of parameters. We only
consider the scenarios presented in this paper, i.e., when
MTBF(*calc*) < MTBF(*mem*) < MTBF(*fs*).

As expected, if either MTBF(*calc*) or MTBF(*mem*)
increases, the optimal pattern requires to take more frequent
in-memory checkpoints and verifications, and the algorithm
implies taking more frequent rollbacks to these checkpoints,
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Results for Scenario 1.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Results for Scenario 2.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Results for Scenario 3.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Results for Scenario 4.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Scenario 2, $x = 4$ hours. Impact of the coefficient of silent error probabilities.</figcaption></figure>

leading to higher overheads. More interestingly, both parameters seem to impact the slowdown in a similar manner, showing that the method can adapt to very variable systems, where either components are more or less reliable.

We have allowed in this figure MTBF(*mem*) and MTBF(*calc*) to become highly unreliable, up to 50 times lower than MTBF(*fs*), meaning there would be a computational or memory failure every 5 minutes. This is to show the robustness of the approach, even though we do not expect such systems to be used in production. At values that we hope to be more realistic, the slowdown is below 1.25, making the method efficient for the benefits of obtaining a reliable result on an unreliable machine.

## Conclusion

The main contribution of this work is the design and analysis of a new flexible approach to protect numerical iterative algorithms against three types of errors: fail-stop, memory and computation errors. We introduce a hierarchical pattern and compute the mean time per iteration within that pattern, enabling to optimize it. Finally, we instantiate the model with different scenarios based on the Preconditioned Conjugate Gradient (PCG) algorithm, demonstrating that the use of the proposed pattern allows us to significantly improve the expected execution time, compared to a naive approach where we would checkpoint at each iteration.

To the best of our knowledge, while simpler patterns have been introduced for simpler models in the literature, this is the first time that all three error types are dealt with, at the price of a very complicated derivation for the expected time per iteration. This work is the first study combining all three error types and opens several future research directions.

In particular, we plan to extend the study to the case with imperfect verifications, which greatly complicates the analysis. Also, we would like to instantiate the pattern to other iterative algorithms, so as to further demonstrate the usefulness of this generic pattern-based approach.

## Acknowledgements

This work was supported in part by the PEPR NumPEX (https://numpex.org). We thank the anonymous reviewers for their comments and suggestions, which helped improve the final version of the paper.

## References

1. Agullo E, Altenbernd M, Anzt H, Bautista-Gomez L, Benacchio T, Bonaventura L, Bungartz H, Chatterjee S, Ciorba FM, DeBardeleben N, Drzisga D, Eibl S, Engelmann C, Gansterer WN, Giraud L, Göddeke D, Heisig M, Jézéquel F, Kohl N, Li XS, Lion R, Mehl M, Mycek P, Obersteiner M, Quintana-Ortí ES, Rizzi F, Rüde U, Schulz M, Fung F, Speck R, Stals L, Teranishi K, Thibault S, Thönnes D, Wagner A and Wohlmuth BI (2022) Resiliency in numerical algorithm design for extreme scale simulations. *Int. J. High Perform. Comput. Appl.* 36(2): 251–285. DOI:10.1177/10943420211055188. URL https://doi.org/10.1177/10943420211055188.

2. Agullo E, Cools S, Yetkin EF, Giraud L, Schenkels N and Vanroose W (2020) On soft errors in the conjugate gradient method: Sensitivity and robust numerical detection. *SIAM Journal on Scientific Computing* 42(6): C335–C358. DOI: 10.1137/18M122858X.

3. Agullo E, Giraud L and Poirel L (2019) Robust preconditioners via generalized eigenproblems for hybrid sparse linear solvers.
SIAM Journal on Matrix Analysis and Applications 40(2): 417-439. DOI:10.1137/17m1153765. URL http://dx.doi.org/10.1137/17M1153765.

4. Ashraf RA and Engelmann C (2018) Analyzing the impact of system reliability events on applications in the titan supercomputer. In: *2018 IEEE/ACM 8th Workshop on Fault Tolerance for HPC at eXtreme Scale (FTXS)*. pp. 39–48. DOI: 10.1109/FTXS.2018.00008.

5. Atchley S, Zimmer C, Lange J, Bernholdt D, Melesse Vergara V, Beck T, Brim M, Budiardja R, Chandrasekaran S, Eisenbach M, Evans T, Ezell M, Frontiere N, Georgiadou A, Glenski J, Grete P, Hamilton S, Holmen J, Huebl A, Jacobson D, Joubert W, Mcmahon K, Merzari E, Moore S, Myers A, Nichols S, Oral S, Papatheodore T, Perez D, Rogers DM, Schneider E, Vay JL and Yeung PK (2023) Frontier: Exploring exascale. In: *Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis*, SC ’23. Association for Computing Machinery. ISBN 9798400701092.

6. Bautista-Gomez L, Benoit A, Cavelan A, Raina S, Robert Y and Sun H (2016) Coping with recall and precision of soft error detectors. *J. Parallel and Distributed Computing* 98: 8–24.

7. Bautista-Gomez L, Zyulkyarov F, Unsal O and McIntosh-Smith S (2016) Unprotected computing: A large-scale study of dram raw error rate on a supercomputer. In: *SC '16: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis*. pp. 645–655. DOI:10.1109/SC.2016.54.

8. Bautista-Gomez L, Zyulkyarov F, Unsal O and McIntosh-Smith S (2016) Unprotected computing: A large-scale study of dram raw error rate on a supercomputer. In: *SC '16: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis*. pp. 645–655. DOI:10.1109/SC.2016.54.

9. Benoit A, Cavelan A, Le Fèvre V, Robert Y and Sun H (2017) Towards Optimal Multi-Level Checkpointing. *IEEE Transactions on Computers* 66(7): 1212–1226. URL https://inria.hal.science/hal-02082416.

10. Benoit A, Cavelan A, Robert Y and Sun H (2016) Optimal resilience patterns to cope with fail-stop and silent errors. In: *IPDPS’2016, the 30th IEEE International Parallel and Distributed Processing Symposium*. IEEE Computer Society Press.

11. Benson AR, Schmit S and Schreiber R (2014) Silent error detection in numerical time-stepping schemes. *Int. J. High Performance Computing Applications*.

12. Boeuf JP (2017) Tutorial: Physics and modeling of hall thrusters. *Journal of Applied Physics* 121(1). DOI:10.1063/1.4972269. URL http://dx.doi.org/10.1063/1.4972269.

13. Bosilca G, Delmas R, Dongarra J and Langou J (2009) Algorithm-based fault tolerance applied to high performance computing. *J. Parallel Distrib. Comput.* 69(4): 410–416.

14. Bronevetsky G and de Supinski B (2008) Soft error vulnerability of iterative linear algebra methods. In: *ICS*. ACM.

15. Cappello F, Geist A, Gropp W, Kale S, Kramer B and Snir M (2014) Toward exascale resilience: 2014 update. *Supercomputing frontiers and innovations* 1(1).

16. Chen Z (2013) Online-ABFT: An online algorithm based fault tolerance scheme for soft error detection in iterative methods. In: *Proc. PPoPP*. pp. 167–176.

17. Chen Z (2013) Online-abft: An online algorithm based fault tolerance scheme for soft error detection in iterative methods. *ACM SIGPLAN Notices* 48: 167–176. DOI:10.1145/2442516. 2442533.

18. Cher CY, Gupta MS, Bose P and Muller KP (2014) Understanding soft error resiliency of blue gene/q compute chip through hardware proton irradiation and software fault injection. In: *SC '14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis*. pp. 587–596. DOI:10.1109/SC.2014.53.

19. Chorin AJ (1968) Numerical solution of the navier-stokes equations. *Mathematics of Computation* 22(104): 745–762. URL http://www.jstor.org/stable/2004575.

20. Daly JT (2006) A higher order estimate of the optimum checkpoint interval for restart dumps. *FGCS* 22(3): 303–312.

21. Davis TA (2006) *Direct Methods for Sparse Linear Systems*. Society for Industrial and Applied Mathematics. ISBN 9780898718881. DOI:10.1137/1.9780898718881. URL http://dx.doi.org/10.1137/1.9780898718881.

22. Di S and Cappello F (2016) Adaptive Impact-Driven Detection of Silent Data Corruption for HPC Applications. *IEEE Transactions on Parallel and Distributed Systems* 27(10): 2809–2823. DOI:10.1109/TPDS.2016.2517639.

23. Fasi M, Langou J, Robert Y and Ucar B (2016) A backward/forward recovery approach for the preconditioned conjugate gradient method. *J. Computational Science* 17: 522–534.

24. Golub GH and Loan CFV (2013) *Matrix computations*. 4 edition. Johns Hopkins.

25. Halko N, Martinsson PG and Tropp JA (2011) Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. *SIAM Review* 53(2): 217–288.

26. Herault T and Robert Y (eds.) (2015) *Fault-Tolerance Techniques for High-Performance Computing*, Computer Communications and Networks. Springer Verlag.

27. Heroux M and Hoemmen M (2011) Fault-tolerant iterative methods via selective reliability. Research report SAND2011-3915 C, Sandia Nat. Lab.

28. Hestenes MR and Stiefel E (1952) Methods of conjugate gradients for solving linear systems. *Journal of Research of the National Bureau of Standards* 46(6): 409–436.

29. Huang KH and Abraham JA (1984) Algorithm-based fault tolerance for matrix operations. *IEEE Trans. Comput.* 33(6): 518–528.

30. Hwang AA, Stefanovici IA and Schroeder B (2012) Cosmic rays don’t strike twice: understanding the nature of DRAM errors and the implications for system design. *SIGARCH Comput. Archit. News* 40(1): 111–122.

31. Jolivet P (2006) *Méthodes de décomposition de domaine. Application au calcul haute performance*. PhD Thesis, Université de Grenoble.

32. Kohl N and Rüde U (2022) Textbook efficiency: Massively parallel matrix-free multigrid for the stokes system. *SIAM J. Scientific Computing* 44(2): C124–C155. DOI:10.1137/20M1376005.

33. Li S, Di S, Zhao K, Liang X, Chen Z and Cappello F (2021) Resilient error-bounded lossy compressor for data transfer. In: *Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis*, SC ’21. New York, NY, USA: Association for
Computing Machinery. ISBN 9781450384421. DOI:10.1145/
3458817.3476195. URL https://doi.org/10.1145/
3458817.3476195.

34. Liu RT and Chen ZN (2018) A large-scale study of failures on petascale supercomputers. *Journal of Computer Science and Technology* 33(1): 24–41. DOI:10.1007/s11390-018-1806-7. URL https://jcst.ict.ac.cn/en/article/doi/10.1007/s11390-018-1806-7.

35. Lyons RE and Vanderkulk W (1962) The use of triple-modular redundancy to improve computer reliability. *IBM J. Res. Dev.* 6(2): 200–209.

36. Marshall J, Hill C, Perelman L and Adcroft A (1997) Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling. *Journal of Geophysical Research: Oceans* 102(C3): 5733–5752. DOI:10.1029/96jc02776. URL http://dx.doi.org/10.1029/96JC02776.

37. Meurant G (2022) Detection and correction of silent errors in the conjugate gradient algorithm. *Numerical Algorithms* 92. DOI:10.1007/s11075-022-01380-1.

38. Michalak SE, DuBois AJ, Storlie CB, Quinn HM, Rust WN, DuBois DH, Modl DG, Manuzzato A and Blanchard SP (2012) Assessment of the impact of cosmic-ray-induced neutrons on hardware in the roadrunner supercomputer. *IEEE Transactions on Device and Materials Reliability* 12(2): 445–454. DOI: 10.1109/TDMR.2012.2192736.

39. Mohr M, Rüde U, Wohlmuth B and Bunge HP (2023) Challenges for mantle convection simulations at the exa-scale: Numerics, algorithmics and software. In: *Impact of Scientific Computing on Science and Society*. Cham: Springer, pp. 75–92.

40. Moody A, Bronevetsky G, Mohror K and De Supinski BR (2010) Design, modeling, and evaluation of a scalable multi-level checkpointing system. In: *SC’10: Proceedings of the 2010 ACM/IEEE International Conference for High Performance Computing, Networking, Storage and Analysis*. IEEE, pp. 1–11.

41. O’Gorman T (1994) The effect of cosmic rays on the soft error rate of a DRAM at ground level. *IEEE Trans. Electron Devices* 41(4): 553–557.

42. Rojas E, Meneses E, Jones T and Maxwell D (2019) Analyzing a five-year failure record of a leadership-class supercomputer. In: *2019 31st International Symposium on Computer Architecture and High Performance Computing (SBAC-PAD)*. pp. 196–203. DOI:10.1109/SBAC-PAD.2019.

00040.

43. Sao P and Vuduc R (2013) Self-stabilizing iterative solvers. In: *Proceedings of the Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems*. New York, NY, USA: Association for Computing Machinery. ISBN 9781450325080. DOI:10.1145/2530268.2530272. URL https://doi.org/10.1145/2530268.2530272.

44. Shantharam M, Srinivasmurthy S and Raghavan P (2012) Fault tolerant preconditioned conjugate gradient for sparse linear system solution. In: *ICS*. ACM.

45. Shoji F, Matsui S, Okamoto M, Sueyasu F, Tsukamoto T, Uno A and Yamamoto K (2015) Long term failure analysis of 10 peta-scale supercomputer. In: *HPC in Asia Poster, ISC*.

46. Snir M, Wisniewski RW, Abraham JA, Adve SV, Bagchi S, Balaji P, Belak J, Bose P, Cappello F, Carlson B, Chien AA, Coteus P, Debardeleben NA, Diniz PC, Engelmann C, Erez M, Fazzari S, Geist A, Gupta R, Johnson F, Krishnamoorthy S, Leyffer S, Liberty D, Mitra S, Munson T, Schreiber R, Stearley J and Hensbergen EV (2014) Addressing failures in exascale computing. *Int. J. High Perform. Comput. Appl.* 28(2): 129–173.

47. Spillane N, Dolean V, Hauret P, Nataf F, Pechstein C and Scheichl R (2013) Abstract robust coarse spaces for systems of pdes via generalized eigenproblems in the overlaps. *Numerische Mathematik* 126(4): 741–770. DOI:10.1007/s00211-013-0576-y. URL http://dx.doi.org/10.1007/s00211-013-0576-y.

48. Sridharan V, DeBardeleben N, Blanchard S, Ferreira KB, Stearley J, Shalf J and Gurumurthi S (2015) Memory errors in modern systems: The good, the bad, and the ugly. In: *20th Int. Conf. on Architectural Support for Programming Languages and Operating Systems (ASPLOS)*. ACM, pp. 297–310.

49. The TOP500 team (2024) HPCG. https://top500.org/lists/hpcg/2024/06/.

50. Young JW (1974) A first order approximation to the optimum checkpoint interval. *Comm. of the ACM* 17(9): 530–531.

51. Ziegler J, Nelson M, Shell J, Peterson R, Gelderloos C, Muhlfeld H and Montrose C (1998) Cosmic ray soft error rates of 16-Mb DRAM memory chips. *IEEE Journal of Solid-State Circuits* 33(2): 246–252.

52. Ziegler JF, Curtis HW, Muhlfeld HP, Montrose CJ and Chin B (1996) IBM experiments in soft fails in computer electronics. *IBM J. Res. Dev.* 40(1): 3–18.