Olivier Dalle

ON REPRODUCIBILITY AND TRACEABILITY OF SIMULATIONS

Reproducibility of experiments is the pillar of a rigorous scientific approach. However, simulation-based experiments often fail to meet this fundamental requirement. In this paper, we first revisit the definition of reproducibility in the context of simulation. Then, we give a comprehensive review of issues that make this highly desirable feature so difficult to obtain. Given that experimental (in-silico) science is only one of the many applications of simulation, our analysis also explores the needs and benefits of providing the simulation reproducibility property for other kinds of applications. Coming back to scientific applications, we give a few examples of solutions proposed for solving the above issues. Finally, going one step beyond reproducibility, we also discuss in our conclusion the notion of traceability and its potential use in order to improve the simulation methodology.

INTRODUCTION

A scientific method defines the "principles and procedures for the systematic pursuit of knowledge involving the recognition and formulation of a problem, the collection of data through observation and experiment, and the formulation and testing of hypotheses" (Merriam-Webster Dictionary 2012.) Computer simulations have long proven to be useful for working on the latter activities, to collect data from a computerized representation of the system under study (called the simulation model), and to test hypotheses by varying the parameters and inputs of the model.

Advanced and commercial simulation software are generally designed for making the life of the experimenter as easy as possible. Unfortunately, by simplifying the process of running an experiment, this software often hides important parameters and experimental design questions. For example, as reported by [START_REF] Perrone | Simulation of large scale networks II: modeling and simulation best practices for wireless ad hoc networks[END_REF], these hidden details end up missing in the published description of the experiment, which prevents a third party from reproducing the experiment with another simulation software, or even with a different version of the same simulation software.

The pressure to publish or perish, combined with the short time allotted for reviewing, and in some cases the unavailability of the source code, almost established as a standard the acceptance of a simulation result without any attempt at reproducing the actual simulation experiments. [START_REF] Ince | The case for open computer programs[END_REF] claim, in their editorial article for the journal Nature, that " [they] have reached the point that, with some exceptions, anything less than release of actual source code is an indefensible approach for any scientific results that depend on computation, because not releasing such code raises needless, and needlessly confusing, roadblocks to reproducibility." Going further than raising the issue, some journals such as Science already mandate that authors supply the source code along with their paper submission, while others, such as Biostatistics, have appointed an "Associate Editor for Reproducibility" [START_REF] Niemeyer | Nature Editorial: If you want reproducible science, the software needs to be open source[END_REF].

In this paper, we discuss reproducibility in the case of simulation-based experiments. We first give the multiple definitions of reproducibility in this particular context, including its degenerate forms, in section 2. Then, we identify and discuss a number of issues that need to be solved in order to achieve this property, in section 3. We discuss the various usages of reproducibility depending on the type of the simulation application, in section 4, and then we propose some solutions in section 5. Finally, in section 6, we present some related works, before we conclude in section 7.

DEFINITIONS

As a generic science concept, reproducibility designates the ability to confirm the results of a previous experiment by means of another similar experiment. In physics, for example, the criteria for reproducibility and repeatability are standardized with respect to results of measurements (NIST 1994). However, the previous definitions do not apply well to the particular context of in-silico experiments. For example, the conditions for the repeatability of results measurements, as defined by U.S. National Institute of Standards and Technology (NIST, op. cit.), requires the same observer, the same instrument used in the same conditions and the same location. These conditions and definitions clearly need to be adapted to the particular case of simulations.

Applied to in-silico experiments (ie. computer simulations), we consider reproducibility as the ability to produce similar data by running a similar simulation. The similarity of data does not mean that data must be strictly identical, but that they must exhibit the same properties according to some metric of interest, and usually in accordance with the goal of the study (some of the conditions defined by NIST could certainly apply here, such as the 95% similarity threshold). For example, data samples collected with two similar experiments may have different numerical values but have the same statistical properties or belong to the same curve with different sampling coordinates.

Before proceeding further with the definition(s) of reproducibility of simulations, we need first to introduce a few elements of terminology to set up the landscape of this paper. These definitions are intended to clarify the content of this paper rather than the body of knowledge of Modeling & Simulation at large or the many variants of the general terms that could be found in literature:

• A simulation is short for simulation run, which is the process of executing a computer program that runs a computer simulation;

• A computer simulation is composed of the following elements:

a model defines the structure and behavior of the system that is the subject to the study; the previous system can take many forms including computer and networking systems, biological systems, transportation networks, factory plants, crowds, battlefields, etc.

a scenario defines a possible history of the system defined by the model; the scenario can be produced either using a dedicated means (ie. a scripting language or some formatted inputs provided to the computer simulation program) or using the same elements as those used for definition of the model (a scenario is a model of the exogenous stimuli applied to the system);

a simulator refers to a reusable simulation engine and its Application Programming Interface (API); the engine is reusable for the simulation of many models and scenarios; an instrumentation defines which variable(s) of interest to observe during the simulation execution and what computation(s) to apply at run-time to the data samples produced by these variables (including their logging to result files).

When referring to computer simulations, we suggest that reproducibility should be considered according to the four levels of requirement presented hereafter: L1. deterministic, identical computation: this level of compliance is also known as the repeatability property; it designates the ability to re-execute the exact same simulation history as previously in terms of computation, including the inherently non-deterministic computations. In the case of a discrete-event simulation, for example, repeatability means that the simulations produce the same series of events and process them in the same order. In the case of distributed simulations, in particular, repeatability implies that parallel activities (eg. event occurring at the same simulation time) must always be processed in the same order, which requires some additional synchronization [START_REF] Fujimoto | Parallel and Distributed Simulation Systems[END_REF]. Although not strictly necessary, repeatability usually comes with the idea of reusing the same simulation source code in order to produce the same computation; this level implicitly requires an identical instrumentation;

L2. non-deterministic, identical computation: this level of compliance is a relaxed form of the previous one in which the non-deterministic situations are not resolved. This lighter form may either result from the lack of additional synchronization, or in some extreme cases, from the inability to fully control the experimental conditions, such as with real-time distributed simulation (see McLean and Fujimoto 2000, for example). Note that it is not considered good practice to produce significantly different results due to non-determinism in the model. Indeed, in practice, if two events occur at the same time and their processing order is not specified by the model, the order of their execution should not influence the simulation result, given that the result should only be correlated to the model, never to the simulator. For example, if two customers A and B arrive and get served at the same time at two different counters, then simulating customer A before customer B or the opposite should not change the total amount of transactions processed by that bank at the end of the day. Assuming that different instrumentation produces different outputs, which have resulted from different computations, achieving reproducibility at this level requires instrumentation to be identical;

L3. identical scenario and instrumentation: this level consists in reproducing the simulation experiment based on the detailed specifications of its scenario and model rather than re-using the same program source code. By neglecting some implementation details (eg. the sequence of random numbers produced by a random number generator), reproducibility at this level may lead to the production of different data, even though the instrumentation is assumed to be identical. Considering that the scenario should have been designed to observe a phenomenon, then the various possible executions of that scenario should reproduce the same phenomenon, and therefore the data, although different, should have identical meaning from a statistical point of view;

L4. similar scenario and instrumentation: this level is the most general and relaxed form of reproducibility. In general science, experiments may be considered similar even if they do not apply to the exact same subjects (qualitatively or quantitatively), especially when it is impossible to reproduce the exact same experimental conditions. Although computer simulation offers a fine control on the experimental conditions, it is still relevant to consider similar cases.

Compared to other forms of experiments, Modeling & Simulation has two distinctive features with respect to reproducibility:

• ability to replay: the ability for computers to reproduce in a deterministic way the exact same computation (the repeatability property discussed above at level L1), is a distinctive feature of simulation. For example, simulation allows for a replay of the same history of the system under study, but with a different instrumentation, in order to make different observations. Not only can the risk of interference (the Heisenberg principle) be avoided with computer-based experiments, but with the help of proper techniques, the instrumentation can even be totally separated from the model and scenario [START_REF] Dalle | An instrumentation framework for component-based simulations based on the Separation of Concerns paradigm[END_REF]. This leads to a different form of reproducibility, in which the experiment is reproduced, but with different outputs.

• fully virtual nature: when building in-silico experiments using simulation, the whole experiment, including the system under study, is artificially created. While this virtual nature offers a great flexibility, it also introduces new risks, that makes it a distinctive feature of simulation. Indeed, in other forms of experiments, the physical reality of the system under study means that at least that part of the system can be trusted and considered error-free, while the experiment process and the observation should be considered error-prone.

In simulation, not only the same sources of error exist, but the system itself can be modeled incorrectly and be a source of error. This situation leads to another form of reproducibility, in which the experiment is reproduced, but with a different system, that intrinsically produces different outputs.

ISSUES FOR REPRODUCIBILITY

In this section we discuss issues for reproducibility that are assumed to have some solution. Indeed situations may exist in which reproducibility is not possible due to some intrinsic property of the simulation. For example, as will be discussed later in section 4, the fact that a Virtual Reality simulation is driven by user interactions in real time makes it difficult to reproduce and is therefore not considered as an issue worth to address.

If we except the case of repeatability, issues for reproducibility fall into two main categories: human factors and technical issues.

Human factors

A number of human factors may impede reproducibility. Although negligence is certainly a primary concern, we won't elaborate further on this particular aspect and will assume that human-related operations are handled by operators always willing to contribute to the best extent of their current knowledge. Other human factors left to consider include the following:

• "genuine" ignorance: Being unaware of some hidden parameters, such as which random number generator is being used by default by the simulation engine or library;

• insufficiently detailed publications: Despite the fact that this issue is widely acknowledged as being major, it may be considered a result of the publication process as a whole, making responsibility very difficult to attribute. Some details may be skipped by authors either inadvertently or on purpose, because the publication space is limited. Additionally, there is very little incentive to provide reproducible content because this aspect is almost never checked by reviewers, either because the reviewing deadlines are too tight, or because publishers do not give that instruction;

