{'fraction_non_alphanumeric': 0.025029618381751738, 'fraction_numerical': 0.009044035441939628, 'mean_word_length': 5.299416618489515, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Resiliency in Numerical Algorithm Design for Extreme Scale Simulations

This work is based on the seminar titled ``Resiliency in Numerical Algorithm Design for Extreme Scale Simulations'' held March 1-6, 2020 at Schloss Dagstuhl, that was attended by all the authors. Naive versions of conventional resilience techniques will not scale to the exascale regime: with a main memory footprint of tens of Petabytes, synchronously writing checkpoint data all the way to background storage at frequent intervals will create intolerable overheads in runtime and energy consumption. Forecasts show that the mean time between failures could be lower than the time to recover from such a checkpoint, so that large calculations at scale might not make any progress if robust alternatives are not investigated. More advanced resilience techniques must be devised. The key may lie in exploiting both advanced system features as well as specific application knowledge. Research will face two essential questions: (1) what are the reliability requirements for a particular computation and (2) how do we best design the algorithms and software to meet these requirements? One avenue would be to refine and improve on system- or application-level checkpointing and rollback strategies in the case an error is detected. Developers might use fault notification interfaces and flexible runtime systems to respond to node failures in an application-dependent fashion. Novel numerical algorithms or more stochastic computational approaches may be required to meet accuracy requirements in the face of undetectable soft errors. The goal of this Dagstuhl Seminar was to bring together a diverse group of scientists with expertise in exascale computing to discuss novel ways to make applications resilient against detected and undetected faults. In particular, participants explored the role that algorithms and applications play in the holistic approach needed to tackle this challenge.


Introduction
Numerical simulation is the third pillar in science discovery at the same level as theory and experiments. To cope with the ever demanding computational resources needed by complex simulations, the computational power of high performance computing systems continues to increase by using an ever larger number of cores or by specialized processing. On the technological side, the continuous shrinking of transistor geometry and the increasing complexity of these devices affect their sensitivity to external effects and thus diminish their reliability. A direct consequence is that High Performance Computing (HPC) applications are increasingly prone to errors. Therefore the design of resilient systems and numerical algorithms that are able to exploit possible unstable HPC platforms has became a major concern in the computational science community. To tackle this critical challenge on the path to extreme scale computation an holistic and multidisciplinary approach is required that needs to involve researchers from various scientific communities ranging from the hardware/system community to applied mathematics for the design of novel numerical algorithms. In this article, we summarize and report on the outcomes of a Dagstuhl seminar held March 1-6, 2020, 1 on the topic Resiliency in Numerical Algorithm Design for Extreme Scale Simulations. We should point out that, although error and resiliency was already quoted by J. von Neumann in his first draft report on EDVAC [260, P.1, Item 1.4], it became again a central concern for the HPC community in the late 2000' when the availability of the first exascale computers was envisioned for the forthcoming decades. In particular, several workshops were organized in the IESP (International Exascale Software Project) and EESI (European Exascale Software Initiative) framework [52]. The hardware/system resilience community has previously defined terminology related to how faults, errors, and failures occur on computing systems [18]. In this article our focus is less on the cause of an error (or the underlying fault), and more on how an error presents itself at the algorithmic level (or layer), impacting algorithms and applications. We thus simplify the terminology often used in the hardware resilience and fault-tolerance community by not using terms like soft error or hard error, and generally do not concern ourselves with the reproducibility of an error (e.g., transient, intermittent or permanent). This abstraction keeps the algorithmic techniques discussed herein general and applicable to a variety of fault models, current architectures, and hopefully of use in future technologies.
To this end, we broadly categorize errors presenting themselves to the algorithmic layer as either detected or undetected. Note that this categorization does not mean an error is undetectable but rather that when it reached the algorithmic layer it was not detected by earlier layers (e.g., hardware, operating system or middleware/system software). This suggests the algorithmic layer has the opportunity to detect a previously undetected error and, if possible, to deploy mitigation methods to make the algorithm resilient; effectively transforming an undetected error at the algorithmic layer into a detected error. This in turn may result in a failure if the algorithm is unable to handle it. For example, an undetected data corruption which results in an application accessing an incorrect memory address may be detectable by the algorithm but it may not be possible for the algorithm to recover what the original memory address was and it may be forced to fail. If the algorithm could not detect the corruption before accessing the memory region, this would conventionally end in a failure (e.g., SIGSEGV issued by the operating system).
Many computing-intensive scientific applications that are dependent on HPC performance upgrades can end up with disrupted schedules because of lack of resilience. A typical example is related to current efforts towards exascale numerical weather prediction [31,32]. On one side, regular upgrades in weather forecast models in operations at weather centres and their spatial resolution have gone hand in hand with expanding computational resources. On the other side, scientific and socioeconomic significance of forecasts crucially hinges on tight time-bound computing schedules and timely forecast dissemination, most notably for high-impact weather events. Current disk-checkpointing schedules still take up acceptable portions of forecast runtimes, but are hardly sustainable -indeed, they already saturate file systems bandwidth. In addition, many weather forecast codes feature preconditioned iterative solvers of linear systems with several hundred thousand unknowns, many thousand times per run. Such components represent vulnerable points in a context of increasingly frequent detected and undetected errors. Novel low-overhead solutions to enhance algorithmic fault-tolerance or provide higher-level system resilience are therefore in high demand in this and other fields where nonlinear dynamics is simulated.
In this article we take a different approach at the classification of errors in HPC systems. In general, we try to divide errors in two main groups, those that are detected and corrected by the hardware/system (which is the focus of Section 2) and those that are detected and sometimes corrected by the numerical algorithms (Section 3). However, the HPC resilience ecosystem is not black and white, but it rather shows a wide palette of greys in between, with multiple fault tolerance tools implemented at the middleware level that are assisted by the applications/algorithms and vice-versa. Figure 1 shows this wide range of different error classifications depending on how much effort is needed at the application/algorithmic level in order to detect/correct the error.  The first category we observe in the leftmost leaf of the tree (blue color) is the case of errors that are both detected and transparently corrected by the hardware/middleware but without any intervention of the applications/algorithms. The clearest example would be a detectable and correctable error in the memory generated by a single bit flip. These types of errors are transparently corrected by the system without any knowledge at the application/algorithmic level that such error mitigation occurred. Other examples could be process replication, system-level checkpointing, process migration, among many others (see Section 2.1).
The second category is the case of errors that are detected at the hardware/system level and are mitigated at the system/middleware level (not at the algorithmic level) but with assistance from the application/algorithm (green color). The most clear example of this is application-based checkpointing libraries, which handle all or most of the data transfers between the compute nodes and the Parallel File System (PFS) independently from the application, but it gets hints from it to know what datasets need to be checkpointed and when should the checkpoint happen.
Other relevant examples are fault tolerant message passing programming models and resilient asynchronous tasks. We divide these sections in those approaches that require just a minor addition in the application code versus those that require a complete change in the programming paradigm (see Section 2.2).
The other leaves of the tree (red color) correspond to those errors that cannot be corrected or mitigated at the hardware/system level and have to be mitigated by changing the algorithm or numerical methods to be able to tolerate those errors. We observe three different types of algorithms in this branch of the tree.
The first type of algorithms focuses on the mitigation of errors that have been detected (first red leaf from left to right), we call them error-aware algorithms. Please note that these algorithms are not in charge of detecting the errors but only of mitigating them. Also, it is important to notice that these algorithms do not depend on how the error was actually detected; it could be hardware/middleware detection as well as algorithmic detection, in the end the process of detection is irrelevant for the mitigation algorithm (see Section 3.2).
The second type of algorithms are those dedicated to the detection of errors that were not detected at the lower levels (fourth leaf from left to right). A good example would be Silent Data Corruption (SDC) errors that pass invisibly through the hardware but then can be caught at the algorithmic level using some numerical techniques (e.g., checksum). These algorithms do not try to mitigate the error per se but only detect it. Once the error has been detected, it can be passed to a error aware algorithm in order to attempt a correction/mitigation (see Section 3.1).
Finally, there also exist algorithms that can operate, tolerate and absorb errors without ever being aware that the error ever occurred (last leaf to the right); we called these, error oblivious algorithms. These are somehow similar to the very first (blue) category, in that the errors are transparently corrected/absorbed, see Section 3.3.
In the following sections we discuss algorithmic and application approaches to address these two categories of errors and distinguish how the approaches vary or are similar. Broadly speaking, the report is divided into two parts. In Section 2 and Section 3 we discuss the state-of-the-art in the areas of infrastructure and algorithms, while in Section 4 we propose possible areas of interest in future research.

System infrastructure techniques for resilience
In this section we describe the state-of-the-art of hardware and system level error detection and mitigation. As previously mentioned, we divide these methods in two categories, the ones that mitigate the error in a completely transparent fashion, and those that require assistance from the algorithmic/application level. The following subsection, Section 2.1, concentrates on the methods falling in the first category. The second category is explored in Section 2.2.

Detected and transparently corrected errors
A wide range of errors can be detected and immediately corrected by various layers in the system, i.e., these errors become masked or absorbed and higher level layers do not have to be involved. The detection/correction mechanisms have an extra cost in terms of storage, processing and energy consumption.
Hardware reliability At the hardware level several techniques exist to detect and correct errors. Most common examples are Error Correcting Codes (ECC) to detect and correct single bit-errors, Cyclic Redundancy Checks (CRC) error correction for network packets or RAID-1 (or higher) for I/O systems. A more comprehensive discussion of these features can be found in the report "Towards Resilient EU HPC Systems: A Blueprint" by Radojkovic et al. [207].
Operating system reliability Operating Systems (OS) have certain capabilities to interact with architectural resilience features, such as ECC and machine check exceptions. OSs are mostly concerned with resource management and error notification. However, some advanced OS resilience solutions exist such as Mini-ckpts [92]. It is a framework that enables application survival despite the occurrence of a fatal operating system failure or crash. It ensures that the critical data describing a process is preserved in persistent memory prior to the failure. Following the failure, the OS is rejuvenated via a warm reboot and the application continues execution effectively making the failure and restart transparent. The mini-ckpts rejuvenation and recovery process is measured to take 3 s to 6 s and has a failure-free overhead of 3 % to 5 % for a number of key HPC workloads.
System-level checkpoint/restart Berkeley Lab Checkpoint/Restart (BLCR) [125] is a system-level checkpoint/restart solution that transparently saves and restores process state. In conjunction with a Message Passing Interface (MPI) [95] implementation, it can transparently save and restore the process states of an entire MPI application. An extension of BLCR [257,262,263] includes enhancements in support of scalable group communication for MPI membership management, reuse of network connections, transparent coordinated checkpoint scheduling, a job pause feature, and full/incremental checkpointing. The transparent mechanism for job pause allows live nodes to remain active and roll back to the last checkpoint, while failed nodes are dynamically replaced by spares before resuming from the last checkpoint. A minimal overhead of 5.6% is reported in case migration takes place, while the regular checkpoint overhead remains unchanged.
The hybrid checkpointing technique [111] alternates between full and incremental checkpoints: At incremental checkpoints, only data changed since the last checkpoint is captured. This results in significantly reduced checkpoint sizes and overheads with only moderate increases in restart overhead. After accounting for cost and savings, the benefits due to incremental checkpoints are an order of magnitude larger than the overheads on restarts.
Silent Data Corruption (SDC) protection FlipSphere [93] is a tunable, transparent Silent Data Corruption (SDC) detection and correction library for HPC applications. It offers comprehensive SDC protection for application program memory using on-demand memory page integrity verification. Experimental benchmarks show that it can protect 50 % to 80 % of program memory with time overheads of 7 % to 55 %.
Proactive fault tolerance using process or virtual machine migration Proactive fault tolerance [88,187,264] prevents compute node failures from impacting running applications by migrating parts of an application, i.e., tasks, processes, or virtual machines, away from nodes that are about to fail. Pre-fault indicators, such as a significant increase in temperature, can be used to avoid an imminent failure through anticipation and reconfiguration. As computation is migrated away, application failures are avoided, which is significantly more efficient than checkpoint/restart if the prediction is accurate enough. The proactive fault tolerance framework consists of process and virtual machine migration, scalable system monitoring and online/offline system health analysis. The process-level live migration supports continued execution of applications during much of process migration and is integrated into an MPI execution environment. Experiments indicate that 1 s to 6.4 s of prior warning are required to successfully trigger live process migration, while similar operating system virtualization mechanisms require 13 s to 24 s. This error oblivious approach complements checkpoint/restart by nearly cutting the number of checkpoints by half when 70% of the faults are handled proactively.
Resiliency using task-based runtime systems Task-based runtime systems have appealing intrinsic features for resiliency due to the fault isolation they provide by design as they have a view of the task flow and dynamically schedule task on computing units (often to minimize the time to solution or energy consumption). Once an error is detected and identified by the hardware or the algorithm, the runtime system can limit its propagation through the application by reasoning about the data dependencies among tasks [176]. For example, one can envision the scenario where an uncorrectable hardware error is detected triggering the runtime system to dynamically redistribute the tasks to the remaining resources available.
Task-based runtime systems can also limit the size of the state needed to be saved to enable restarting computations, when an error is encountered [171,172,254]. In classical checkpoint-restart mechanisms, the size of the checkpoint can become very large for large-scale applications, and managing it can take up a significant portion of the overall execution. A task-based runtime system simplifies the identification of points during the application execution when the state size is small, since only task boundaries need to be considered for saving the state. Further, identification of idempotent tasks can greatly help task-based runtimes to further reduce the overheads by completely avoiding data backups specific to those tasks. Recent works on on-node task parallel programming models suggest that a simple extension of the existing task-based programming framework enables efficient localized recovery [200,236,238].
The checkpointing itself can also be achieved completely asynchronously [171,172,254]. The runtime allows tasks to read data being saved, and only blocks those tasks that attempt to overwrite data being saved. Since the runtime system knows which data will soon be overwritten by some tasks, it can prioritize the writing of the different pieces so as to have as little impact on the execution as possible. At the restarting point, the runtime also has all information to be able to achieve a completely local recovery. The replacement node can restart from the last valid checkpoint of the previously-failed node, while the surviving nodes can just replay the required data exchanges.
With the recent emergence of heterogeneous computing systems utilizing Graphics Processing Units (GPU), the task programming model is being used to offload computation from the Central Processing Unit (CPU) to the GPU. VOCL-FT [202] offers checkpoint/restart for computation offloaded to GPU using OpenCL [118]. It transparently intercepts the communication between the originating process and the local or remote GPU to automatically recover from ECC errors experienced on the GPU during computation. Another preliminary prototype design extends this concept in the context of OpenMP [42] using a novel concept for Quality of Service (QOS) and a corresponding Application Programming Interface (API) [89]. While the programmer is specifying the resilience requirements for certain offloaded tasks, the underlying programming model runtime decides on how to meet them using a QOS contract, such as by employing task-based checkpoint-restart or redundancy.
Resilience via complete redundancy The use of redundant MPI processes for error detection has been widely analyzed in the last decade [55,70,208,273]. Modular redundancy incurs high overhead, but offers excellent error detection accuracy and coverage with few to no false positive or false negatives.
Complete modular redundancy is typically too expensive for actual HPC workloads. However, it can make sense for certain subsystems such as parts of a PFS. The Meta Data Service (MDS) of a networked PFS is a critical single point of failure. An interruption of service typically results in the failure of currently running applications utilizing its file system. A loss of state requires repairing the entire file system, which could take days on large-scale systems, and may cause permanent loss of data. PFSs such as Lustre [67] often offer some type of active/standby fail-over mechanism for the MDS. A solution [128] for the MDS of the Parallel Virtual File System offers symmetric active/active replication using virtual synchrony with an internal replication implementation. In addition to providing high availability, this solution is taking advantage of the internal replication implementation by load balancing MDS read requests, improving performance over the non-replicated MDS.
Resilience via partial redundancy Partial redundancy has been studied to decrease the overhead of complete redundancy [85,234,239,240]. Adaptive partial redundancy has also been proposed wherein a subset of processes is dynamically selected for replication [108]. Partial replication (using additional hardware) of selected MPI processes has been combined with prediction-based detection to achieve SDC protection levels comparable with those of full duplication [37,38,188]. A Selective Particle Replication approach for meshfree particle-based codes protects the data of the entire application (as opposed to a subset) by selectively duplicating 1 % to 10 % of the computations within processes incurring a 1 % to 10 % overhead [57].
Resilience via complete and/or partial redundancy RedMPI [91] enables a transparent redundant execution of MPI applications. It sits between the MPI library and the MPI application, utilizing the MPI Profiling Interface (PMPI) to intercept MPI calls from the application and to hide all redundancy-related mechanisms. A redundantly executed application runs with r * m MPI processes, where m is the number of MPI ranks visible to the application and r is the replication degree. RedMPI supports partial replication, e.g., a degree of 2.5 instead of just 2 or 3, for tunable resilience. It also supports a variety of message-based replication protocols with different consistency. Not counting in the need for additional resources for redundancy, results show that the most efficient consistency protocol can successfully protect HPC applications even from high SDC rates with runtime overheads from 0 % to 30 %, compared to unprotected applications without redundancy. Partial and full redundancy can also be combined with checkpoint/restart [85]. Non-linear trade-offs between different levels of redundancy can be observed when additionally using checkpoint/restart, since computation on non or less redundant resources is significantly less reliable than computation on fully or more redundant resources.
Interplay between resilience and dynamic load balancing Scheduling of application jobs at the system level contributes to exploiting parallelism by placing and (dynamically) balancing the batch jobs on the local site resources. The jobs within a batch are already heterogeneous; yet, current batch schedulers rarely co-allocate, and most often only allocate, computing resources (while network and storage continue to be used as shared resources). Dynamic systemlevel parallelism can arise when certain nodes become unavailable (due to hard and permanent errors) or recover (following a repair operation). This can be exploited during execution by increasing opportunities for system-level co-scheduling in close proximity of jobs that exhibit different characteristics (e.g., co-scheduling a classical computeintensive job in close proximity to a data-intensive job) and by dynamic resource reallocation to jobs that have lost resources due to failures or to waiting jobs in the queue.

Detected errors mitigated with assistance
In this section we focus on correction methods that need assistance from the upper layers in order to achieve resilience and correctness. It is important to note that there are multiple methods that offer assisted fault tolerance but some of them involve a few additional lines of code while others require rewriting the whole applications using a specific programming model. Therefore, we will divide this section into subsections depending on the programming and/or redesign effort that is required.

Correction with incremental redesign
As explained in Section 2.1, it is possible to perform system-level checkpointing without any feedback from the application or the algorithm or any upper layer. The issue with system-level checkpointing is that the size (and therefore the time and energy cost) of checkpointing is much larger than what is really required to perform a restart of the application. Thus, application-level checkpointing is an attempt to minimize the size of checkpoints to the minimum required for the application to be able to restart.
Performance modeling and optimization of checkpoint-restart methods Research on simulation tools assessing the performance of certain checkpoint-restart strategies is presented in various publications [15,73,87,165]. Different theoretical approaches are used and tools are developed that either simulate a fictional software or wrap an actual application.
A lot of work has been done to examine and model the performance of multilevel checkpointing approaches [28,34,156,274]. Here, the parallel distribution of the snapshots as well as the target storage system are considered as objectives for performance optimization. Asynchronous techniques are considered, such as non-blocking checkpointing, where a subset of processes are dedicated to manage the creation and reconstruction of snapshots [69,219]. As a measure to saving storage and speeding up I/O, data compression is another subject that is considered in the literature as, e.g., by Di and Cappello [74], and in one of the case studies in Section 3.2.2.
Resilient checkpointing has been considered with the help of nonvolatile memory, as for instance implemented in PapyrusKV [154], a resilient key-value blob-storage. Other resilient checkpointing techniques include the selfcheckpoint technique [245], which reduces common redundancies while writing checkpoints, or techniques reducing the amount of required memory through hierarchical checkpointing [182], or differential checkpointing [153].
Message logging Message logging is a mechanism to log communication messages in order to allow partial restart as for example examined by Cantwell et al. [51]. While improving on basic checkpointing strategies, message loggingbased approaches can themselves entail large overheads because of log sizes. The checkpointing protocol developed by Ropars et al. [213] does not require synchronization between replaying processes during recovery and limits the size of log messages. Other approaches combine task-level checkpointing and message logging with system-wide checkpointing [237]. This protocol features local message logging and only requires the restart of failing tasks. It is also possible to combine message logging with local rollback and User Level Failure Mitigation (ULFM) (Section 2.2.2) to improve log size [173].
Multilevel checkpointing libraries Current HPC systems have deep storage hierarchies involving High-Bandwidth Memory, Dynamic Random-Access Memory, Non-Volatile Memory, Solid-State Drives and the PFS, among others. Multilevel Checkpointing libraries offer a way to leverage the different storage layers in the system through a simple interface. The objective is to abstract the storage hierarchy to the user, so that one does not need to manually take care of where the data is stored or the multiple data movements required between storage levels. Each level of checkpointing provides a different trade-off between performance and resilience, where usually lower levels use close storage that offers higher performance but limited resilience, and higher levels rely on stable storage (e.g., PFS), which is more resilient but slower. Mature examples of multilevel checkpoint libraries are SCR [183], FTI [28], CRAFT [226] and VeloC [189]. Both SCR and FTI provide support via simple interfaces for storing application checkpoint data on multiple levels of storage, including RAM disk, burst buffers, and the parallel file system. Both SCR and FTI provide redundancy mechanisms to protect checkpoint data when it is located on unreliable storage and can asynchronously transfer checkpoint data to the parallel file system in the background while the application continues its execution. In addition, FTI also supports transparent GPU checkpointing. Finally, VeloC is a merge of the interfaces of both FTI and SRC. Note that some of these libraries offer the option for keeping multiple checkpoints so that the application can roll-back to different points in the past if necessary.
Containment Domains Containment Domains (CDs) provide a programming construct to facilitate the preservation-restoration model, including nesting control constructs, and durable storage [248]. The following features are attractive for large-scale parallel applications. First, CDs respect the deep machine and application hierarchies expected in exascale systems. Second, CDs allow software to preserve and restore states selectively within the storage hierarchy to support local recovery. This enables preservation to exploit locality of storage, rather than requiring every process to recover from an error, and limits the scope of recovery to only the affected processors. Third, since CDs nest, they are composable. Errors can be completely encapsulated, or escalated to calling routines through a well-defined interface. We can easily implement hybrid algorithms that combine both preservation-restoration and data encoding.
Use cases include an implementation of a parallel resilient hierarchical matrix multiplication algorithm using a combination of ABFT (for error detection) and CDs (for error recovery) [16]. It was demonstrated that the overhead for error checking and data preservation using the CDs library is exceptionally small and encourages the use of frequent, fine-grained error checking when using algorithm based fault tolerance.
Application versioning Global View Resilience (GVR) [64] accommodates APIs to enable multiple versioning of global arrays for the single program, multiple data programming model. The core idea is the fact that naive data redundancy approaches potentially store wrong applications states due to the large latency associated with error detection and notification. In addition to multiple versioning, GVR provides a signaling mechanism that triggers the correction of application states based on user-defined application error conditions. Use cases include an implementation of resilient Krylov subspace solvers [275].
Mitigating performance penalties due to resilience via dynamic load balancing Detected and corrected errors induce variation in the execution progress of applications when compared to error-free executions. This can manifest itself as load imbalance. Many application-level load balancing solutions have been proposed over the years and can help to address this problem. We mention here a few available packages.
Available load balancing software includes Zoltan [252] that requires users to describe the workload across processes as a graph and offers an object oriented interface. Further we mention Dynamic Loop Scheduling for Load Balancing (DLS4LB) [54], a recently developed library for MPI applications that contains a portfolio of self-scheduling based algorithms for load balancing. StarPU [254] proposes support for asynchronous load-balancing [172] for task-based applications. The principle is to let the application submit only a part of its task graph, let some of it execute on the platform and observe the resulting computation balance. A new workload distribution can then be computed and the application is allowed to submit more of the task graph, whose execution can be observed as well. OmpSs [80] is an effort to extend OpenMP in order to support asynchronous execution of tasks including a transparent interface for hardware accelerators such as GPUs and FPGAs. OmpSs is built on top of the Mercurium compiler [251] and the nanos++ runtime system [249].
HCLib [271] is a task-based programming model that implements locality-aware runtime and work-stealing. It offers a C and C++ interface and can be coupled with inter-process communication models, such as MPI. Charm++ [151] features an automatic hierarchical dynamic load balancing method that overcomes the scalability limitation of centralized load balancing as well as the poor performance of completely distributed systems. Such a technique can be triggered dynamically after a failure hits the system and the workload needs to be redistributed across workers.

Correction with major redesign
The correction of some detected errors might have a strong impact of the algorithm that has to implement the mitigation. The mitigation design can be made more affordable if some components of the software stack have already some appealing features to handle such situations.
Resilience support in the Message Passing Interface (MPI) Most MPI implementations by default are designed to terminate all processes when errors are detected. However, this termination occurs irrespective of the scope of the error, requiring global shut-down and restart even for local errors in a single process. This inherent scalability issue can be mitigated if MPI keeps all survived processes to continue and/or if restart overheads are reduced. The MPI community has proposed several recovery approaches, such as FA-MPI [127] or MPI-ULFM [41] to enable alternatives of global shut-down, as well as better error handling extensions, like MPI Reinit [159], to reduce overhead and impact of failures. Among these approaches, MPI-ULFM is the most advanced and well known. It provides a flexible low-level API that allows application specific recovery via new error handling approaches and dynamic MPI communicator modification under process failures, although with significant complexities for the application developer using the new APIs. Several approaches have been proposed to mitigate this complexity by creating another set of library APIs built atop of MPI-ULFM [51,99,100,225,253]. However, as of now, in part due to its complexity when used on real-world applications and limited support in system software, MPI-ULFM as a whole has not been adopted in the MPI standard and hence is not readily usable for typical HPC application programmers. Nevertheless, various aspects of ULFM are in the process of standardization and will provide more mechanisms in MPI to build at least certain fault tolerant applications, starting with the upcoming MPI 4.0 standard.
Resilience abstractions for data-parallel loops Data-parallel loops are widely encountered in N -body simulations, computational fluid dynamics, particle hydrodynamics, etc. Optimizing the execution and performance of such loops has been the focus of a large body of work involving dynamic scheduling and load balancing. Maintaining the performance of applications with data-parallel loops running in computing environments prone to errors and failures is a major challenge. Most self-scheduling approaches do not consider fault-tolerance or depend on error and failure detection and react by rescheduling failed loop iterations (also referred to as tasks). A study of resilience in self-scheduling of data-parallel loops has been performed using SimGrid-based simulations of highly unpredictable execution conditions involving various problem sizes, system sizes, and application and systemic characteristics (namely, permanent node failures), that result in load imbalance [241]. Upon detecting a failed node, re-execution is employed to reschedule the loop iterations assigned to the failed node.
A robust Dynamic Load Balancing (rDLB) approach has recently been proposed for the robust self-scheduling of independent tasks [180]. The rDLB approach proactively and selectively duplicates the execution of assigned chunks of loop iterations and does not depend on failure or perturbation detection. For exponentially distributed permanent node failures, a theoretical analysis shows that rDLB is linearly scalable and its cost decreases quadratically with increasing system size. The reason is that increasing the number of processors increases the opportunities for selectively and proactively duplicating loop iterations to achieve resilience. rDLB is integrated into a dynamic loop scheduling library (DLS4LB, see Section 2.2.1) for MPI applications. rDLB enables the tolerance of up to (P − 1) process failures, where P is the number of processes executing an application. For execution environments with performance-related fluctuations, rDLB boosts the robustness of Dynamic Loop Self-scheduling (DLS) techniques by a factor up to 30 and decreases application execution time up to 7 times compared to their counterparts without rDLB.
Resilience extension for performance portable programming abstractions With the increasing diversity of the node architecture of HPC systems, performance portability has become an important property to support a variety of computing platforms with the same source code while achieving a comparative performance to those programmed with the platform specific programming models. Today, Kokkos [82] and Raja [30,250] accommodate modern C++ APIs to permit an abstraction of data allocation and parallel loop execution for a variety of runtime software and node architectures. This idea can be extended to express the redundancy of data and computation to achieve resilience while hiding the details of the data persistence and redundant computation. Recently, the resilient version of Kokkos was proposed for a natural API extension of Kokkos' data (memory space) and parallel loop (execution space) abstractions to (1) enable resilience with minimal code refactoring for the applications already written with Kokkos and (2) provide common interface to call any external resilience libraries such as VeloC [189]. The new software will be released in a special branch in https://github.com/kokkos/kokkos.
The resilience abstraction idea has also been applied to task parallel programming models such as Charm++ [151], HClib [271], HPX [150], OmpSs [80] and StarPU [254] to integrate a variety of resilient task program execution options such as replay, replication, algorithm-based fault tolerance and task-based checkpointing. Task-based programming models indeed have a very rich view over the structure of the application computation, and notably its data, and have a lot of control over the computation execution, without any need for intervention from the application. Replaying a failed task consists of issuing it again with the same input, discarding the previous erroneous output, and replicating a task consists of issuing it several times with different output buffers and comparing the result. Dynamic runtime systems can then seamlessly introduce replay and replication heuristics, such as trying to run different implementations and/or computation units, without the application having to be involved beyond optionally providing different implementations to be tried for the same task.
The task graph view also allows for very optimized checkpointing [171,172,254]. In the task-based programming model, each checkpoint is a cut in the task graph, which can be expressed trivially within the task submission code, and only the data of the crossing edges need to be saved. Even better, the synchronization between the management of checkpoint data and application execution can be greatly relaxed. The transfer of the data to the checkpoint storage can indeed be started as soon as the data is produced within the task graph, and not only once all tasks before the checkpoint are complete. A checkpoint is then considered complete when all its pieces of data have been collected. It is possible that tasks occurring after the checkpoint may run to completion before the checkpoint itself is completed. All in all, this allows for a lot more time for the data transfers to complete, and lessens the I/O bandwidth pressure.
Software engineering approaches for resilience by design Resilience design patterns [141,142] offer an approach for improving resilience in extreme-scale HPC systems. Frequently used in computer engineering, design patterns identify problems and provide generalized solutions through reusable templates. Reusable programming templates of these patterns can offer resilience portability across different HPC system architectures and permit design space exploration and adaptation to different (performance, resilience, and power consumption) design trade-offs. An early prototype [14] offers multi-resilience for detection, containment and mitigation of silent data corruption and MPI process failures.

Numerical algorithms for resilience
In this section, we focus on the handling of errors at the algorithmic level. We see three different classes of problems to tackle here: (i) detection of un-signaled errors (mostly bit flips and other instances of silent data corruption, Section 3.1), (ii) correction of errors that have been signaled but could not be corrected at the hardware or middleware layer (by error aware algorithms, Section 3.2), (iii) design of error oblivious algorithms that deliver the correct result even in the presence of (not too frequent) errors (Section 3.3).
In addition to correctness in the presence of errors, an important challenge in all our considerations is efficiency in terms of algorithm runtime. In this context, additional algorithmic components such as work stealing and asynchronous methods (where missing data are simply an extreme case of delay) have to be considered. We mention these methods when describing methods that can make use of such runtime optimizing measures.

Error detecting algorithms
In this section, we focus on mechanisms to numerically detect errors that have not been detected by the underlying system or middleware. We have identified several techniques that allow us to (likely) notice the occurrence of an error at several layers of numerical algorithms. Table 1 gives an overview of some detection techniques and the algorithmic components or numerical methods where they are applicable.

Exceptions
Exceptions are a way a program signals that something went wrong during execution. We consider the case where exceptions are caused by data corruption that can, for example, lead to division by zero or out-of-range access. Most programming languages support a way of handling exceptions. The algorithm programmer can register an exception handler that gets called whenever an exception occurs. If the error is recoverable, the exception handler will specify how best to continue afterwards. If the error is not recoverable, the program will be aborted. Exceptions are a straight- forward way to detect certain types of errors and can be applied to all numerical algorithms. However, they obviously only see a small subset of all possible errors and it is not trivial to decide when to use exceptions handlers in the light of a trade-off between correctness, robustness and runtime efficiency.

Checksums
Checksums could be used at the hardware or middleware layer to detect errors, but here we will discuss checksums as employed on the algorithmic layer where we have a more detailed knowledge about the existence of numerical or algorithmic invariants. Checksum techniques have been used in various numerical algorithms. We list some examples below.
BLAS 2/3: Checksum encoding matrices, introduced by Huang and Abraham [137] requires (i) adding redundant data in some form (encoding), (ii) redesign of the algorithm to operate on the respective data structures (processing), and (iii) checking the encoded data for errors (detection). We ignore the recovery phase here and refer to Section 3.2.
Checksums are used in FT-ScaLAPACK [267] for dense matrix operations such as MM, LU and QR factorization and more recently in hierarchical matrix multiplication [16]. Wu et al. give a good survey of checksum deployment in dense linear algebra [268].
Gauss-Seidel/Jacobi/SOR and multigrid: In [179], checksums are used to detect errors in the Jacobi smoother, the restriction and interpolation operators of a multigrid method solving a two-dimensional Poisson equation.
Krylov subspace methods: Tao et al. propose a new checksum scheme using multiple checksum vectors for sparse matrix-vector multiplication, which is shown to be generally effective for several preconditioned Krylov iterative algorithms [247]. Also [1,227] use checksums for protection within the conjugate gradient (CG) algorithm.
FFT: Checksum can also be used in Fast Fourier Transforms (FFT)s similarly as in matrix-vector multiplication. Liang et al. [167] develop a new hierarchical checksum scheme by exploiting the special discrete Fourier transform matrix structure and employ special checksum vectors. Checksums are applicable to many important kernels such as matrix-matrix multiplication, but are costly. In addition, it can be difficult to specify a suitable threshold for 'equality' in the presence of round-off errors. For many numerical calculations such as scalar products, checksums are not applicable at all.

Constraints
In some applications, constraints for different types of variables are known. Examples are positivity constraints, conservation laws for physical quantities or known bounds for internal numerical variables.
Krylov subspace methods: Resilience was already of importance in the early days of digital computers. In the original PCG paper [132], Hestenes and Stiefel noticed that the reciprocal value of α (the step length) is bounded above (repectively, below) by the reciprocal of the smallest eigenvalues (respectively the inverse of the largest eigenvalue) of the matrix. The inequality involving the largest eigenvalue (for which in practice it may be cheaper to get an approximation) was used to equip PCG with error detection capabilities in [1].
Partial differential equations: Checking for bounds can be associated with minimal or extremely high cost depending on whether extra information has to be computed (such as eigenvalues of matrices) or not. Reliability is, in general, an issue as only those errors leading to violation of these constraints can be detected. An example of the use of problem-informed constraints can be found in [186]. In this work, the authors derive a priori bounds for the discrete solution of second-order elliptic PDEs in a domain decomposition setting. Specifically, they show that the bounds take into account the boundary conditions, are cheap to compute, general enough to apply to a wide variety of numerical methods such as finite elements or finite differences, and provide an effective way to handle faulty solutions synthetically generated.

Technical error information
In many numerical large scale applications, the main computational task involves the approximate computation of integrals, algebraic systems, systems of ODEs or PDEs. For all these problems, various types of error information such as residuals, differences between iterations, round-off error estimates and discretization error estimates can be used as indicators of errors either by their size or by monotonicity criteria. We give several examples from literature for different classes of numerical algorithms.
Krylov subspace methods: Round-off error bounds can be used in Krylov subspace methods. They fit in the general framework of round-off error analysis [133] and have been considered in the context of Krylov subspace methods in finite precision arithmetic [169,178].
Vorst and Ye proposed a residual gap bound [256] (bound for the norm of the residual gap between the true and the computed residuals) based on round-off error analysis that was later used as a criterion for actual error detection in [1] when bit flips occur. The detection of errors in Krylov methods via violation of orthogonality is proposed in [63].
Multigrid: Calhoun et. al [50] apply a residual/energy norm-based error detection for algebraic multigrid. They use two criteria: (i) the reduction of the residual norm as a weak criterion and (ii) the reduction of the quadratic form when solving the linear system Ax = b for symmetric positive matrices.
The quadratic for E calculated at level i during the down-pass of a V-cycle should be less than the energy calculated at level i during the down-pass of the next V-cycle.
When using the full approximation scheme residual norm reductions can also be verified at each level in the hierarchy of a multigrid-cycle. The structure of the full approximation scheme additionally provides smart recovery techniques utilizing its lower resolution approximations [11].
Time-stepping: For iterative time-stepping with spectral deferred corrections, monitoring the residual of the iteration can be used to detect errors in the solution vectors [119]. In the context of parallel-in-time integration with parareal, consecutive iterates are considered in [191] to detect errors in the solution vector. In [35], an auxiliary checking scheme in contrast to the original base scheme is used to detect and correct errors during implicit and explicit time-integration. Estimating the local truncation error with two different methods is used in [121] to implement a resilient, high-order Runge-Kutta method. This "Hot Rod" approach is then also used for error correction.

Multi-resolution
Multi-resolution means that information is available at different resolution levels, in terms of spatial discretization (PDE), time discretization (ODE and PDE), order of discretization (PDE in space and time), matrix dimensions (numerical linear algebra, multigrid), frequencies, and so on. This leads to a certain redundancy -not an artificially introduced, but an inherently available one. This redundancy can be used to detect discrepancies or anomalies and, hence, errors that could not be detected by the system. There are numerous examples for the mentioned problem classes, we outline one example in more detail here.
Sparse grids / Combination technique: Sparse grids [48] are one particular class of multi-resolution methods. There, via the use of hierarchical bases, certain structures often seen in d-dimensional data can be exploited to alleviate the curse of dimensionality, without a significant loss of accuracy. Sparse grids have been successfully used in a wide range of problem classes where spatial discretization plays a role, such as interpolation [145], quadrature [47,109,110], solvers for PDEs [124,129], or machine learning tasks [104,105,201] (e.g., classification, regression, clustering, or density estimation). One particular incarnation of sparse grid methods is the so-called combination technique [117]. There, based on an extrapolation-style approach, a linear combination of a specific set of full, but very coarse-grid solutions is used to get a sparse fine-grid solution. The various coarse grid solutions can be obtained in a completely independent way, using (parallel) standard solvers. This opens the way to (1) a natural two-level parallelization and to (2) an easy and cheap detection of system undetected errors: Since we actually compute solutions for the same problem on different (i.e., differently discretized) grids anyway, we can use these to detect anomalies -just by comparing the available solutions. And the detection leads immediately to a mitigation strategy (see Section 3.2.2), since we can easily exchange coarse grids in case of errors, just by changing the combination pattern [8,9,123,130,134,192]. Therefore, this is an example for a smart algorithm that is able to do both detection and mitigation.
Further examples are mentioned in Section 3.1.4 as multi-resolution typically comes with corresponding error estimates based on differences between solutions at different resolution levels: multigrid and parallel time stepping.

Redundancy
Redundancy is a strategy for error detection that can be applied to all of the numerical algorithms mentioned in Table 1. It covers two approaches. In the first approach computational resources may be replicated twice or thrice. Such instances are called DMR [144,265] or TMR [223,261]. In the second approach the computations are repeated twice or thrice on the same resource [17,259]. An advantage of this approach is the flexibility at the application level.
Note that the first approach costs more in space or resources, the second approach costs more in time.
The redundancy based error detection technique described in [33] relies on in-depth analysis of application and platform dependent parameters (such as the number of processors and checkpointing time) to formalise the process of both resource and computation replication. It provides a closed-form formula for optimal period size, resource usage and overall efficiency.
Ainsworth et. al [5] use replication of fault-prone components as an error detection technique in a multigrid method. Also error detection in the time stepping methods from [35] mentioned in Section 3.1.4 can be interpreted as redundancy based error detection.
The main disadvantage of replication is its cost in terms of performance, although recomputing only some instructions instead of the whole application lowers the time redundancy overhead [193]. However, redundancy in some calculations should in particular be considered as a possible strategy for error detection as in modern supercomputers the cost of arithmetic operations tends to decrease compared to communication time.

Error aware algorithms
In this section, we look at error correction techniques within an application. We assume that the application has been notified that part of the algorithm's data is corrupted or lost. In that context, mitigation or containment actions have to be undertaken at the algorithmic design level, where the appropriate actions depend on the data detection granularity and how the notification mechanism was activated. It is possible to design both lossy and lossless mitigation procedures that are tailored to the numerical algorithms under consideration.
In Section 3.2.1 we give a brief literature overview of ideas that can be used to complement numerical mitigation or containment procedures. Then, in Section 3.2.2 we offer a more detailed discussion of some recent successful attempts by presenting a few case studies in the context of the solution of Partial Differential Equations (PDE).

Error aware algorithms for the solution of linear systems
A wealth of literature already exists on various, mostly isolated ideas and approaches that have appeared over time. Checkpoint-restart methods are the most generic approaches towards resilience for a broad spectrum of applications, see Section 2.2.1 for an introduction. We first describe a general mental model to design resilient numerical algorithms independent of actual machine specifications that lead to what is nowadays referred to as Local-Failure Local-Recovery (LFLR) techniques. Then we move to 'classical' algorithm-based fault tolerance, which originally was developed to detect and correct single bit flips on systolic architectures devoted to basic matrix computations, see Section 3.1.2. Finally, we discuss a range of ideas and techniques not covered by the case studies below.
Local-failure local-recovery As far back as a decade ago, an abstract framework was developed to separate algorithm design from unclear machine specifications, see also Section 2.2.1. The idea of a selective reliability model as introduced by Hoemmen [45,135] is machine-oblivious and highly suitable for algorithm design for machines with different levels of (memory) reliability. It has led to the concept of Local-Failure Local-Recovery (LFLR) [253]. This model provides application developers with the ability to recover locally and continue application execution when a process is lost. In [253], Teranishi and Heroux have implemented this framework on top of MPI-ULFM (Section 2.2.2) and analyzed its performance when a failure occurs during the solution of a linear system of equations.
Original algorithm-based fault tolerance with checksums The term Algorithm-Based Fault Tolerance (ABFT) was originally coined in conjunction with protecting matrix operations with checksums to handle bit flips [136], mostly assuming exact arithmetic calculation for detection and mitigation. (See Section 3.1.2 for a more detailed discussion on checksums). The main drawback of checksums is that only limited error patterns can be corrected and its robust practical implementation in finite precision arithmetic can be complicated to tune to account for roundoff errors. A second drawback is that the checksum encoding, detection and recovery methods are specific to a particular calculation. A new scheme needs to be designed and proved mathematically for each new operation. A further drawback is to tolerate more errors, more encoded data is needed, which may be costly both in memory and in computing time.
ABFT concepts have been extended to process failures for a wide range of matrix operations both for detection and mitigation purposes [44,62,79,147,269] and general communication patterns [149]. ABFT has also recently been proposed for parallel stencil-based operations to accurately detect and correct silent data corruptions [58]. In these scenarios the general strategy is a combination of checkpointing and replication of checksums. In-memory checkpointing [147] can be used to improve the performance. The main advantage of these methods is their low overhead and high scalability.
In practice, the significance of a bit flip strongly depends on its location, i.e., which bit in the floating point representation is affected. Classical ABFT has been extended to take into account floating point effects in the fault detection (checksums in finite precision) as well as in the fault correction and to recover from undetected errors (bit flips) in all positions without additional overhead [181].
Iterative linear solvers Iterative linear solvers based on fixed point iteration schemes are, in general, examples of error oblivious algorithms, as described in Section 3.3. The convergence history of the scaled residual norm observed within the iterative scheme often resembles the curves displayed in Figure 2. In this case the iterative scheme is a multigrid method, as in [115,138]. The peaks in the residual occur after data has been lost and when the iterations are allowed to restart with some form of replacement of the lost data. In the simplest case, the lost data may just be re-initialized with the value of zero, and recovery techniques to obtain better solutions are discussed in Section 3.2.2.
It can be seen that, depending on when in the course of the iteration a small portion of the approximate solution suffers from an error, we observe a delay in convergence, directly proportional to an increase in runtime. In the case where errors appear too often, the solver might not recover and other mitigation actions might have to be considered.  [195]. In contrast to restarting, a number of algorithm based recovery strategies have been proposed, including approximate or heuristic interpolation methods [2]. An approach of exactly recovering the state of the iterative solver before the node failure has been investigated for the Preconditioned Conjugate Gradient (PCG) and related methods [164,196]. This also includes studying scenarios with multiple simultaneous node failures [197] and scenarios where no replacement nodes are available [194].
Approximated recovery and restart in sparse numerical linear algebra For matrix computations, eigensolvers or basic kernels such as iterative linear system solvers, some recovery ideas rely on forming a small dimensional linear algebra problem where the inputs are the still valid data and the unknowns are the lost/corrupted ones. The outcome of this procedure is subsequently used to replace the lost/corrupted data and the numerical algorithm is somehow started again from that meaningful initial guess. The recovery procedure is tailored to the actual numerical algorithm. As an example, consider a fixed point iteration scheme for a linear system and suppose the lost data are entries of the iterate vector, the most dynamically evolving data in this computational scheme. Matrix entries of the iteration scheme related to the lost data, as well as some neighbouring entries, serve to build the left-hand side of a linear problem (either a linear system or a least-square problem) while the right-hand side is built from valid data. The solution of this small problem is then used to replace the corresponding lost entries of the iterate vector. The complete, updated vector is taken as a new initial guess when restarting the fixed point iteration. If the data is not corrupted too often the classical convergence theory still applies and because the new initial guess incorporates updates from the calculations performed before the error was detected, the global convergence rate is not strongly affected. The method described in adaptive recovery techniques for extreme scale multigrid in Section 3.2.2 is an example application of this technique.
For numerical schemes based on nested subspace search, such as Krylov subspace methods, closely related techniques have been successfully applied both for eigensolvers and linear solvers that further exploit the sparsity structure of the matrices to reduce the computational cost associated with the recovery procedure. At the cost of a light checkpoint performed once when starting the linear solver (mostly the matrix and the right-hand side vector in case of linear system solution) this mitigation approach has no overhead if the data is not corrupted during the solution computation. We refer to [2,3,161] for some illustrations on those numerical remedies in a parallel distributed memory framework and to [146] where these ideas are exploited for a lower granularity of data loss in a task-based runtime system. See Section 2 for references relevant to task-based runtime systems.
We also note that these ideas can be extended to hybrid iterative/direct numerical schemes, that have a domain decomposition flavor, where the recovery procedure can be enriched with additional features of the parallel numerical scheme such as redundancy or properties of the preconditioners [4]. They can also be extended to the time domain in the context of multilevel parallel-in-time integration techniques [229].

Error aware algorithms for the solution of partial differential equations
The ideas introduced above in Section 3.2.1 are application agnostic but naturally apply to linear systems arising from the discretization of a PDE. In that latter case, more information from the underlying PDE can be closely tailored to intrinsic features of solvers such as multigrid. In this section we discuss some research works on mitigation and containment that exploit the properties of PDEs to aid the recovery techniques. We also present some mitigation processes that are only relevant in the PDE setting.
Adaptive recovery techniques for extreme scale multigrid Some of the most efficient solvers of PDE, such as parallel geometric multigrid methods [114,143], can be based on the exchange of ghost layers in a non-overlapping domain partitioning. This automatically leads to a redundancy in interface data between subdomains that in turn permits the design of an efficient two-step recovery strategy for iterative solvers. This is of particular interest in largescale parallel computations. When each subdomain is large, then the ratio between the data on its surface and the volume data in its interior becomes small.
When a processor fails, the information within one or several subdomains is lost. For the recovery and continued solution, the redundant ghost layer information is used in a first step, to recover locally either Dirichlet-or Neumann-type data for the subdomains. The global problem can then be formulated in two partitions, the outer healthy subdomain and the inner faulty subdomain, where the recovery must reconstruct the lost data. Both subproblems must be bi-directionally coupled via the interface and the corresponding ghost layers of unknowns.
After re-initialization, the corrupted and reinitialized data could pollute the solution globally, meaning that the locally increased error in the faulty domain can spread globally and thus also affect the healthy subdomain. In order to avoid this pollution, the communication between the healthy and faulty sub-problems is interrupted during the second step of the recovery process. In the second step, we continue with the original iterative solver restricted to the healthy sub-problem and select a suitable one for the faulty one. After some number of asynchronous iteration steps both subproblems are reconnected, see [138], and the global iterative solver strategy is resumed. Note that the reconnecting step is mandatory for the convergence of the iterative solver. The tearing step separating the subdomains is mandatory to preserve the accuracy of the dynamic data in the healthy sub-problem, and without this step the corrupted data from the faulty sub-domain pollutes the global solution. Of critical importance for the performance of the method are the accuracy of the faulty sub-problem solver at re-connection time and the time spent in the recovery mode. In the faulty domain, the lost data can be initialized with 0, or, alternatively, compressed checkpointed data can be used as described in the following section on compression techniques for checkpoint-restart. Note, however, that with straight-forward compression techniques, compressed checkpoint data will only be useful to recover the low frequency components in the faulty domain. If the local recovery is performed with multigrid, then the low frequencies are in any case cheap to recover, as compared to the cost of recomputing the lost high frequency components.
The accuracy within a multigrid strategy can be easily controlled by a hierarchical sum of weighted residuals [216]. The overhead cost for the a-posterior error indicator is quite small compared to the overall solver cost. Having an estimate for the algebraic error in both sub-problems at hand, the re-connection step is determined automatically. To speed up the time which is spent in the recovery, a so-called 'superman strategy' is applied [138], see also Figure 3 for an illustration. More resources compared to the situation before the fault are allocated to the faulty sub-problem. A short recovery phase in combination with carefully selected re-coupling criteria then guarantees a highly efficient fault-tolerant solver.
Of special interest is a massively parallel multigrid method as base solver. In combination with the tearing and intersection approach for the recovery, it results in a hybrid approach. In case of a Stokes-type system, yielding after discretization a saddle point problem, the strategy can either be applied on the positive definite Schur complement for the pressure or, as it was done in [139], on the indefinite velocity-pressure system. In that case an all-at-once multigrid method with an Uzawa-type smoother acting on both solution components turns out to be most efficient, see [78]. Numerical and algorithmic studies including multiple faults and large-scale problems with more than 5 · 10 11 degrees of freedom and more than 245000 cores have been demonstrated [138,139]. The automatic re-coupling strategy is found to be robust with respect to the fault location and size and also handling multiple fault. In many scenarios a complete recovery can be achieved with almost no increase in runtime and while maintaining excellent parallel efficiency.
Adaptive mesh refinement, load balancing, and application level checkpointing Adaptive Mesh Refinement (AMR) functionality and load balancing require similar data linearization-and transfer functionality as is needed for application level checkpointing. This is exploited in the waLBerla framework [24,221,222] that features an object oriented design for composing coupled multiphysics simulation software. waLBerla's load balancing is based on routines to transfer simulation data between processors so that functionality to serialize, pack, send, and unpack all relevant data is already available as a by-product of the AMR functionality. Note that the waLBerla software architecture imposes this structure for Eulerian mesh based data as well as for Lagrangian particle-based models and it canonically extends to coupled Eulerian-Lagrangian multiphysics models. For this to work transparently, the routines for migrating simulation data must be suitably encapsulated. Then this functionality can be used to write user level checkpoints either on disk or in memory. Note that writing checkpoints will inevitably imply overheads in memory consumption and communication time, but that restoring a checkpoint is cheap, since it initially only requires re-activating the redundantly stored data. This is especially true when in-memory checkpointing is used as explored and analyzed in [156]. The simple restoration of checkpointed data may of course lead to load imbalance, but the functionality to redistribute load is also available as part of the parallel AMR functionality. In this sense, user-level checkpointing can be based in a natural, efficient, and straightforward way on the functionality of parallel AMR algorithms combined with load balancing functionality.
Compression techniques to accelerate checkpoint-restart for Kryloy-MG solvers Compressed checkpointing is a possibility to improve the efficiency of classical checkpoint-restart schemes, both in terms of the overhead to generate the checkpoints and to recover the data if an error occurs. The added efficiency mainly comes from a reduced memory footprint which is beneficial for communication and storage. It is particularly efficient if the compression method is tailored to the target application. As an example, in-memory compressed checkpoints combined with LFLR (see Section 3.2.1) for iterative linear solvers, e.g., multigrid preconditioners in Krylov schemes, are described below.
Lossy Compression: As already mentioned in Section 3.2.2, paragraph 'Approximated recovery and restart', initially only the dynamical data, i.e., the approximate solution, are protected. Lossy compression allows a balance between the accuracy of the discretization error of the assembled system and the numerical error within the solver. Specifically in [10], the SZ library [72,166,246] is employed, which prefers, by construction, structured data ideally associated with a structured grid. Another important feature is that the compression accuracy can be prescribed and adapted to the situation. Unfortunately, a higher compression accuracy usually leads to a lower compression rate and higher compression time, which is crucial in terms of resilience overhead.
Note that multigrid can be interpreted as a lossy compression technique in itself, with a number of mathematical peculiarities that need consideration [115]. Multigrid algorithms use a hierarchy of grids to solve linear systems in an asymptotically optimal way. This hierarchy can be used to restrict, i.e., lossily interpolate, the iterate from fine to coarse grids. Such a lower-resolution representation of the iterate can then be stored as a compressed checkpoint. Conversely, the multigrid prolongation (coarse-to-fine grid interpolation) operator is used to decompress the data. With only small additional computations, the multigrid hierarchy can also be used for error detection.
Recovery: Several recovery techniques can be devised [10]. As a baseline approach checkpoint-restart is mimicked and the global iterate is simply replaced with its decompressed representation, independently of the compression strategy. The second proposed approach follows the LFLR strategy and re-initializes only the local data that is lost on faulty computing nodes by using checkpoint data stored on neighbouring computing nodes. Contrary to the first approach, this is mostly local and only needs minimal communication to receive a remotely stored backup. In particular, the recovery procedure itself does not involve the participation of other processes except those sending the checkpointed data. As a worst-case fallback when the backup data is not sufficient, a third recovery approach is established, which is still mostly local. Here, an auxiliary problem is solved iteratively with boundary data from the neighbouring computing nodes. This is similar to the adaptive recovery techniques for extreme scale multigrid from above or the approximated recovery and restart of Section 3.2.1. An auxiliary problem is constructed, either by domain decomposition overlap or the operator structure, and solved with an initial guess based on the checkpoint data to accelerate the iterative recovery phase. Experiments show that this approach can almost always restore the convergence speed of the fault-free scenario independently of the used backup technique, only the number of additional local recovery iterations varies. For more details, we refer to [10].
Resilience with sparse grids Resilience can be added on various abstraction levels of the algorithm. For PDE problems one traditionally adds resilience on the level of linear algebra operations, on the solver level for linear/nonlinear equations, or on the time-stepping algorithm. However, this may in some cases not be coarse-grained enough to minimize the overhead of resilience techniques, especially when errors occur rarely. In [123,129,130,134,192,199] the authors demonstrate a fault-tolerant framework for solving high-dimensional PDEs that applies fault tolerance on top of the individual PDE solver. The framework boosts the scalability of black-box PDE solvers while making it simultaneously resilient to faults by applying the sparse grid combination technique. In this technique the PDE simulation is distributed over many coarse grids, which can be processed in parallel. At regular intervals the results of these grids are combined to obtain the final sparse grid result. In presence of faults the affected grids can be neglected and an alternative combination scheme is calculated via an optimization routine. If too many grids are lost, the last combination result serves as an in-memory checkpoint to recompute the required grids. In [192] it is shown that this lossy recovery provides very good results even with high error frequencies. At the same time the parallel efficiency is only slightly affected.
Adaptive mesh refinement Adaptive refinement techniques in combination with finite element methods are well established for fault-free computations. In terms of fault tolerance, this means that in addition to the assembled linear system, the geometric mesh structure must be protected. This requires the reconstruction of the data structures containing the mesh hierarchy. For the use of multigrid or multilevel methods, we also need to recover multiple levels of adaptive grid refinement after a fault has occurred. The recovery process must take into account the intra-grid as well as the inter-grid data dependencies.
We refer to [233] for a parallel adaptive multigrid method that uses a sophisticated dynamic data structures to store a nested sequence of meshes and the iterative evolving solution. Stals demonstrates that it is possible to implement a fault recovery procedure that builds on the original parallel adaptive multigrid refinement algorithm [232] in the case of a fail-stop fault. It is assumed that a copy of the coarsest grid can always be accessed after a fault has occurred, i.e., it is stored off the processor. The challenge in recovering an adaptively refined grid is that the mesh distribution changes during any load balancing procedures, i.e., the local information that was available during the original refinement process will have been modified or removed. Nevertheless it is demonstrated that the neighbouring healthy processors contain enough intact information so that the necessary structure can be recovered to pass into the refinement routine.
In the case of uniform refinement, the original multilevel grid is recovered. In the case of an adaptively refined grid, enough of the structure is recovered to re-establish the correct communication pattern allowing the solution process to run to completion, but potentially with reduced accuracy. The neighbouring healthy processors will only contain enough information to guide the refinement around the edge of the recovered subgrid. Further refinement within the interior of the recovered subgrid may be required to improve the accuracy of the solution.
These techniques were implemented with minimal disruption to the original code. An example of one the few necessary modifications is that in the original code, communication was used to ensure that the elements were refined in the appropriate order to avoid degenerate grids. In the resilient version of the the code that communication had to be removed as the refinement was restricted to the faulty processor.

Error oblivious algorithms
In this section, we give examples of algorithms that are error oblivious in the sense that they can recover without assistance from errors that do not occur too frequently. For example, many fixed point iterative solvers are able to execute to completion if, e.g., a bit flip error occurs in the solution vector. However, every error likely increases the execution time of the algorithm. We thus define two quality criteria for error oblivious algorithms and use to assess the examples in the remainder of this section: (i) correctness, and (ii) efficiency in terms of execution time.
Finding an algorithm that fulfills (i) and can also compete against error aware algorithms as described in Section 3.2 remains an open problem.
Error oblivious usually means that an error slowly 'leaves the system' during several iterative sweeps over the data. Error mitigation in error aware algorithms, on the other hand, requires specific measures to correct the error, and can only be applied when the error has been detected on a hardware, middleware or algorithmic layer, but removes the disturbance of the calculation process by the error immediately.
We do not expect the error oblivious algorithms to be impervious to all types of errors. An iterative method may be not error oblivious if the error changed the matrix entries. This concept is defined as selective reliability, see Section 3.2.1.

Gossip based methods
A potentially interesting alternative in large-scale parallel environments that does not require any explicit error detection mechanisms utilizes gossip-based methods and their inherent resilience properties. Such algorithms by nature build up redundancy in the system and can thus can efficiently recover automatically from various types of faults/errors without any need to explicitly detect them. In particular, Gansterer et al. have studied and extended the resilience of gossip-based aggregation and reduction methods [56,101,190]. Based on these building blocks, they have developed and analyzed several more complex resilient numerical algorithms, such as orthogonalization methods [101,102], eigensolvers [235], and least squares solvers [205].
While the strong resilience properties and execution-time robustness of these approaches are promising, there is a certain price in terms of basic runtime compared to classical deterministic numerical high performance algorithms. It remains to be investigated whether they can be competitive in a fault-prone, but otherwise classical system with global view and centralized control. Their competitiveness can be expected to increase significantly if some of these classical properties have to be weakened at the extreme scale.

Fixed-point methods
We view fixed-point methods as methods that converge globally when certain conditions are satisfied. For example, the Jacobi iterative schemes will converge for any initial guess if the matrix is diagonally dominant. Fixed-point based iterative methods are by design resilient to bit flips. However, the convergence delay can be significant. Anzt et al. [12,13] propose techniques improving the cost-robustness with little overhead.
A class of numerical algorithms that by design have properties attractive for resilience are asynchronous iterative methods [23,29,39,40,97,230,231,242]. In order to avoid misunderstandings, we point out that this class of methods is unrelated to the idea of asynchronous dynamic load balancing [155] as addressed in Section 2.1. Instead, asynchronous iterative methods, stemming from the concept of chaotic iterations [60], are fixed-point methods that seek the solution of a problem by independently updating subdomains -which can be subdomains in the geometric sense, subsets, or individual components of the solution approximation -according to some fixed-point linear or nonlinear iterative scheme. A particularity of the asynchronous methods is that the independent updates neither adhere to a specific update order, nor synchronize in terms of a handshake with other updates, but still converge globally in the asymptotic sense. In particular, these methods are robust with respect to some subdomains being updated at a much lower pace as each update just uses the most recent non-local information available. In that sense, asynchronous solvers can have good performance in unreliable environments where messages can be dropped or processes can become unresponsive for limited time. Also, in cases where messages are corrupted (and corruption can be detected), an asynchronous solver can simply drop such a message. In cases where processes remain unresponsive, a mechanism is still needed to recover that process and its state, but the remaining processes can continue computing unchanged. Therefore, asynchronous methods are somehow error oblivious.
With the increasing cost of global synchronizations, and the attractive properties concerning fine-grained parallelization and resilience against communication delays and errors, asynchronous methods have gained attention in particular for numerical computations [243]. Chow et al. [65,66] developed an asynchronous algorithm for generating incomplete factorizations, Coleman et al. [68] further improved this algorithm by employing measures that reduce the runtime overhead when encountering errors. More general is the idea of asynchronously updating subdomains in Schwarz decompositions. In particular asynchronous restricted additive Schwarz methods and asynchronous optimized Schwarz methods have been identified to combine algorithm-inherent resilience with scalability on pre-exascale hardware architectures [83,103,112,175,270].
Independently, asynchronous multilevel methods have been proposed and analyzed under the name Fully Adaptive Multigrid method [214]. Here the multigrid smoothing process is executed asynchronously so that it can be employed for concurrent operations on different levels of the mesh hierarchy. The iteration is executed in a Southwell style [228] and is controlled by efficient hierarchical error estimators [216]. The parallel implementation [215] will automatically correct errors. More recently, asynchronous methods have been proposed for nonlinear multi-splitting [244] and eigenvalue computations like Google's Pagerank algorithm [157]. More recently, also the idea of asynchronously solving coarse-grid error correction equations has been investigated, leading to an asynchronous multigrid algorithm [266]. While case studies reveal attractive properties, these newly developed asynchronous iterative methods (such as asynchronous multigrid) are not fixed-point iterations, and developing a convergence theory for those algorithms remains a challenge.

Krylov subspace solvers
A comprehensive overview about the use of selective reliability with Krylov methods in the presence of bit flips is given in James Elliott's PhD thesis [86]. Elliott evaluates the CG and GMRES solvers with the algebraic multigrid preconditioner, see also [10] for a more recent study. Coleman et al. [68] consider Krylov subspace solvers in combination with the incomplete ILU algorithm ParILUT. In [84] Elliot et al. investigate the effect of bit flips on the convergence of GMRES and propose strategies for minimizing the numerical impact.
The authors of [31] present a monotonicity-based fault detection and correction procedure for a Generalized Conjugate Gradient Krylov solve and perform tests with manual fault injection. While the solver manages to converge even with large amounts of corrupted data, the basic recovery procedure speeds up convergence with minimal detection and correction overhead.
In [218] the authors use a slightly different terminology and call their method numerically self-stabilizing, a term which originates in the context of distributed systems [77]. They introduce two error oblivious [77] iterative linear solvers: one for the steepest descent and one for conjugate gradient. In the latter case, they consider necessary conditions for conjugate gradient to converge. Those conditions are borrowed from non-linear conjugate gradient [276] and are maintained in a correction step (typically performed every other ten iterations). The correction step does not explicitly correct errors, but re-computes quantities such as the residual at regular intervals. Therefore, we classify these methods as error oblivious instead of error aware.

Domain decomposition
In [116] Griebel and Oswald use probabilistic analysis to model the effect of errors on the convergence of the classical overlapping Schwarz algorithm. They conclude that this method does indeed converge in the presence of errors. Glusa et al. [113] mention that asynchronous domain decomposition methods are by definition fault-tolerant. In [184,210,211], the authors discuss resiliency of a task-based domain decomposition preconditioner for elliptic PDEs. By leveraging the domain decomposition approach, the problem is reformulated as a sampling problem, followed by a regression-based solution update. The regression is formulated and implemented such that it is oblivious to corrupted samples. The authors combine this algorithmic approach with a server-client implementation based on ULFM, see Section 2.2.2. They show promising results of this approach in terms of resiliency to missing tasks, corrupted data and hardware failure.

Time stepping
In [119], iterative time-stepping using spectral deferred corrections are shown to be error oblivious at the cost of more iterations for the affected time-step. With error-estimators in place, time-integration techniques like Runge-Kutta methods will repeat the calculation of a time-step with smaller step sizes, if errors in the solution vectors are relevant [61]. This type of algorithms is resilient against errors in the solution vector of the new time step. Repeating the new time step with a reduced time step size is not the optimal measure in case of an error where repeating the step with the same time step size would be more efficient, but it leads to correct results.

Future directions
In the final section we focus on the future direction of resilient algorithms. We highlight what changes need to be made to current infrastructures to support the goals proposed by algorithm and application developers. Furthermore, we list those algorithms that are likely to come to the forefront as resiliency plays a more important role in the cost-benefit analysis of extreme scale simulations. And we mention some numerical methods that are yet to be fully explored in the context of resilient algorithms.

Systems in support of resilient algorithms
We propose that resiliency will only be obtained by a multilayered approach incorporating operating systems, file systems, communication, programming models, algorithms, applications and education. In terms of the layers covered by infrastructure, the goal is to increase systems and delivered performance while keeping the detectable errors in the upper algorithm based layers constant. We refer the reader to the recently published report by Radojkovic et al. [207] for an overview of the needs of the next generation HPC systems.

Error correcting codes
Poulos et al. [204] propose hardware ECC assistance that can pass error syndrome information through to an application and use this to fix detected errors. When an ECC hardware error occurs that results in a Detectable, but Uncorrectable Error (DUE), the ECC hardware generates a syndrome which is a byproduct of the error detection. For many ECC schemes, a syndrome that corresponds to a DUE can be used to generate a list of possible corrections, one of which is taken to be the original uncorrupted data. In this work, the authors show that this set is relatively small, meaning that the set of potential values for an application to search for their correct answer (before corruption) is also small. They also study the error value distribution and show that for certain classes of problems it can be easy to identify obviously wrong answers. For the application studied in [204], work was done to correct a hydrodynamics application using conservation laws and average of neighbor cells. This work requires changes to the hardware error reporting techniques and modification to the operating system to determine which application observed the DUE and pass it to an interrupt handler.

Improving checkpoint/restart
Independent of any additions, changes or new developments in the algorithmic or the system area, checkpoint/restart will remain a necessary component for any system. For one, no other technique can provide the needed resilience against full system outages; further, checkpoint/restart is also needed for developers to deal with limited job execution times and possible migration between systems or debugging purposes at large scale.
Improving classical checkpoint/restart for homogeneous systems Observing the necessity of checkpoint/restart makes it critical to further optimize, enhance and support efficient checkpoint/restart mechanisms-even on classical, homogeneous systems-and provide users with library based solutions for core checkpoint/restart functionality. In particular, the following avenues should be pursued to optimize checkpoint/restart.
• Use additional algorithmic techniques to be able to reduce checkpoint frequency.
• Reduce data to be written to disk by eliminating redundancy and possibly compressing checkpoint information. Note that suitable data compression will typically require user-level knowledge, suitable interfaces must be provided. • Overlap/Offload checkpoint operations to allow for asynchronous checkpoint/restart operations. • Integrate checkpoint/restart with novel programming approaches to minimize checkpointable state. • Keep the restart requirements local to the neighbour nodes of the failed node.
• Localize checkpoint data to own or localized nodes. This could be supported by local non-volatile memory, as targets for checkpoint data. While this has the potential to reduce communication, as it avoids remote data transfers, it may require additional hardware support to retrieve data from non-functional nodes, e.g., by accessing data through fast JTAG-like interfaces.
• Exploit user-level knowledge for serializing, packing, compressing data, see e.g. how existing AMR functionality [156] can be exploited for efficient checkpointing in Section. 2.2.1.
Checkpoint/Restart for heterogeneous systems In addition to classical checkpoint/restart for homogeneous systems, node-local checkpoint/restart support for heterogeneous systems will help containing error and failure propagation. Such support may be provided transparently to the application by the underlying infrastructure, such as GPU drivers or task-based environments, or exposed in the programming model, such as OpenMP Offload [76].

Scheduler and resource management
Support for resilience, especially at the workflow-level, has a direct impact on resource management in HPC systems and hence requires new developments in this area as well.
Node-level parallelism With increasing node-level parallelism, the impact of OS noise (typically caused by unpredictable interrupts) becomes even more important. Therefore, dedicated node-level resources are needed to exclusively run the OS and minimize the impact of OS noise on the multi-threaded application running on the other cores.
Adaptive system and application load balancing The batch scheduler needs to adaptively balance the system load onto the available resources, via seamless application migration. While the application needs to adapt to the capabilities of the newly allocated resources, if different from the original allocation, without incurring performance penalties. The former has typically been implemented via checkpointing and process migration [174]. The latter has typically been implemented for applications that can adjust their granularity, e.g. from finer to coarser, depending on resource availability either triggered by the application or the system [46]. When exposing and expressing parallelism in applications, in addition to accounting for and matching the multiple levels of hardware parallelism (nodes, sockets, cores), the decomposition granularity needs to be flexible to support evolvability and malleability and allow for adaptive load balancing at the application and system levels.
Adaptive resource management The batch scheduler in conjunction with the distributed runtime system employed by the application (e.g., MPI, Charm++, HPX) needs to support resources errors/failures and recover them without terminating the applications in the process. This approach should work both with rigid and moldable applications as well as with evolving and malleable applications.

Programming models with inherent resiliency support
Certain applications and algorithms may naturally be resilient against errors. This may lend them as natural candidates for asynchronous parallel execution (via asynchronous many-task programming). While this mitigates the challenges associated with bulk synchronous parallel execution, asynchronous parallel execution may influence, in the presence of silent errors, the convergence rate of the numerical algorithms and might lead to incorrect results.
Programming model and runtime support for resilience can offer transparent handling of errors and failures or can assist the application in handling them. Consistent programming model support for resilience based on realistic error/failure models is needed to properly handle such events with low overhead. Higher-level abstractions for programming resilient applications are needed to help with error/failure handling complexities and to offer reuse of concepts and codes across applications.

Future directions for the solution of partial differential equations
In this section, we focus on discretizations for linear and non-linear partial differential equations as well as solvers for the resulting discrete and sparse systems of equations. We introduce a list of algorithmic properties that we found are, or can be, contributing to the resilience of the algorithms described in Section 3. Table 2 lists these properties and indicates where we found relevant examples of how they can foster resilience for either linear or non-linear solvers or for spatial or time discretization. In the following subsections we describe these examples in more detail and highlight the several (mutually related) properties that could be of interest in the context of resilient algorithms.

Redundancy and replication
A failure that is not fixed by the system (hardware and middleware) typically results in a loss or corruption of data.
To tackle this problem, redundancy techniques can be used to detect and recover from data corruption and data loss. The performance of these algorithms is usually measured in the amount of memory and computational overhead they entail, the detection rate of errors, the rate of false-positives they achieve, and the accuracy of the recovery.
Optimizing these performance indicators should be of main concern for future algorithm design. One existing class of algorithms that apply redundancy are multiresolutional techniques such as multigrid and the sparse grid combination technique described in Section 3.2. They inherently add redundancy through the hierarchical structure. Sparse grid combination techniques calculate the same solution on different anisotropic grids. The coefficients of the combinations of the components grids can be recalculated if one or more nodes are lost due to faults. This redundancy of the component grids allows the algorithm to obtain an alternative approximation of the solution. However, if a component grid is distributed on too many nodes, then the approximation will fail if a fault occurs on any one of those nodes. Another class of algorithms add redundancy through recomputation with different models and configurations such as in ensemble or multifidelity techniques. A more straight-forward approach is to directly add redundancy through replication of certain algorithmic paths, cf. the following subsection on recalculation techniques.
Depending on the underlying architecture, replication can be a competitive option to increase detected and undetected error robustness. If computation speed significantly outpaces memory access and communication, each operation can be executed multiple times while the data is still accessible in the RAM. This can be used for redundancy-based sanity checks of low-level operations or even for checksum-like approaches.
Overlapping data in parallel algorithms can serve as a starting point for mitigation, albeit not for detection. In the case studies explored in Section 3.2.2, these are applied to elliptic PDEs, though an extension to other models should be feasible. Furthermore by even increasing the ghost layer size and thereby adding extra redundancy, other reconstruction possibilities might become possible. This could already be taken into account during the domain partitioning process.

Hierarchy and mixed precision
Hierarchical discretizations have proven to be advantageous in various respects. Related notions are multi-resolution or multi-level discretizations, but also (recursive) sub-structuring in the engineering nomenclature of the Finite Element Method (FEM). Built into the hierarchy are problem-inherent information and structures that are well-suited for modern hierarchy-based solvers. In FEM, for example, hierarchical bases carry information about both location and frequency, which leads to a special built-in redundancy that can be exploited for error detection (see Section 3.1). Therefore, from a resilience perspective, hierarchy should be a core paradigm for discretization design. This applies irrespectively of whether the hierarchical bases are formulated in the spatial (h) or the order (p) sense.
From a solver perspective, multigrid methods for elliptic and parabolic PDE problems are relevant approaches towards resilient numerical algorithms. They inherently act on different granularities, representations, scales, and levels and can be used to quantify differences between these levels. For local recovery, local multigrid methods are highly efficient, especially when they can be accelerated with the superman strategy [138]. Additionally the low-resolution duplicates can be used for some kind of approximate recovery or minimal rollback like re-application of the smoother on a specific level in a multigrid scheme. Detection of errors within multigrid is often possible due to algebraic relations or on the basis of hierarchical multi-grid-inherent error estimates [11,139,216], which hold true inside such schemes. As stated in Section 4.3.1, the inherent redundancy incorporated in these algorithms is also beneficial.
Mixed-precision arithmetics are typically used within the numerical solver parts to speed-up computations. However, the discretization can enable the flexibility to store data at varying precision. Examples for this are hierarchical approaches such as hierarchical bases, where a function value is stored as a hierarchical surplus only. As another example, the usage of wavelets in multiresolutional analysis can serve. In both cases, contributions of higher levels typically require less accuracy, as only the most significant bits contribute to the overall point values.

Error control
For many numerical methods, a wide range of classical a priori and a posteriori error estimation techniques are available, see among many others [6,22,122,152,160,206,216], which constitute the basis of many adaptive numerical algorithms.
Adaptive time discretization methods are the state of the art for ODE solvers, while, for PDE solvers, spatial adaptivity techniques are also widely used. Local time step adaptation is feasible in the framework of so called local time stepping or multirate approaches, where different components of the system can have different time step sizes, see [43,53,94,106,217,220,224], which are however still far from mainstream for most applications. For PDE solvers, local spatial adaptivity techniques are also very common [20,21], but their incorporation in operational applications is still a research topic, see e.g. [36,163,185,203,255] for developments concerning oceanography and numerical weather forecasting.
The error estimations on which all these methods rely on also constitute the basis of an error detection mechanism, since some undetected errors, like bit flips on significant floating point digits, will result in errors exceeding the allowed error tolerances. To some extent, these techniques are also examples of ABFT or error oblivious approaches, since bit flips and other silent errors occurring during the computation of the solution at the next time step or on a refined mesh could be automatically corrected by the repeated computations triggered by the error threshold violation. Furthermore, silent errors in the data at the current time or mesh level could be identified by the failure of the time step or mesh refinement to correct the error.
Combined with other ABFT strategies, adaptive discretization strategies based on error estimators can be a powerful and so far rather underrated tool for protecting a simulation from undetected errors in the solution vectors. On the other hand, error estimators should not be used as a black box for resiliency purposes. Indeed, errors can lead to severe over-resolution or, potentially, even under-resolution in space or time and the error estimators themselves could be affected by undetected errors.
As seen in Section 3.1.4, some iterative solvers for the solution of linear systems have invariants, such as monotonicity for Krylov solvers. These properties can be put to good use in devising resilience strategies, for example activating an additional restart of the Arnoldi procedure as soon as an increase in the residual norm is observed.
The idea of interval arithmetic is to compute bounds of intervals that always contain the exact result [7,158]. Probabilistic methods for rounding error estimation [71,96,98,198,258] require several executions of arithmetic operations with different perturbations or different rounding modes (for instance three executions for Discrete Stochastic Arithmetic [81]). With both approaches, the comparison of several computed results enables one to control rounding errors (or detect and mitigate actually wrong results).

Locality, asynchronicity and embarassingly parallelism
One important aspect of resilient algorithms is error confinement as global dependencies propagate errors to other processors and complicate recovery. Locality-emphasizing numerical algorithms achieve this by limiting dependencies to local areas or completely removing them. Consequently, error mitigation can be limited to a local subdomain. Typical examples for these schemes are domain decomposition, which splits the domain into several subareas, and classical discretization schemes such as finite elements, finite differences and finite volumes.
As mentioned in Section 3.2.1, domain decomposition schemes such as additive Schwarz methods, or substructuringinspired FETI [90] or also the fully adaptive multigrid method [214] are naturally asynchronous and resilient to message loss. In this context, we use the term asynchronous primarily in the sense of reducing the time synchronicity in parallel computations -from communication-avoiding schemes via a reduction of synchronization points up to vastly decoupled schemes. Using this inherent property, a failure in a subdomain would result in a message loss that does not hinder convergence in other subdomains, because a global wait for a message update and synchronization are not necessary. In addition, asynchronous methods may better adapt to heterogeneous processors and networks than their synchronous counterparts as it has been shown in the context of Grid computing [19,59]. Both the localized and asynchronous approaches, achieve their impact through a decoupling of computations. Going further in this direction leads to embarrassingly or nearly embarrassingly parallel algorithms. These represent a group of algorithms where it is relatively easy to decouple subproblems in time or space. The subproblems can therefore be calculated completely independently, and errors do not propagate to other subproblems. Examples of such methods are Monte Carlo simulations and computations with the sparse grid combination technique. Since it is expected that only a few tasks will encounter errors and the scheduling is automatically balancing the load, the overall execution time does not suffer too much. Future algorithmic design should therefore aim at increasing asynchronicity and locality to move towards embarrassingly parallel problems.

Stochastic
Stochastic methods can be superior to deterministic methods when it comes to resilience. Stochastic methods do not require the program to take a deterministic path, faulty parts can be neglected or exchanged easily by other results. A popular example are Monte Carlo methods where we sample randomly in the computation domain and can simply neglect failed samples. Ensemble methods are examples where different instances or models of a concrete problem setting are computed. Even if one of these computation fails, the ensemble computation can still return a -maybe slightly less accurate -result. Stochastic elements can therefore help the future algorithm design to reduce the dependencies on specific results of the computation. These methods, however, need to be evaluated not just by highlighting their resilience properties, but also taking into account the cost of a single run: if a single run is expensive to complete, simply discarding it might be impractical.

Iterative methods
Iterative solvers may be viewed as inherently more robust than direct solvers because they do not compute their solution using a pre-defined sequence of numerical operations as direct solvers typically do. Indeed, by their nature, they perform a sequence of operations to update and improve their current approximation. If an error is encountered during computation, the probability of deleting this error or at least its effect may be higher than in a direct solver. Especially fixed-point-based methods (domain decomposition, relaxation, . . . ) may be viewed as inherently resilient as they have the property to always converge to the correct solution independent of the initial state (global convergence). Some errors may induce a low influence on convergence speed and can thus be safely ignored. In other cases, a restart -optionally with recovery techniques -may be employed to ensure both resilience and efficiency in terms of runtime.

Low memory footprint -matrix-free
The classical approach to represent linear operators as sparse matrices produces large amounts of static data which has to be restored upon failure. Checkpoint-restart approaches feature high memory cost, naturally multiples of the storage needed for the solution vector. Algorithmic alternatives to checkpoint-restart require possibly complicated or costly re-assembly. Matrix-free methods do not represent the operators as static data in the first place. Therefore, large sparse matrix data structures do not have to be restored upon failure as they are computed on the fly anyway. Extreme-scale applications will benefit from matrix-free approaches due to their low memory footprint, also in terms of runtime, (due to high memory access cost) and higher limits for the overall problem size [25,27].
In addition to saving memory and, therewith, reducing the risk of memory corruption, matrix-free methods can also be combined with automatic code generation [162] in a stencil-based approach, i.e., for finite difference methods on uniform structured grids. In such cases, the matrix entries may be 'hard wired' into code, such as 5-point stencils for Laplace's equation. Automatic code generation provides a means to increase resiliency in the code generator or domain specific language and, thus, facilitate resilience aware software development.
For finite element methods, one can use local assembly kernels [26]. Here, the trade-off between computation and storage and, in the future, resilience is relevant in particular for higher order elements.

The final mile: towards a resilient ecosystem
The future directions described above will provide critical enhancements towards providing resilient computation for numerical simulations. Alone, however, they are insufficient, as they must be embedded in the larger ecosystem and in the efforts to make that ecosystem support such novel resilience approaches. This requires another set of crucial developments.

Tools to support resilience software development
Developers will need the right tools to support their algorithmic efforts. These tools, as they exist today, are often designed without faults and errors in mind and, therefore, do not sufficiently support the development of resilient systems. In particular, we identified three areas in which enhanced tool support for resiliency is needed: a) introspection to help track errors and failures along with their root causes, b) validation through controlled fault scenarios to enable targeted testing of new error mitigation features, and c) transformation to transparently add error and failure checks into codes.
Tools for introspection Introspection is critical to ensuring early error detection and the timely activation of correction and mitigation mechanisms throughout the various layers of the software ecosystem.
System Monitoring: Knowing about the health state of a system requires monitoring it and understanding its behavior. Future work needs to focus on scalable system monitoring, real-time analyzes of system monitoring data, and autonomous decision making on corrective actions for self-aware resilient systems. In order to gain a deeper understanding, types of monitored data should be homogenized across system and sites, and, if possible, sanitized logs should be available to the community.
Application and Data Structures Monitoring: Applications need to automatically monitor their performance and correctness with the use of tools. The tools can be developed in abstraction, at the compiler-level, or at the runtime-level.
Tools for validation Currently, there are no standard tools to test the correctness and performance of resilient algorithms under undetected errors and fault. This is due to a lack of fault injection tools that reflect realistic situations. DeBardeleben et al. [120] have developed a hardware error simulator tool to understand the behavior of numerical algorithms under faulty hardware with a great accuracy, but this approach cannot evaluate the execution time of resilient algorithms at scale. Vendors provide fault injection tools [126,148] for better execution efficiency, compromising the accuracy of the hardware behavior. Compiler approaches or other in-house error injections [49,107] could allow the program to execute as efficiently as the original binary, but the correctness is further compromised. There are also tools that can analyze an application's vulnerability very quickly but do not actually produce the application's faulty output. One technique for this, DiSCvar [177], uses algorithmic differentiation and exposes how changes to each variable impact output results. It is important to note that these techniques do not actually produce that corrupted output. Hence, they are very fast but they may not be useful to developers looking to explore precisely how corruption changes their application. It is likely that a combination of these techniques, which identify most critical regions of an application coupled with fault injection at those locations, may serve as a good compromise between the two techniques.
Any novel approaches that fill the gap between the accuracy and execution efficiency of error injections will facilitate the code development of resilient algorithms, and the new tools should be built with the existing continuous integration infrastructure. Such tools likely require hardware knowledge that is considered intellectual property by the semiconductor vendors. However, efforts which explore this space using open hardware technologies (RISC-V, Sparc, etc.) can shed light on this space but may be of varying usefulness when application developers look to understand how their applications will perform on hardware that has not been fault injected at the register transfer or microcode level.
Tools for code transformation Compilers are able to generate binaries with resilience capability as suggested in the work by [209]; the generated binary instruments redundant computation, register allocations to enable error detection and correction during program execution. The recent work by Lin [170] leverages LLVM to generate SIMD instructions to perform redundant computation and verification. Source-to-source code transformation has been proposed to enable triple modular redundancy in loops [168] and automatic instrumentation of checkpointing [212]. Similarly, this idea can be extended to redundant threading for error mitigation, facilitated with OpenMP-like programming language extension [140]. These approaches automatically introduce resilience with some performance penalty, preventing the users from selective adaptation of resilience for performance optimization, and these redundant computations are benefited from the memory hierarchy, preventing doubling (or tripling) of the execution time.
In addition to such specific systems that support the addition of resilience to existing codes, automated generation of code, e.g., via Domain Specific Languages (DSL) can help with the transparent support of resilient computation. Examples for this can be stencil generators, as already discussed in Section 4.3.7.