• business related limitations: As mentioned in introduction, availability of the source code seems necessary to ensure reproducibility. However, providing an Open Source simulation software is not always an option. Firstly, some commercial simulators are commercial products with a business model that is not designed for Open Source. Furthermore, as will be explained in greater details in the next section, simulation has many applications that are not aimed toward producing science and, therefore, do not necessarily require a reproducibility based on the source code availability. Secondly, when using an open source simulator, some additional elements, such as the models and scenario used to produce the results, may not be fully disclosed at source level, for intellectual property reasons;

• manipulation errors: These errors result from the manual handling of some of the tasks in the simulation study work-flow. Trying to build a detailed typology of such errors is pointless (and certainly endless). However, assuming that their probability of occurrence is in proportion with the amount of human intervention required in simulations, human involvement should be minimized as much as possible and replaced by automation, as will be later discussed in section 5.

Technical issues

A variety of technical issues may also arise, including the following ones:

• software bugs: Because of software bugs in simulation engines and libraries, perfectly valid experiments and models can produce erroneous results that are unreproducible in newer versions of the simulator once the bugs are fixed. Being able to reproduce these erroneous results with old versions, however, is still useful, because that can help to explain other discrepancies or track erroneous results in time. Unfortunately, the various versions of a simulation software, especially the ones with bugs, are not always available;

• software availability: Multiple reasons may explain software unavailability, including losing track of old versions, commercial reasons (eg. discontinued or unaffordable products) or confidentiality and competitive reasons (eg. providing a good simulator of a system to a competitor may not be desirable).

• floating points numbers: Although very useful, floating point numbers (FP) are reputed to be an endless source of complications, as explained by [START_REF] Goldberg | What Every Computer Scientist Should Know About Floating-Point Arithmetic[END_REF]. For example, in the popular networking simulator ns-2, the time of occurrence of events is encoded using FP and it is easy to observe discrepancies when scheduling repetitive tasks during a simulation: Due to FPs errors, repeatedly adding a little delay does not produce the same date as multiplying an incremental delay by the number of iterations. Indeed, although semantically identical, ns-2 considers the dates 0.99999999999 and 1.0000000001 as different. This kind of error is likely to introduce a bias in the results given that all events occurring at time 0.9 -9 end up having more priority than those occurring at time 1.0 -01. Such a bias may be unreproducible with another simulator in which time is encoded using integers;

• computer and Operating Systems evolutions: Ensuring the long term reproducibility of a simulation may be a challenging issue, given the continuous evolution of computer systems and architectures. Depending on the language used for writing the simulation, a significant dependency may exist with the operating system or additional libraries.

USAGE EXAMPLES & APPLICATIONS

In this section, we discuss the various usages of reproducibility depending on the type of the simulation application.

In section 4.1, we first use the example of a network performance evaluation as a typical case of scientific study.

Then through the following sections, 4.2 to 4.5, we discuss Road Traffic simulations, Training, Gaming and Virtual Reality as examples of non scientific studies. The summary of these various cases is then presented in section 4.6. These examples show that reproducibility at various levels may be considered and even have value for simulations beyond the limited scope of scientific studies.

Network performance evaluation

This type of simulation is normally used as part of a scientific study and usually consists in assessing the performance of a protocol or network architecture subject to a typical traffic workload. The simulation consists in reproducing the protocol(s) operations while the simulated traffic is flowing through the network. Multiple simulations are necessary to estimate the metrics of interest, such as the level of Quality of Service achieved end-to-end, or the average size of the packet queues in the network routers. Parallel simulations may be employed for large or complex networks, which explain the need for reproducibility at level L1. A special form of network simulation, called emulation, consists in coupling real hardware or software with the simulation. In this case, the non-determinism due to the real-time operation of the network may limit the reproducibility to level L2. Compared to other types of simulation studies, a network simulation study usually requires to develop new models.

Reproducibility at Level L3 is sometimes used as a technique of verification and validation: The same simulation is run using different simulators to assess the correctness of the new model implementation. Since network simulation is mostly used for scientific studies, reproducibility at level L4 is relevant.

Road traffic

This type of simulation is typically used in an operational context: For example, it may be used by a traffic management center, to compute the dynamics of traffic within a road network over a period of time. The road network is modeled in terms of lanes and geographical elements such as crossings, traffic lights, traffic circles, etc. The traffic is usually defined by an Origin-Destination matrix that describes the starting and ending points of each vehicle on the road network. A possible goal of this type of simulation is to identify congestion points either when the traffic changes (eg. following an event) or when the network changes (eg. a street or lane is closed). Parallel simulation may be employed, which explains the possible need for reproducibility at level L1. Reproducibility at level L2 does not seem to be necessary, even though road traffic simulators may be fed with real-time traffic information, because in this case the simulation still remains deterministic. Reproducibility at level L3 is not necessary in operational context because the simulator is often used to make decision in real time and because writing new models for this type of simulation is not common, once the initial development stage is complete. Reproducibility at level L4 is not necessary because scientific studies are not the main use for this type of simulation.

Training

This type of simulation is designed to put users in operational situations according to some predefined scenario; they are used to train operators in particular tasks or resolve situations using a virtual copy of their working environment. For example, such simulator can be used to train nuclear plant operators, train or subway conductors, or to train plane and helicopter pilots to fly with no visibility using their instruments. Parallel simulation is not expected on a large scale and therefore issues related to reproducibility at level L1 are irrelevant. Repeatability at level L2 is a highly desirable feature, since repeating identical actions is an actual process of learning. Repeatability at level L3 is of marginal use, because it is not common practice to offer trainees several simulators for the same training purpose, however, a given simulator may be subject to the release of new versions, which then would require reproducibility at level L3. This type of simulator produces interactivity rather than outputs, and therefore reproducibility at level L4 is not applicable.

Gaming

Although very similar to the previous one, this type of interactive simulation has different requirements. Reproducibility at level L1 is irrelevant for the same reasons as for Training; reproducibility at level L2 is possible as part of the game design (many games follow strictly a pre-established scenario), but is not mandatory; Reproducibility at level 3 is not relevant, because most games have only one (commercial) implementation; and reproducibility at level L4 is not applicable.

Virtual Reality

This type of simulation is another form of interactive simulation in which a view on the virtually simulated world is computed in real-time, depending on the user's physical actions. In many such simulations, the user is a passive observer and the virtual world is either static or has its own dynamics, independent of the user's actions. Since the view generated in real-time by each simulation is driven by the user's action in real-time, each simulation is meant to be unique and therefore reproducibility at levels L1 or L2 is irrelevant. Reproducibility at level L3 is relevant only with more advanced forms of VR simulations in which the user becomes active and fully interacts with the VR world. These advanced forms may then be considered as either Gaming or Training depending on their intended use. Since reproducibility is impossible at level L2, it must be provided at level L3 if required. For example, a VR-based flight-combat simulator obviously requires the VR world to react to users' actions; it may be used for Training, and therefore requires the simulation to follow a pre-existing scenario. Since the simulation cannot be reproduced at level L2 because of the user's real-time interactions, however, only the simulation scenario can be reproduced. Reproducibility at level L4 is not applicable to VR simulations.

Summary

SOLUTIONS FOR REPRODUCIBILITY

In this section we discuss possible solutions for reproducibility issues in scientific studies based on simulations.

In section 2, we distinguished among four levels of reproducibility, which in turn suggested specific solutions at corresponding levels, which we presented in sections 5.2 to 5.4. However, some issues are level-independent, which we will first discuss in section 5.1.

Level independent issues

Human errors are an example of level-independent issues. As already suggested in introduction, a pragmatic way to solve these issues is to reduce human intervention by relying on automation as much as possible. The Simulation Automation Framework for Experiments project [START_REF] Perrone | SAFE: Simulation Automation Framework for Experiments[END_REF], is an example of the ongoing development of such a solution for the ns-3 network simulator.

Solutions for Repeatability (levels L1 & L2)

In the introduction we mentioned several works on protocols to solve the non-determinism issues of parallel execution [START_REF] Fujimoto | Parallel and Distributed Simulation Systems[END_REF] or real-time interactions [START_REF] Mclean | Repeatability in Real-Time Distributed Simulation Executions[END_REF]. If we except these determinismrelated issues, repeatability is mainly a question of ensuring the long-term execution of a program. This issue is challenging given that computer architectures and operating systems and their accompanying libraries are both continuously evolving. The following two pragmatic solutions can be proposed for this evolution issue:

• use a platform-independent solution: platform independence can be provided at various levels, but the most efficient is certainly at the language level, assuming the language is stable due to its large user community. This level of independence can be achieved by adopting a simulation formalism, such as Petri Nets or DEVS [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF], especially in association with one their well maintained software platform (eg. CPNTools, [START_REF] Jensen | Coloured Petri Nets -Modeling and Validation of Concurrent Systems[END_REF] or in accordance with a standard (eg. DEVS standardization, [START_REF] Wainer | Standardizing DEVS model representation[END_REF]).

• ensure the availability of the computing solution: virtual appliances provide a timely and straight-forward implementation of this solution. A virtual appliance is a virtual machine fully configured to run a particular application or service, including the application software and all of its dependencies. Since the virtual appliance is a self-contained snapshot of a functional system, it only requires the continuous support of virtual machine. This solution can also solve the issue related to the availability of the source code: Since with a virtual appliance, the execution is fully handled by the virtual machine at the binary level, the source code is not necessary.

Solutions for level L3

At this level, the challenge is to reproduce the model and scenario of a previous simulation study. For this purpose, an obvious requirement is the ability to retrieve the corresponding models and scenario specifications, which become more and more challenging to retrieve as time goes on. Long-term archives of the model, the scenario, and to a larger extent of any simulation experiment is a critical issue when using simulation for scientific purposes. A solution to this problem could be to place the responsibility in the hands of the publishers and editors of scientific publications. As already mentioned in introduction, a few publications already enforce such a policy. However, this policy comes at a cost and the real difficulty is therefore to offer enough incentives to both overcome this overhead and deeply change the practices. A summary of the changes required to enforce such a policy is given hereafter:

• The sources of the code required to reproduce the simulations should be provided with each paper submission that contains a scientific result obtained by means of simulations;

• The previous requirement should be enforced by editors and publishers;

• This additional material should be part of the publication and, as such, archived and redistributed by the publishers along with the printed material of the publication;

• The review process should include a review of the provided material and at least one test of reproducibility.

Enforcing such a policy comes with a number of positive and negative side-effects listed hereafter:

-Publishers would need to adapt their archival and diffusion means to support the provided source code;

-Editors would have to account for this new form of submitted material in their decisions; for example, it might be necessary to ensure that the reviewers are sufficiently familiar with the type of source code provided and able to check the reproducibility;

-The peer-review process becomes more complex and takes more time;

+ As part of the publication, the provided source code can be cited;

+ The general availability of the source code makes the reproduction and extension of previous published works much easier, which is a clear benefit to the whole scientific community;

+/-This policy enforces the use of Open Source software, which may not be acceptable for commercial or strategic reasons, but would certainly be largely supported within the scientific community;

Note that the source code required to reproduce the simulations may extend well beyond the code of the models and experiments. Depending on the dependencies existing among the results of the simulations and other elements, it may be necessary to include the code of the simulator engine and libraries, as well as those of the operating system. However, given the additional amount of space required for the storage of these elements, we can assume that publishers would encourage reusing previously published ones, or would put restrictions on attempts at publishing new ones (eg. request a justification from the authors explaining why they need a new configuration instead of the one already published).

Solutions for level L4

At this level we assume that we want to reproduce the results of a previous experiment, but not necessarily using the source code. The issues about using Floating Points numbers and insufficiently detailed publications are typically observed at this level. Obviously, avoiding Floating Point as much as possible is a way of solving the former, while the latter leads again to the publication process, as discussed in the previous section. However, the requirement on the publications becomes lighter at this level, compared to the level L3: if we assume that providing the source code is no longer necessary, then any sufficiently detailed specification of the experiment may be accepted. Unfortunately, this lighter requirement would make the task of reviewers more difficult than it was in the previous level. Not asking to do this, however, would leave the issue of insufficiently detailed publications unsolved.

RELATED WORKS

The idea of reproducible research applied to computations seems to be first presented in [START_REF] Schwab | Making Scientific Computations Reproducible[END_REF], which opened a new research direction (see [START_REF] Fomel | Introduction to: Special Issue on Reproducible Research[END_REF]. [START_REF] Rahmandad | Reporting Guidelines for Simulation-based Research in Social Sciences[END_REF] discuss the concepts of "Minimum Model Reporting Requirement" and "Preferred Model Reporting Requirement" as a means of ensuring reproducibility in simulation-based research in Social Science, at level L3, even though the authors do not explicitly distinguish multiple levels. An interested reader will also find a number of related works in various domains in the previous paper, including biochemistry, biology, computational research, and mathematics.

However, many published works based on simulation still fail to meet the minimum conditions to ensure reproducibility. For example, [START_REF] Kurkowski | MANET Simulation Studies: The Incredibles[END_REF] surveyed the results of MANET simulation studies published between 2000 and 2005 in the ACM MobiHoc Symposium: 75% of these papers used simulation but less than 15% only of the simulations were repeatable and only 7% addressed such important issues as initialization bias.

This failure to ensure reproducibility is one of the reasons that concerns about credibility are raised in simulation. About this general issue of credibility: [START_REF] Pawlikowski | On Credibility of Simulation Studies Of Telecommunication Networks[END_REF] surveyed over 2200 publications on telecommunications networks in proceedings of the IEEE INFOCOM and such journals as the IEEE Transactions on Communications, the IEEE/ACM Transactions on Networking, and the Performance Evaluation Journal and their conclusion was that "[...] the majority of recently published results of simulation studies do not satisfy the basic criteria of credibility." The two main criteria they consider in order to assess credibility were the use of appropriate pseudo-random numbers generators and the appropriate analysis of simulation outputs.

This credibility issue is surprising given that the many steps followed, and decisions taken in simulation studies have been quite extensively covered in the literature and, in particular, in the simulation reference books [START_REF] Banks | Discrete-Event System Simulation[END_REF], Law 2007, etc.) These steps and decisions are usually described by means of a typical work-flow that starts with the definition of the goal of the study and proceeds, with possibly several iterations, through the major steps of modeling, validation, verification, and experiment planning, and concludes with an analysis of the results produced by the simulation(s). Although some of these work-flows depicted in literature are quite detailed (eg. [START_REF] Wainer | Discrete-Event Modeling and Simulation: a Practitioners approach[END_REF], they all choose the same activity perimeter that starts with the goal definition and ends with data analysis. This already large, but still limited perimeter, explains why the issues of reproducibility are seldom covered: they occur at a higher level.

Indeed, reproducibility is part of a larger scientific process that extends well beyond the boundaries of Modeling & Simulation, both in terms of scope and duration. In terms of scope, for example, a physical experiment can be reproduced using simulations in order to confirm the physical results. In terms of duration, for example, a simulation could be reproduced a long time after a study was completed (and possibly published), either as part of another study, or to confirm published results.

CONCLUSIONS

In this paper, we discussed the issues and current status of simulation reproducibility. We first made explicit the multiple levels of reproducibility that can be achieved using computer simulations. We explored their potential use in various forms of simulation, and found a few in which they might be relevant. Then, we reviewed a number of issues that need to be solved in order to achieve reproducibility in the particular case of scientific applications and we discussed some solutions to solve them at the various levels we identified.

Our analysis reached similar conclusions to [START_REF] Ince | The case for open computer programs[END_REF], which advocate the use of Open Source software for scientific publications. We also discussed the publication process itself and how it should evolve in order to ensure the reproducibility of simulations based on Open Source software. Although these recommended change would have significant impact on the current practices for all the actors involved, we strongly believe that they would result in greater benefits for the whole scientific community.

Last but not least, solving the reproducibility issues also opens new perspectives, such as allowing for the traceability of simulations results.

Traceability is the ability to trace the dependencies of results and publications to the software elements that were used to produce them. In conjunction with reproducibility, this ability is of particular interest when a bug is discovered, because it can help identify (trace) which results and publications are impacted by the bug. Since the simulation software and models are continuously subject to improvements and new versions, ensuring traceability rapidly becomes a problem of version and dependency management. In these cases, the most appropriate level of reproducibility is level L2, so that the code of the impacted simulation can be re-executed after the bug is fixed. However, in order to cope with more complex forms of errors, it is certainly desirable to allow the reexecution of whole experiments, which may involve a number of complex tasks. Solutions have emerged in the last few years, such as the Apache Foundation's Maven project for management of dependencies and versions or the MyExperiment and Bonita projects for the management of work-flows [START_REF] Dalle | Should Simulation Products Use Software Engineering Techniques or Should They Reuse Products of Software Engineering? Part 1[END_REF].

Therefore, used conjointly, these two properties of reproducibility and traceability would allow for the proper detection and correction of errors in published results obtained using simulations.

Table 1 :

 1 Table1gives a summary of the usage (and benefits) for reproducibility in terms both of the level of reproducibility and the type of the simulation application. It shows that despite the main use by and stringent requirements of the scientific field, there are still some potential applications of reproducibility in other types of simulations. Nevertheless, in the remaining sections of this paper, we will only concentrate on reproducibility in the context of scientific studies. Usage for each level of reproducibility depending on the application type

	Simulation		Level of Usage
	application	Level L1	Level L2	Level L3		Level L4
	Network perfor-	yes, if parallel	yes, if RT emul.	yes, verification	yes, 3d party vali-
	mance evaluation	exec.				dation
	Road traffic	possibly, if paral-	no	no		no
		lel exec.
	Training	no, irrelevant	yes	marginal,	with	N/A
				new versions
	Gaming	no, irrelevant	possibly, by de-	no, irrelevant		N/A
			sign
	Virtual reality	no, irrelevant	no, irrelevant	yes, if active user	no

ACKNOWLEDGMENTS

This work has been co-funded by the DISSIMINET Associated Team initiative (INRIA Sophia Antipolis and Carleton University) and the French Agence Nationale de la Recherche (ANR) SONGS Project. I would like to address my special acknowledgements to Joe Peters and the Faculty of Applied Sciences at the Simon Fraser University for hosting me while writing this paper and for their funding support through an Ebco Eppich Fellowship grant. Last but not least, I would like to kindly acknowledge the three anonymous reviewers of this paper for their insightful comments and suggestions, and my wife for her careful and professional English corrections.


P Gosselet

C Rey

J Pebrel Lmt

Total and selective reuse of Krylov subspaces for the resolution of sequences of nonlinear structural problems

Keywords: Krylov solvers, multiresolution, model reduction

published or not. The documents may come

Total and selective reuse of Krylov subspaces for the resolution of sequences of nonlinear structural problems

Pierre Gosselet, Christian Rey, Julien

Introduction

Accelerating the convergence of Krylov iterative solvers [START_REF] Saad | Iterative methods for sparse linear systems[END_REF] is an old issue which has returned to the spotlight because of the increasing number of applications for which these are preferred to direct solvers today. Traditional approaches aim at improving the condition number by using frameworks in which efficient preconditioners exist (e.g. domain decomposition methods [START_REF] Klawonn | FETI and Neumann-Neumann iterative substructuring methods: connections and new results[END_REF][START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF]), or for which good initialization vectors [START_REF] Gosselet | On the initial estimate of interface forces in FETI methods[END_REF], relevant augmentation subspaces [START_REF]Conjugate gradient method with preconditioning by projector[END_REF][START_REF] Saad | Analysis of augmented Krylov subspace methods[END_REF][START_REF] Chapman | Deflated and augmented Krylov subspace techniques[END_REF] or suitable block strategies (see [START_REF] Aliaga | A Lanczostype method for multiple starting vectors[END_REF] for a very general block-Lanczos algorithm) are available. For instance:

• For 3D elasticity problems, domain decomposition methods come with "physical" augmentation associated with the global equilibrium of floating substructures (rigid body motions), which makes the methods scalable [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF][START_REF] Mandel | Balancing domain decomposition for problems with large jumps in coefficients[END_REF]; for plate and shell problems, additional augmentation through "corner modes" [START_REF] Farhat | The two-level FETI method for static and dynamic plate problems -part I: An optimal iterative solver for biharmonic systems[END_REF][START_REF] Farhat | The two-level FETI methodpart II: Extension to shell problems. parallel implementation and performance results[END_REF][START_REF] Le Tallec | A Neumann-Neumann domain decomposition algorithm for solving plate and shell problems[END_REF] is required.

• For structures with repeated patterns, block strategies are possible [START_REF] Gosselet | A domain decomposition strategy to efficiently solve structures containing repeated patterns[END_REF].

1

• For restarted algorithms, one can use deflation or augmentation [START_REF] Morgan | A restarted GMRes method augmented with eigenvectors[END_REF][START_REF] Erhel | Restarted GMRes preconditioned by deflation[END_REF], or block techniques [START_REF] Baker | On improving linear solver performance: A block variant of GMRes[END_REF].

• For problems with multiple right-hand sides, deflation [START_REF] Farhat | Extending substructure based iterative solvers to multiple load and repeated analyses[END_REF][START_REF] Saad | h. A deflated version of the conjugate gradient algorithm[END_REF][START_REF] Erhel | An augmented conjugate gradient method for solving consecutive symmetric positive definite linear systems[END_REF] is a rather classical approach.

The problem with these techniques is that they require some a priori information which is seldom available, except in specific cases. Many recent works present theoretical and practical comparisons of the numerous algorithms which have been developed in connection with these ideas [START_REF] Tang | Theoretical and numerical comparison of various projection methods derived from deflation, domain decomposition and multigrid methods[END_REF][START_REF] Simoncini | Recent computational developments in Krylov subspace methods for linear systems[END_REF].

Multiresolution approaches form a general framework in which numerical information is available to accelerate the convergence of Krylov solvers. Multiresolution refers to situations in which the solution of a mechanical problem cannot be achieved through the resolution of a single linear system. For example, calculating the solution of a nonlinear or time-dependent problem or exploring a design of experiment during an optimization procedure requires the resolution of sequences of linear systems. Multiresolution is more general than using multiple right-hand sides because the matrices themselves are likely to change from one system to another (multiple right-hand and left-hand sides). Thus, the problem consists in solving a k-indexed family of large, sparse, linear n × n systems of the form:

A (k) x (k) = b (k) (1)
Although different, the systems are assumed to be similar to one another. This similarity can be defined in several ways: in terms of rank, by the fact that rank(A (k) -A (k-1) ) rank(A (k) ); or in a spectral sense by the fact that the eigenspaces remain stable from one system to another; or in terms of the Krylov subspaces generated [START_REF] Carpraux | Stability of the Krylov bases and subspaces[END_REF]. The first case can be dealt with easily, even with direct solvers, by using the Sherman-Morrison formula; the second case requires augmentation strategies in order to eliminate the most penalizing part of the spectrum and improve the active condition number 1 [START_REF] Rey | A Rayleigh-Ritz preconditioner for the iterative solution to large scale nonlinear problems[END_REF][START_REF] Gosselet | On a selective reuse of Krylov subspaces in Newton-Krylov approaches for nonlinear elasticity[END_REF][START_REF] Parks | Recycling Krylov subspaces for sequences of linear systems[END_REF][START_REF] Wang | Large-scale topology optimization using preconditioner Krylov subspace methods with recycling[END_REF]; and the last case calls for preconditioning techniques [START_REF] Rey | An acceleration technique for the solution of non-linear elasticity problems by domain decomposition[END_REF][START_REF] Risler | Iterative accelerating algorithms with Krylov subspaces for the solution to large-scale nonlinear problems[END_REF].

While most of the studies of Krylov methods for multiresolution (often referred to as the recycling of Krylov subspaces) are set in the framework of GMRes/MinRes [START_REF] Parks | Recycling Krylov subspaces for sequences of linear systems[END_REF][START_REF] Wang | Large-scale topology optimization using preconditioner Krylov subspace methods with recycling[END_REF], we chose to work on the specific case of the resolution of symmetric, positive definite systems using conjugate gradients (CGs), in which the convergence is under control and related to easily calculated spectral properties [START_REF] Van Der Sluis | The rate of convergence of conjugate gradients[END_REF]. In earlier works, the authors developed efficient preconditioners based on previous Krylov subspaces [START_REF] Rey | An acceleration technique for the solution of non-linear elasticity problems by domain decomposition[END_REF][START_REF] Risler | Iterative accelerating algorithms with Krylov subspaces for the solution to large-scale nonlinear problems[END_REF] which took advantage of the conjugation properties of CGs, but did not extract the most interesting part of the information available in the Krylov subspaces, and they proposed augmentation techniques using Ritz vectors [START_REF] Rey | A Rayleigh-Ritz preconditioner for the iterative solution to large scale nonlinear problems[END_REF]. Typically, these works were aimed at nonlinear mechanical systems solved by Newton-Raphson linearization and FETI or BDD domain decomposition [START_REF] Gosselet | A domain decomposition method for quasi incompressible formulations with discontinuous pressure field[END_REF][START_REF] Léné | Some strategies to compute elastomeric lamified composite structures[END_REF].

The recycling of Krylov subspaces can also be analyzed from the model reduction point of view. Since Krylov solvers satisfy Petrov-Galerkin conditions, they share many common points with strategies based on Karhunen-Loeve expansion [START_REF] Meyer | Efficient model reduction in non-linear dynamics using the karhunen-loeve expansion and dual-weighted-residual methods[END_REF][START_REF] Ryckelynck | On the "a priori" model reduction: Overview and recent developments[END_REF][START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF]. These similarities are well-known [START_REF] Freund | Krylov-subspace methods for reduced-order modeling in circuit simulation[END_REF][START_REF] Freund | Model reduction methods based on Krylov subspaces[END_REF][START_REF] Heres | Combining Krylov subspace methods and identification-based methods for model order reduction[END_REF]. But our objective is not to develop reduced models of mechanical systems in order to perform fast but coarse analyzes; it is to define, improve and reuse reduced models in order to carry out calculations both rapidly and accurately.

In this paper, we undertake a more in-depth investigation of augmentation using a selection of post-processed Ritz vectors. In Section 2, we begin with a detailed presentation of the theoretical framework of the augmented preconditioned conjugate gradient method; then, in Section 3, we propose a first reuse algorithm in a multiresolution framework (TRKS); in Section 4, we improve this algorithm by proposing a procedure for selecting the "best" Ritz vectors (SRKS and "cluster"); finally, in Section 5, we propose an evaluation of the method in the case of nonlinear mechanics and parametric problems, using domain decomposition methods [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF] to define efficient preconditioners.

2 The augmented preconditioned conjugate gradient method

Algorithm and properties

Let us consider the linear problem

Ax = b, (2)
where A is an n × n symmetric positive definite matrix, and let us study the resolution of this system using the augmented preconditioned conjugate gradient algorithm. With M being the n × n symmetric positive definite matrix of the preconditioner, we introduce the following notations: i = 0 . . . m the iteration number

x i the i th approximation r i = b -Ax i = A(x -x i ) the i th residual (3)
With no loss of generality, the presentation can be limited to the case of a zero initial guess x 00 = 0. (Otherwise, one can set b ← b -Ax 00 .) Let C be a subspace of R n of dimension n c , and let Matrix C = [c 1 , . . . , c nc ] be a basis of C. The search principle of the augmented left-preconditioned conjugate gradient is:

find x i ∈ K i (M -1 A, C, M -1 r 0 ) such that r i ⊥ K i (M -1 A, C, M -1 r 0 ) (4)
where K i (M -1 A, C, M -1 r 0 ) is the augmented Krylov subspace associated with preconditioned operator M -1 A and augmentation subspace C:

K i (M -1 A, C, M -1 r 0 ) = span M -1 r 0 , . . . , (M -1 A) (i-1) M -1 r 0 ⊕ C (5)
A classical implementation relies on the definition of a convenient initialization and projector pair (x 0 ,P ):

x = x 0 + P y C T r 0 = 0 ⇐ x 0 = C(C T AC) -1 C T b C T AP = 0 ⇐ P = I -C(C T AC) -1 C T A (6)
One should note that since AP = P T AP = P T A augmentation preserves symmetry. One should also note that P C = 0. The system to be solved is:

AP y = (P T AP )y = r 0 = P T b (7)
The C-augmented, M -preconditioned conjugate gradient technique (APCG) implemented by projection is presented in Algorithm 1. (For the sake of simplicity, the methods will be described assuming exact arithmetic, even though they are compatible with more realistic full reorthogonalization [START_REF] Lingen | Efficient Gram-Schmidt orthonormalisation on parallel computers[END_REF].)

Algorithm 1: APCG(A, M , C, b) Calculate AC, (C T AC) -1 ; (P = I -C C T AC -1 C T A); x 0 = C(C T AC) -1 C T b; r 0 = b -Ax 0 = P T b; z 0 = P M -1 r 0 , w 0 = z 0 ; for j = 1, . . . , m do α j-1 = (r j-1 , w j-1 )/(Aw j-1 , w j-1 )
x j = x j-1 + α j-1 w j-1 r j = r j-1 -α j-1 Aw j-1 z j = P M -1 r j+1 w j = z j -β j w j-1 β j = (Aw j-1 , z j )/(w j-1 , Aw j-1 ) end

The following basic relations hold:

(r i , z j ) = 0, i = j (w i , Aw j ) = 0, i = j (8)
With W i = [w 0 , . . . , w i-1 ] and Z i = [z 0 , . . . , z i-1 ] being two bases of K i (P M -1 A, z 0 ), the projector enables the spaces to be divided orthogonally:

K i (M -1 A, C, M -1 r 0 ) = K i (P M -1 A, z 0 ) ⊥ A ⊕ C (9)
Of course, in the absence of optional constraints (C = 0, P = I), APCG reduces to standard preconditioned conjugate gradients PCG(A,M , b); if, in addition, M -1 = I, it becomes a standard conjugate gradient algorithm CG(A, b).

Let us recall a first result which was proven in [START_REF]Conjugate gradient method with preconditioning by projector[END_REF] for the case of nonpreconditioned augmented conjugate gradients. Proposition 1. Let V = Range(P )

• APCG(A, I, C, b) is equivalent to CG(P T AP |AV , P T b) in the sense that both generate the same residuals. x i , the i th APCG approximation, is connected to y i , the i th CG approximation, by x i = x 0 + P y i .

• APCG(A, I, C, b) does not break down; it converges, and its asymptotic convergence rate is governed by the condition number κ P T AP |AV κ (A).

Consequently, augmentation strategies never decrease the asymptotic convergence rate. The following corollary is straightforward:

Corollary 1. Let D = [d 1 , . . . , d m d ] be a set of m d linearly independent vec- tors such that E = [C, D] is a full column rank matrix. Let P E = I - E(E T AE) -1 E T A
, and let V E be the range of P E .

• APCG(A, I, E, b) is equivalent to APCG(P T AP , I, D, P T b) in the sense that both generate the same residual. x E i , the i th approximation of APCG(A, I, E, b), is connected to x D i , the i th approximation of APCG(P T AP , I, D, P T b), by

x E i = x 0 + P x D i . • The asymptotic convergence rate is governed by κ P E T AP E |AV E κ P T AP |AV κ (A).
In conclusion, an increase in the size of the augmentation can only improve the asymptotic rate of convergence. (In the worst case, it leaves it unchanged.) Now let us focus on the effect of preconditioning. Since M is a symmetric positive definite matrix, it can be factorized in Cholesky form M = LL T (where L denotes a lower triangular matrix with positive diagonal coefficients). Let us introduce the notation:

 = L -1 AL -T ; b = L -1 b ; x = L T x Ĉ = L T C ; P = I -Ĉ( ĈT  Ĉ) -1 ĈT  (10)
Then, the following equivalence between preconditioned and non-preconditioned augmented conjugate gradients holds:

Proposition 2. APCG(A, M , C, b) is equivalent to APCG( Â, I, Ĉ, b) with r = L -1 r = ẑ = L T z, ŵ = L T w, α = α and β = β. Its asymptotic convergence rate is governed by κ P T Â P| Â V κ Â . Proof. Since P = L T P L -T , we obtain directly x0 = Ĉ( ĈT Â Ĉ) -1 ĈT b = L T x 0 , r0 = b -Â x0 = L -1 r 0 = ẑ0 = L T z 0 and ŵ0 = L T w 0 .
By induction, it follows that αj-1 = (r j-1 , ẑj-1 )/( Â ŵj-1 , ŵj-1 ) = α j-1 , rj = L -1 r j , βj = ( Â ŵj-1 , ẑj )/( Â ŵj-1 , ŵj-1 ) = β j and ŵj = L T w j . Proposition 1 provides the inequality concerning the asymptotic convergence rate.

Putting these propositions together, APCG(A, M , C, b) is equivalent to CG(L -1 P T AP L -T , L -1 P T b). All these results lead us to propose an efficient augmentation by analogy with an equivalent, simpler system solved by classical conjugate gradients.

Interpretation and choice of the augmentation

From a "constraint" point of view, the projection guarantees the C-orthogonality of the residual throughout the iterations (C T r j = 0). For example, in the FETI domain decomposition method, the residual is the displacement jump between the subdomains; in the case of shell and plate problems, matrix C is introduced to enforce the continuity of the displacement at the corner points [START_REF] Farhat | The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems[END_REF]. In the BDD domain decomposition method, matrix C is associated with the rigid body motions of floating substructures and, therefore, local Neumann problems in the preconditioner are always well-posed [START_REF] Mandel | Balancing domain decomposition[END_REF]; for shell and plate problems, the matrix is enriched by corner mode corrections [START_REF] Le Tallec | A Neumann-Neumann domain decomposition algorithm for solving plate and shell problems[END_REF]. In both cases, matrix (C T AC) -1 , called a coarse grid matrix, plays a crucial role in the scalability of these methods.

From a "spectral" point of view, augmentation can be used to decrease the active condition number ("active" referring to eigenelements solicited by the right-hand side) and, thus, improve the asymptotic convergence rate. This is called a deflation strategy [START_REF] Erhel | An augmented conjugate gradient method for solving consecutive symmetric positive definite linear systems[END_REF][START_REF] Chapman | Deflated and augmented Krylov subspace techniques[END_REF], which boils down to building matrix C by using (approximate) eigenvectors associated with the lowest eigenvalues. Obviously, when C consists of the n c eigenvectors associated with the lowest eigenvalues (λ 1 . . . λ nc . . . λ n ), the condition number decreases strictly: κ P T AP AV = λn λn c < λn λ1 . From a "model reduction" point of view, subspace C represents a "macro" (or coarse) space in which the macro part of the solution is calculated directly during the initialization while the "micro" part of the solution, when required, is obtained during the iterations.

Estimation of computation costs

With regard to the numerical cost of augmentation, the main operations for the construction of the projector are: (i) the block product AC (and assembly with neighbors for domain decomposition methods), (ii) the block dot-product (C T AC) (plus an all-to-all sum for domain decomposition methods), and (iii) the factorization of the fully-populated coarse matrix (C T AC). Then, the application of the projector consists simply of (i) one block dot-product ((AC) T x) (plus an all-to-all exchange), (ii) the resolution of the coarse problem, and (iii) the matrix-vector product (Cα).

Thus, provided that the number of columns of matrix C is small, the main cost is related to the calculation of AC. One must bear in mind that block operations (on "multivectors") are comparatively much faster than single vector operations, especially when the matrices are sparse (because data fetching is factorized). In a domain decomposition context, product AC corresponds to the resolution of Dirichlet or Neumann problems in substructures, which makes the simultaneous treatment of many columns very efficient (and minimizes the number of exchanges). One must also remember that a conjugate gradient iteration involves a preconditioning step which may be expensive. (The cost is comparable to that of an operator product in optimal domain decomposition methods.) Thus, the additional cost of augmentation relative to the cost of one iteration depends on many parameters (the size of the problem, the number of augmentation vectors, the number of subdomains, the preconditioner chosen...). Typically, in the examples presented in this paper, we found that, using an optimal preconditioner, the CPU cost of between 4 and 7 augmentation vectors (depending on the hardware configuration) cost no more than one CG iteration.

A question which is not addressed in this paper is the verification of the full-rank property of matrix C, which affects the quality of the factorization of matrix (C T AC). Strategies to correct a dependence among the columns of matrix C due to inexact arithmetic can be found in [START_REF] Aliaga | A Lanczostype method for multiple starting vectors[END_REF].

Total reuse of Krylov subspaces

In this section, we show how it is possible to define efficient augmentation strategies in a multiresolution context. Let us consider the sequence of linear systems:

A (k) x (k) = b (k) , k = 1, . . . p (11)
where A (k) is an n × n symmetric positive definite matrix and b (k) is the righthand side. Each linear system is solved using an augmented preconditioned conjugate gradient algorithm APCG(

A (k) , M (k) , C (k) , b (k)
). Let m (k) be the number of iterations which is necessary to reach convergence, and let

W (k) m = w (k) 0 , . . . , w (k) m (k) -1 (12)
be a basis of the associated Krylov subspace.

As explained in the previous section, augmentation never increases the condition number which governs the asymptotic convergence rate. More precisely, the presence of active eigenvectors of the current preconditioned problem in C (k) may increase the efficiency of the iterative solver significantly. Classical strategies can be used in the case of invariant preconditioned operators (A (k) = A, M (k) = M ) and multiple right-hand sides.

It is more difficult to define efficient strategies in the general case of varying operators with no information available on their evolution. A simple and natural idea is to reuse previous Krylov subspaces. A first algorithm which reuses all the previous Krylov subspaces is Total Reuse of Krylov Subspaces (TRKS) (Algorithm 2), which needs only a few comments:

• Since (according to (8)) C (k) T A (k) W (k)
m = 0, the vectors of the concatenated matrix C (k+1) are linearly independent. Therefore, APCG(A (k+1) , M (k+1) , C (k+1) , b (k+1) ) does not break down and converges.

• The previous Krylov subspaces are fully reused through concatenation without post-processing; the only downside is that the memory requirements increase due to the need to save the Krylov subspaces.

• If the number of columns of matrix C (k) becomes too large, the method may become computationally inefficient, even though the number of iterations decreases considerably. Nevertheless, TRKS probably leads to the best reduction in the number of iterations achievable by reusing Krylov subspaces. Therefore, it can be used as a reference in terms of the reduction of the number of iterations for any other algorithm based on a reuse of Krylov subspaces.

• One possible way to reduce the cost of TRKS without reducing the size of C consists in using approximate solvers, as in the IRKS strategy [START_REF] Risler | Iterative accelerating algorithms with Krylov subspaces for the solution to large-scale nonlinear problems[END_REF].

Algorithm 2: TRKS-APCG Initialize C (0) = C 0 (an n × m 0 full-rank matrix); for k = 0, . . . , p -1 do Solve A (k) x (k) = b (k) with APCG(A (k) , M (k) , C (k) , b (k) ); Define W (k) m = . . . , w (k) j , . . . 0 j<m (k) ; Concatenate: C (k+1) = C (k) , W (k) m end
In order to reduce the cost associated with the total reuse of Krylov subspaces, we propose to work on extracted sub-subspaces, an operation often referred to as the recycling of Krylov subspaces. The objective is to retain the smallest number of independent vectors which achieve the greatest decrease in the number of iterations. Clearly, the most effective approach would be to calculate approximate eigenvectors from the previous Krylov subspaces for the current operator. However, because of the variability of the operators, the extraction of such information would be extremely time consuming and would affect the global efficiency. Conversely, approximate eigenvectors of previous problems can be calculated from the associated Krylov subspaces at nearly no cost. In the following section, we describe an efficient algorithm for the extraction of such approximation vectors along with a simple selection procedure to recycle only a few of these vectors. Of course, the performance of our method depends on the stability of the eigenspaces from one system to another. This topic, especially concerning the lower part of the spectrum, is discussed in [START_REF] Kilmer | Recycling subspace information for diffuse optical tomography[END_REF].

Selective recycling of Krylov subspaces

The standard convergence of conjugate gradients corresponds to an asymptotic convergence rate. Using this property to predict the number of iterations n which is required to reach an accuracy level cg leads to a huge overestimation.

Indeed, one has:

x i -x A x 0 -x A 2(σ 1,n ) i cg ⇒ i n = ln( cg /2) ln(σ 1,n ) with σ r,s = √ κ r,s -1 √ κ r,s + 1 and κ r,s = λ r λ s (13)
This result alone cannot explain the improvement in the convergence rate observed during the iteration process. This superconvergence phenomenon can be explained by a study of the convergence of Ritz values [START_REF] Morgan | A restarted GMRes method augmented with eigenvectors[END_REF] which enables one to define an instantaneous convergence rate [START_REF] Van Der Sluis | The rate of convergence of conjugate gradients[END_REF]. This explanation can be improved by a study of the influence of the distribution of the eigenvalues [START_REF] Notay | On the convergence rate of the conjugate gradients in presence of rounding errors[END_REF][START_REF] Axelsson | On the rate of convergence of the preconditioned conjugate gradient method[END_REF].

The objective of recycling Krylov subspaces is to find the best augmentation space in order to trigger superconvergence quickly. This section is organized as follows: we start with a review of Ritz eigenelement analysis and continue with a brief presentation of the improved convergence results; these results lead to a number of selection strategies, which will be assessed in Section 5.

Ritz analysis: theory and practical calculation

For 0 i < m, Ritz vectors (ŷ i m ) and values (θ i m ) are approximations of the eigenvectors and eigenvalues of the symmetric positive definite matrix Â; their definition is similar to that of the iterates in the conjugate gradient algorithm [START_REF] Carpraux | Stability of the Krylov bases and subspaces[END_REF] find

(ŷ i m , θ i m ) ∈ K m ( Â, v0 ) × R such that Ây i m -θ i m y i m ⊥ K m ( Â, v0 ) (14)
The symmetric Lanczos algorithm [START_REF] Saad | Numerical Methods for Large Eigenvalue Problems[END_REF] enables one to build a particular orthonormal basis of K m ( Â, v0 ), denoted V m . Then, the search principle becomes:

y i m = Vm q i m , V T m  Vm q i m = θ i m q i m (15)
The Lanczos basis Vm makes the Hessenberg matrix Ĥm = V T m  Vm symmetrical and tridiagonal. Vm and Ĥm can be recovered directly from the conjugate gradient coefficients [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]: Vm = . . . , (-1) j rj rj , . . .

0 j<m Ĥm = tridiag(η j-1 , δ j , η j ) 0 j<m with δ 0 = 1 α 0 , δ j = 1 α j + β j-1 α j-1 , η j = β j α j (16)
Since matrix Ĥm is symmetrical and tridiagonal, its eigenelements (θ m j , q m j ) 1 j m can be calculated easily, for example using a Lapack procedure. Let us define We presented Ritz analysis for the equivalent symmetric system described previously because symmetry simplifies the calculation of eigenelements, but the analysis can be transferred back to the left-preconditioned system using the following transformation rules:

Θ m = diag(θ 1 m . . . θ m m ) and Q m = q 1 m , . . . , q m m such that Ĥm = Q m Θ m Q m T . Θ m
V m = L -T Vm = . . . , (-1) j z j (r j , z j ) 1/2 , . . . H m = Ĥm = V T m AV m Y m = L -T Ŷm = V m Q m ( 17
)
The Ritz vectors are the solution of a generalized eigenproblem and satisfy the following orthogonality properties:

Y T m AY m = Θ m and Y T m M Y m = I m (18)
One can show that when m increases the Ritz values converge toward the eigenvalues of Â, and that the convergence is either from above or from below depending on their rank [START_REF] Van Der Sluis | The rate of convergence of conjugate gradients[END_REF][START_REF] Jia | On the convergence of the Ritz values, Ritz vectors and refined Ritz vectors[END_REF]:

θ 1 m θ 1 m-1 θ 2 m . . . θ m-1 m θ m-1 m-1 θ m m (19)
In addition, in the case of clearly distinct eigenvalues, the convergence of a Ritz value results in the convergence of the associated Ritz vector.

Relation between the convergence of conjugate gradients and the convergence of the Ritz values

In [START_REF] Van Der Sluis | The rate of convergence of conjugate gradients[END_REF], the superconvergence phenomenon is explained by the convergence of the Ritz values through the definition, at each iteration, of a instantaneous convergence rate associated with the part of the spectrum that is not yet approximated correctly by the Ritz values: at a given conjugate gradient iteration, one can find a deflated system (with some of its extreme eigenvalues removed) with similar behavior. Let [λ l , .., λ r ] be the spectrum of the deflated operator.

The equivalent convergence rate is:

x -x i+1 A F i,l,r 2 σ l,r x -x i A ( 20
)
where F i,l,r quantifies the convergence of the l smallest and r largest Ritz values to the extreme eigenvalues:

F i,l,r = max l >l J (i) l,l max r r L (i) r,r J (i) l,l = l j=1 1 - λ l λ j 1 - λ l θ i j -1 L (i) r,r = r j=1 1 - λ n-r λ n+1-j 1 - λ n-r θ i i+1-j -1
Since this result holds for every pair (l, r), the effective convergence rate at Iteration i corresponds to the pair (l, r) which minimizes σ i,l,r = F i,l,r σ l,r . Then, after some iterations, the superconvergent conjugate gradient algorithm behaves very much like a conjugate gradient algorithm augmented by the extreme eigenvectors which are associated with the converged Ritz values. In a multiresolution context, provided the linear systems have similar spectral properties, the Ritz vectors associated with the converged Ritz values obtained for one system should define a viable augmentation space for the subsequent resolutions.

Effect of the distribution of the eigenvalues

The effect of the distribution of the eigenvalues on the convergence of conjugate gradients was studied in [START_REF] Notay | On the convergence rate of the conjugate gradients in presence of rounding errors[END_REF][START_REF] Axelsson | On the rate of convergence of the preconditioned conjugate gradient method[END_REF]. The results take into account the fact that preconditioning often leads to clustered eigenvalues as opposed to uniformly distributed eigenvalues, as can be seen in Figure 1.

In addition to other results, the authors showed that if a spectrum consists of p isolated eigenvalues in the high part of the spectrum, p isolated eigenvalues in the low part of the spectrum and n-2p uniformly distributed central eigenvalues, then the conjugate gradient convergence takes the form:

n ñ = 2p + int   ln ( cg /2) ln σ p+1,n-p - p i=1 ln λ n-p+i 4λ i 1 - λ i λ n-p+i ln σ p+1,n-p   (21)
The convergence rate is approximately equal to the classical convergence rate for the central part, plus one iteration per higher eigenvalue and a little more than one iteration per lower eigenvalue. These results can be combined with the work by Jiao [START_REF] Jia | On the convergence of the Ritz values, Ritz vectors and refined Ritz vectors[END_REF][START_REF] Jia | The convergence of harmonic Ritz values, harmonic Ritz vectors and refined harmonic Ritz vectors[END_REF] on the convergence of Ritz values. In general, since the method is related to the power iteration method, a correct approximation by the Ritz values is obtained first for the highest eigenvalues, then for the lowest part of the spectrum, resulting in superconvergence (which is governed by the asymptotic convergence rate of the reduced spectrum).

Selection procedures

The results of Section 4.2 lead to a first proposal of a selection procedure for converged Ritz vectors: convergence is identified by the stagnation of the Ritz values; if the conjugate gradient algorithm converges at iteration m, the Ritz values are calculated for the previous two states Θ m and Θ m-1 . Once ranked, the m most recent Ritz values Θ m are compared to the m -1 previous values according to the following criteria:

           θ j m has converged if |θ j m -θ j m-1 | |θ j m | ε, 1 j m-1 θ m-j m has converged if |θ m-j m -θ m-1-j m-1 | |θ m-j m | ε, 0 j m -2 ( 22
)
where ε is a user parameter which is easy to adjust since the criterion is generally either very high (before the convergence of the Ritz value) or very small (after convergence). Figure 1 illustrates that property with the simple example of the operator associated with the decomposition of a linear elastic cube into ten subdomains; in that case, the higher half of the spectrum has converged. The principle of the selective recycling of Krylov subspaces (SRKS-APCG) is described in Algorithm 3. Basically, in addition to the memory required by APCG, the SRKS-APCG algorithm requires storage for m n-vectors (z j ) j=1,m . One should note that the selected vectors are normalized by the square root of the associated Ritz value in order to improve the condition number of the coarse matrix. (If operator A remained constant, matrix (C T AC) would be the identity matrix.)

For better computational efficiency, a restart parameter can be introduced in order to limit the size of the augmentation space associated with parameter Algorithm 3: SRKS-APCG Initialize C (0) = C 0 (full column rank matrix); for k = 0, . . . , p -1 do

• Solve A (k) x (k) = b (k) with APCG(A (k) , M (k) , C (k) , b (k) );
• Define V m = . . . , (-1) j zj (rj ,zj ) 1/2 , . . .

0 j<m ; • Define H m = tridiag(η j-1 , δ j , η j ) 0 j<m δ 0 = 1 α0 , δ j = 1 αj + βj-1 αj-1 , η j = √ βj αj ; • Compute eigenelements (Q m , Θ m ) of H m (θ 1 m . . . θ m m ); • Compute Y m = V m Q m = y 1 m , . . . , y m m ; • Extract H m-1 = tridiag(η j-1 , δ j , η j ) 0 j<m-1 ; • Compute eigenvalues (θ j m-1 ) of H m-1 ; for j = 1, . . . , m -1 do C = C, y j m √ |θ j m | if |θ j m -θ j m-1 | ε|θ j m |; C = C, y j+1 m √ |θ j+1 m | if |θ j+1 m -θ j m-1 | ε|θ j+1 m | end • Concatenate C (k+1) = [C (k) ; C] , C = [0]; • If dim(C (k+1) ) n c lim , then C (k) = C (0)
end n c lim in Algorithm 3. This limit size can be set after a complexity analysis under the assumption that all non-augmented systems would be solved in the same number of iterations. However, we did not use such a restart procedure in our experiments.

In order to be even more selective, we propose a reselection strategy based on a prediction of the efficiency of the retained vectors. Indeed, the results of Section 4.3 in terms of the effect of the distribution of the eigenvalues lead us to retain only the converged Ritz vectors which belong to the external part of the spectrum:

• this is known to be the first part of the spectrum whose approximation by Ritz values is good;

• since the convergence of Ritz vectors is identified by the stagnation of the associated Ritz values, the fact that the external Ritz values are distinct ensures that the Ritz vectors approximate the eigenvectors correctly [START_REF] Jia | On the convergence of the Ritz values, Ritz vectors and refined Ritz vectors[END_REF];

• while choosing vectors in the dense central zone does not modify the shape of the spectrum and does not improve convergence, selecting the external part of the spectrum triggers superconvergence instantly.

In order to select only the external part of the spectrum, we implemented the cluster identification algorithm proposed in [START_REF] Molinari | Multiple temporal cluster detection[END_REF]. This algorithm seeks the piecewise constant distribution which is nearest (in a least squares sense) to the distribution of the distances among the sorted eigenvalues. The only parameter required is the minimum size of the cluster, which we set at one-fifth the number of preselected vectors. As will be shown in the next section, the performance achieved with this reselection algorithm is not outstanding, but some results in terms of gain per augmentation vector are worth considering.

Numerical assessments

We present three numerical experiments. Two concern the evaluation of a structure made of random materials, as is the case in a Monte-Carlo simulation. In the first case, the materials are elastic; in the second case, which is a nonlinear problem, they are elastic-plastic. The last case is a large displacement problem, which raises specific difficulties.

The methods were implemented in the ZEBULON code [START_REF]Z-set user manual[END_REF] and parallelism was introduced using MPI. The calculations were performed on the LMT-Cachan cluster, which consists of dual quadcore and dual hexacore processors connected by a gigabit network. The calculations were always carried out on homogeneous sets of processors which were entirely dedicated to one task which fit entirely in memory, so swapping was not necessary. In each case, we indicate the CPU time which measures the amount of work performed for one subdomain. The Wall Clock Time (WCT), a global measure which is more sensitive to external perturbations induced by the operating system and the presence of other users, was considered to be unreliable in many cases; so we mention it only for the first set of experiments. One should note that the gains calculated with WCT were always greater.

The CPU plots show the total time as well as the time dedicated to augmentation (preparation of the coarse operator, initialization and projections); the difference represents the iterations of the solver.

All the calculations used a dual formulation of the interface problem through domain decomposition (FETI). The convergence was evaluated using the norm of the residual (which corresponds to the displacement gaps at the interfaces) normalized by the condensed right-hand side. Classically for such structural problems, total reorthogonalization was used to enforce the A-conjugation of the search directions. (The case without reorthogonalization is discussed briefly in the first example.)

The case of a sequence of linear systems

We considered a cube (of side 50 mm) with 4 × 4 × 4 = 64 small cubic inclusions (of side 5.5 mm). A slice through this structure is shown in Figure 2. The cube was clamped over one side, and the opposite side plus another side were subjected to uniform pressure. The mesh consisted of 125, 000 linear hexahedral elements for a total of 400, 000 degrees of freedom. Three automatic decompositions (into 12, 48 and 96 subdomains) were performed using the Metis algorithm [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF] (see Figure 3). The resulting interface system contained 54, 000 unknowns for the 12-subdomain decomposition, 96, 000 unknowns for the 48-subdomain decomposition and 133, 000 unknowns for the 96-subdomain decomposition. All the materials were isotropic, linear and elastic, and were characterized by their Young's modulus and Poisson's coefficient. The material properties of each inclusion and of the matrix were chosen randomly following a normal law with a relative standard deviation equal to 10%, leading to a ±23% variation range about the nominal value. The average Young's modulus was 200 MPa for the matrix and 20, 000 MPa for each inclusion, and the average Poisson's coefficient was 0.27 for the matrix and 0.35 for each inclusion. The objective was to perform the calculations for 40 draws of the 130 coefficients. We used a dual formulation (FETI) with both a Dirichlet (optimal) and a lumped preconditioner, leading to 10 -3 and 10 -6 APCG accuracy respectively. We considered the following algorithms: conjugate gradients (cg), total reuse of Krylov subspaces (trks) and selective reuse of Krylov subspaces with two values of the criterion, ε = 10 -6 (srks6) and ε = 10 -14 (srks14). In addition, in the last case (ε = 10 -14 ), we also attempted to further refine the selection by not selecting the converged Ritz values contained in the central cluster (identified by the algorithm proposed by [START_REF] Molinari | Multiple temporal cluster detection[END_REF]); this method is labeled (clust14).

Comparison of the strategies

The results for the 12-subdomain decomposition are summarized in Table 1, which gives the average number of APCG iterations to convergence, the average size of the augmentation space, the final size of the augmentation space, the average CPU and wall clock times per system, from which we also deduced the augmentation time (operator preparation and projection). Computations were conducted one dual hexacore processor (one subdomain per core). When the average and final sizes of the augmentation space are close, this means that For the 12-subdomain decomposition, Figures [START_REF] Carpraux | Stability of the Krylov bases and subspaces[END_REF][START_REF]Conjugate gradient method with preconditioning by projector[END_REF][START_REF] Erhel | An augmented conjugate gradient method for solving consecutive symmetric positive definite linear systems[END_REF] (for an objective of 10 -3 accuracy) and Figures [START_REF] Chapman | Deflated and augmented Krylov subspace techniques[END_REF][START_REF] Erhel | Restarted GMRes preconditioned by deflation[END_REF][START_REF] Farhat | The two-level FETI methodpart II: Extension to shell problems. parallel implementation and performance results[END_REF]) (for an objective of 10 -6 accuracy) give the evolutions of the number of APCG iterations to convergence for each linear system, the dimension of the augmentation space n c and the the CPU time for the resolution of each system, with both lumped and Dirichlet preconditioners.

Without full reorthogonalization, the performance was very poor and led to about twice the number of iterations of the recommended fully reorthogonalized conjugate gradients. This was expected because systems resulting from domain decomposition formulations are known to often require full reorthogonalization [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF]. Furthermore, one should note that the additional iterations carried out in the non-reorthogonalized case led to vector sets which made the Ritz analysis more complex due to the appearance of nonphysical, multiple eigenvalues. The non-reorthogonalized approach was no longer considered in the other examples.

With the TRKS approach, two types of behavior were observed. In the lowaccuracy case (10 -3 ), for both preconditioners (but especially for the lumped preconditioner), the size of the augmentation space reached a plateau, which means that the augmentation space contained almost all the required information; the gains in terms of both the number of iterations (> 90%) and the CPU time (> 55%) were excellent. In the high-accuracy case (10 -6 ), the size of the augmentation space never stabilized; therefore, even though the number of iterations decreased drastically, the CPU time increased. Table 2 gives extended performance results for TRKS which confirm this analysis. The gains are given relative to conjugate gradients. The efficiency of augmentation is defined by the average decrease in the number of iterations per augmentation vector; the higher the required accuracy, the less efficient the TRKS approach. These results justify our decision to select the subspaces so that the dimension of the augmentation space would remain under control.

The SRKS14 approach succeeded in limiting the size of the augmentation space and led to a satisfactory decrease in the number of iterations. As can be seen on the figures, SRKS6 did not stabilize the augmentation space as efficiently and behaved half way between TRKS and SRKS14; therefore, we will choose SRKS14 as our reference algorithm from now on.

The cluster strategy as it stands today gave unsatisfactory results: even though it often led to the best gain per augmentation vector, it seemed to impair the selection of useful vectors and allow much less reduction in the number of iterations than SRKS. After the resolution of many systems, it tended to lead to the same augmentation space as SRKS.

To confirm that hypothesis, we compared the spaces C SRKS and C cluster after the 40 resolutions for the low-accuracy Dirichlet case. We used the following procedure: first, the vectors were orthonormalized using SVD: C = U ΣV T ; then SVD was applied to the concatenated matrix [U SRKS , U cluster ]. A plot of the singular values is shown in Figure 10. Independent spaces would lead to a constant value equal to 1, while for nested spaces the common space would lead to { √ 2, 0} pairs of singular values. One can observe that the spaces are not exactly nested, but come quite close.

In conclusion, the cluster strategy is not mature yet, but it is promising. It was not considered for the following experiments because, due to the larger number of systems involved, it would behave quite similarly to SRKS. 3 shows the relative performance of SRKS14 as a function of the number of subdomains, of the preconditioners and of the accuracy. The efficiency of the augmentation is defined as the decrease in the average number of iterations per augmentation vector. One can observe that the efficiency ranged between 0.5 and 0.85 and was best for the lower accuracy and the improved preconditioner. For these spectra in which there exist no small isolated eigenvalues (which could lead to efficiencies greater than 1), such results are consistent with the theory (see Section 4.3). In the next section, we will see that this moderate efficiency does not preclude significant CPU improvements.

Study of SRKS14 in various configurations

The gains in terms of the number of iterations were relatively stable, typically between 50% and 60% in the high-accuracy case.

Influence of the hardware configuration on the CPU gains

Now, let us study the performance of SRKS14 in terms of CPU time for the same decomposition into 48 subdomains, but using different hardware configurations: One can note that the processors in Configuration A were different from those used in the other cases. In all the cases, the memory was sufficient to avoid swapping. The results are given in Table 4 for the Dirichlet preconditioner and in Table 5 for the lumped preconditioner. One can see that Configurations B,C and D had similar performances and were slower than Configuration A due to the different memory technology. One interesting factor is the ratio of the average CPU cost of an iteration to the average CPU cost of an augmentation vector (the last columns of Table 4 and5). One can see that in Configuration A, 4 augmentation vectors cost no more than one iteration; in the other configurations 7 augmentation vectors cost no more than one iteration. Since we saw that one needs about 1/0.6 1.6 augmentation vectors to save one iteration, the advantage of augmentation is clear. Indeed, we observe a 32% CPU improvement in Configuration A and a 40% to 50% improvement in the other configurations.

Note that when the lumped preconditioner is used the equivalent cost of an iteration is only 2.8 augmentation vectors in Configuration A and 4.5 augmentation vectors in Configuration D (see Table 5). Since the efficiency of the augmentation vectors is less when this inexpensive preconditioner is used (in the high accuracy case), so is the CPU improvement. The ratio of the CPU time per iteration to the CPU time per augmentation vector for SRKS (Column 4 of the previous tables) turned out to be relatively stable for a given machine with a given preconditioner. This is due to the stability of the size of the augmentation space which prevented the cost from soaring (as would happen with TRKS). Thus, the CPU performance can be deduced from the iteration gains and the augmentation efficiency (see Table 3). For instance, the CPU gain for SRKS with the 96-subdomain decomposition was slightly greater than 50%.

The case of a sequence of nonlinear problems

Now let us consider a hexahedral holed plate (10×10×0.2 mm with a center hole of radius 1 mm, see Figure 11) subjected to unidirectional tension (a prescribed normal displacement). The plate was discretized into 61, 000 linear hexahedral elements for a total of 41, 000 degrees of freedom. The structure was divided into 8 subdomains using the Metis algorithm, which resulted in an interface system with 3, 000 unknowns. The problem was solved using one 8-core processor (one subdomain per core). Elastic-plastic behavior with nonlinear isotropic hardening and a Von Mises'-type plasticity criterion was assumed. Denoting σ the Cauchy stress tensor, (u) the symmetric gradient of the displacement field u, and K the Hooke tensor, the material law can be written as:

         (u) = e + p , σ = K : e if f (σ) = 0 then ˙ p = λf ,σ if f (σ) 0 then ˙ p = 0 f (σ) = 3 2 σ : σ -R 0 + Q 1 -e -bλ (23)
The coefficients were assigned a normal law with a 10% relative standard deviation, which implied variations of up to ±23% in the coefficients. The mean values of the material parameters were: E = 200, 000 MPa, ν = 0.3, R 0 = 300 MPa, b = 22 and Q = 170 MPa. The loading was applied in two steps: first, a single increment to reach the elastic limit; then, 16 equal increments in order to multiply the prescribed displacement by 4. The objective of the study was to analyze 21 configurations. Again, the linear solver used was FETI with a Dirichlet or lumped preconditioner. The accuracy objective for the linear systems was set at 10 -6 . (The accuracy must be high for the nonlinear process to run well). Because of the approximations, not all the methods converged in the same number of Newton iterations; on average, one nonlinear analysis required the resolution of 95 tangent systems. Table 6 summarizes the performances of the various methods; Figure 12 shows the evolution of the average number of APCG iterations with the lumped preconditioner during the sequence of linear systems; Figure 13 shows the evolution of the size of the augmentation space; Figures 14 and15 show the evolutions of the average CPU time and wall clock time for the resolution of one linear system along with the evolution of the average augmentation time (operator creation and projection). This test leads to conclusions similar to the previous ones. More specifically, one can observe that the SRKS augmentation space selected after the first nonlinear configuration remained stable. Conversely, since the TRKS augmentation space never reached a plateau, the solutions of the linear systems did not belong to a common space. SRKS was the most efficient method, leading to a 20% CPU gain and a 36% wall clock time improvement.

The case of a large displacement problem

Finally, let us consider the problem of the buckling of a straight heterogeneous beam with a circular cross section (length/diameter ratio equal to 30), clamped at one end and subjected to an axial pressure at the other, with no radial displacement. The heterogeneities consisted of five straight fibers whose stiffness was 1, 000 times that of the matrix. The problem was formulated in the updated Lagrangian framework, assuming linear elastic behavior (characterized by the Young's modulus and Poisson's coefficient) in the current configuration. The beam was discretized into 90, 000 linear hexahedral finite elements for a total of 300, 000 degrees of freedom. It was divided into 10 subdomains using the Metis algorithm, leading to an interface system with 16, 000 unknowns. A single 12core processor was used (1 subdomain per core, leaving 2 inactive cores). The pressure was applied incrementally up to the configuration shown in Figure 18, in which the maximum axial displacement was about 3% of the total length. 12 increments were used, leading to the resolution of about 30 tangent linear systems.

We used a FETI solver with a Dirichlet preconditioner and an "identity" projector. The FETI convergence criterion was set to 10 -6 . Figure 16 shows the evolution of the number of conjugate gradient iterations required for the resolution of each linear system. Figure 17 shows the evolution of the size of the augmentation space. Three algorithms were tested: classical conjugate gradients, total reuse of subspaces, and selective reuse of subspaces (ε = 10 -14 ). Table 7 summarizes the main results.

The following observations can be made: • The performance of TRKS was really impressive because the systems converged in 10 times fewer iterations than with CG, but the size of the associated augmentation space was large (up to 208 vectors) and never ceased to increase;

• SRKS appeared to be efficient: the number of iterations was divided by 3 with a space whose size increased slowly, then reached a plateau;

• With the hardware configuration used, the best results in terms of computation time (a 48% CPU improvement) were achieved with TRKS, but the gain normalized by the number of augmentation vectors was better with SRKS (a 29% CPU improvement). SRKS was truly successful in controlling the dimension of the augmentation space.

In this example, the augmentation proved to be very efficient in terms of the iteration gain per augmentation vector, especially for SRKS (0.85). This was probably because of a specificity of the spectrum of the preconditioned operator due to the use of domain decomposition in large displacements. Indeed, the tangent matrix of a floating subdomain with prescribed [START_REF] Paige | Approximate solutions and eigenvalue bounds from Krylov subspaces[END_REF], but convergence is slower than when all the eigenvalues are positive. The positivity of the preconditioner makes the selection procedure still possible. Moreover, the negative eigenvalues are systematically selected by the procedure, so using augmentation causes the solver to iterate in the subspace in which the operator is positive, leading to a much better convergence rate.

Conclusion

This paper dealt with the resolution of sequences of large linear systems with varying matrices and right-hand sides using conjugate gradients. We proposed several algorithms based on the augmentation of the current Krylov subspace by a selection of previously generated subspaces. The advantage of these methods is that some of the iterations are replaced by the preprocessing of a coarse problem associated with optimized operations. When low accuracy is sufficient, total reuse of the previous subspaces (the TRKS algorithm) appears to lead to satisfactory results. When high accuracy is required, the subspaces are too unstable, which causes the dimension of the TRKS augmentation space to soar. Therefore, we proposed to retain only the part of the subspace generated by the Ritz vectors associated with converged Ritz values of the preconditioned operator (the SRKS algorithm). These vectors can be built very inexpensively. Such an augmentation was found to remain stable throughout the linear systems and to lead to a reduction in the number of iterations which is consistent with the theory. In terms of computation time, the proposed method leads to a variable, but always positive, gain compared to non-augmented systems. We observed CPU time improvements of 20% to 50%, and wall clock time improvements of 40% to 70%.

Up until now, our attempts to improve the selection algorithm by eliminating the converged values in the central part of the spectrum have not led to impressive results. This probably means that the Ritz vectors associated with converged Ritz values contain meaningful information which cannot be removed from the analysis of the current system. A continuation of this work could consist in a better analysis of the accumulation of the augmentation vectors. This can be done by studying the coarse matrix C T AC whose distance to the identity matrix characterizes the variation of the Krylov subspaces. Another objective would be to port some of the ideas presented in this paper to nonsymmetric solvers.

  and Ŷm = Vm Q m are the Ritz values and associated Ritz vectors, which are approximations of the eigenelements of operator  and satisfy: Ŷ T m  Ŷm = Θ m and Ŷ T m Ŷm = I m

2 Figure 1 :

 21 Figure 1: Ritz spectrum and convergence of the Ritz values

Figure 2 :Figure 3 :

 23 Figure 2: A slice through the heterogeneous cube (shear stress) Figure 3: The decomposition into 48 subdomains

Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 4567 Figure 4: Cube 12 subdomains, lumped, 10 -3 accuracy, number of iterations per linear system Figure 5: Cube 12 subdomains, Dirichlet, 10 -6 accuracy, number of iterations per linear system

Figure 8 :Figure 9 :

 89 Figure 8: Cube 12 subdomains, lumped, 10 -3 accuracy, CPU time per linear system Figure 9: Cube 12 subdomains, Dirichlet, 10 -6 accuracy, CPU time per linear system

Figure 10 :

 10 Figure 10: Singular values of [C SRKS , C cluster ] accur. precond. # subdomains iteration gain CPU gain efficiency of augmentation

Figure 11

 11 Figure 11: The holed plate example (plastic strain)

Figure 12 :Figure 13 :Figure 14 :Figure 15 :

 12131415 Figure 12: Plate, lumped -avg. # it. / linear system Figure 13: Plate, lumped -dimension of aug. space

Figure 16 :Figure 17 :Figure 18 :

 161718 Figure 16: Buckling of the heterogeneous beam: number of iterations per linear system Figure 17: Buckling of the heterogeneous beam: dimension of the augmentation space for each linear system

Table 1 :

 1 Performance summary for the cube with inclusions most of the augmentation space was identified with the first systems. For a given configuration (accuracy and preconditioner), the figures in bold in the three columns 'average number of iterations', 'average CPU time' and 'average wall clock time' indicate the best strategy in terms of gain per unit augmentation vector compared to CG.

	12 subdomains

Table 2 :

 2 Relative performance of TRKS

	-3	lumped	12 48	99.4% 97.7%	58.2% 60.1%	0.94 0.97
		Dirichlet	12 48	94.7% 96.2%	55.5% 62.6%	0.54 0.69
	10 -6	lumped Dirichlet	12 48 12 48	89% 90.2% 84.4% 83.3%	-121.2% -162.3% -55.1% -90.7%	0.25 0.37 0.19 0.21
	Table

Table 4 :

 4 CPU performance of SRKS14 for 10 -6 accuracy with the Dirichlet preconditioner

	Configuration CG avg. CPU CPU gain CPU per iteration /
				CPU per augm. vector
	A	9.5	32.6%	4
	B	25.7	41.5%	6.7
	C	29.7	48.9%	7.7
	D	29.9	47.6%	7.8
	Configuration CG avg. CPU CPU gain CPU per iteration /
				CPU per augm. vector
	A	10.7	22.4%	2.8
	D	27.4	33.3%	4.5

Table 5 :

 5 CPU performance of SRKS14 for 10 -6 accuracy with the lumped preconditioner

Table 6 :

 6 Holed plate, performance summary

Table 7 :

 7 Recycling performance for the buckling problem may become non-positive, contrary to the Dirichlet operator which remains positive definite. It is known that a slight lack of positivity of the operator does not prevent reorthogonalized conjugate gradients from converging

	Neumann conditions

"active" referring to the part of the spectrum of the matrix which is solicited by the right-hand side.