User/Programmer education
According to the system log study by [75], many application job failures are triggered by the mistakes of the users such as script errors and program bugs including excessive file and thread creations. This means that better software engineering practices and training of users should be pursued with similar efforts to the deployment of resilience strategies.
The Exascale Computing Project (ECP) by the US DOE has made a substantial investment on educating tools, software engineering and HPC system usage for a variety of the users. Additionally, the scientific and mathematical library teams in the ECP have introduced software engineering policies [272] to improve the software quality, documentation and testing process for better interoperability and composability of multiple library packages. This activity, though not directly relevant to resilience, will gradually help to reduce application errors and failures for large scale HPC systems.

Conclusions
This article presents a snapshot of current research on resilience for extreme scale computing. It has grown out of the Dagstuhl seminar 20101 held March 1-6, 2020, bringing experts from the field together on the topic Resiliency in Numerical Algorithm Design for Extreme Scale Simulations. This seminar became a starting point to develop a synthesis between the system perspective on resilience and the algorithmic perspective.
While resilience is undoubtedly an issue for extreme scale computing, it is less clear what algorithms on the user or application level can contribute to mitigate faults. The seminar provided ample room to discuss these topics and thus became the starting point for this article. Many diverse aspects were found to be relevant, that require a holistic and multidisciplinary approach involving different and complementary scientific communities.
In particular, it clearly appeared that a fundamental distinction lies in whether faults are detected or not, and if they are not automatically detected, whether they are detectable. If they are, algorithms can often be developed to detect errors and in a second stage to correct them. It was found that some algorithms are naturally tolerant against faults or have the intrinsic feature to be error oblivious. They can thus be naturally applied on a system subject to errors.
Besides redundancy and checkpointing as classical techniques to mitigate faults, new algorithm-based resilience techniques have been developed for several classes of numerical algorithms. This includes linear algebra and solvers for partial differential equations, two classes of algorithms that are prominent in many scientific workloads on supercomputers. Some of these mitigation methods show remarkable success in the sense that faults can be compensated algorithmically by recovery procedures with only little extra cost in time or in silicon. On the other hand it also becomes clear that integrating such techniques in a computational infrastructure is still facing many obstacles. This includes the still poorly defined interface between user-level fault mitigation techniques and system level functionality, as, it is, e.g., necessary to reliably and quickly detect a device (core, memory, ...) failure on a large parallel machine.
Despite its breadth, the article is far from being comprehensive. The selection of topics is a subjective overview of current research in the field of resilience for extreme scale computing and it delivers an outlook into possible and promising future research topics and solutions.