# Enrichissement de préconditionneurs pour la résolution de systèmes linéaires symétriques indéfinis

Sylvain Mercier, Serge Gratton, Nicolas Tardieu, Xavier Vasseur

► To cite this version:

Sylvain Mercier, Serge Gratton, Nicolas Tardieu, Xavier Vasseur. Enrichissement de préconditionneurs pour la résolution de systèmes linéaires symétriques indéfinis. 12e Colloque national en calcul des structures, CSMA, May 2015, Giens, France. hal-01515057

HAL Id: hal-01515057

https://hal.science/hal-01515057v1

Submitted on 27 Apr 2017

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

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

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

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Enrichissement de préconditionneurs pour la résolution de systèmes linéaires symétriques indéfinis

S. Mercier<fnref n="1" />,<fnref n="2" />, S. Gratton<fnref n="3" />, N. Tardieu<fnref n="1" />, X. Vasseur<fnref n="2" />

**Résumé** — Dans cette communication est présentée une méthode d'enrichissement de préconditionneurs dans le cadre des méthodes itératives de Krylov pour une séquence de systèmes linéaires symétriques indéfinis. En supposant que les matrices varient légèrement, l'idée est de récupérer de l'information spectrale approchée lors de la résolution du premier système afin d'accélérer la convergence des suivants. Cette méthode est illustrée par deux études industrielles d'EDF nécessitant la résolution de problèmes de type point-selle, effectuées à l'aide du logiciel *Code_Aster*.

**Mots clés** — Préconditionnement, système symétrique indéfini, point-selle.

## 1 Introduction

En mécanique des solides, différents types de non-linéarités peuvent intervenir provenant par exemple du comportement des matériaux ou de grandes déformations géométriques. Pour résoudre ceci dans un cadre implicite, il est courant d'utiliser une méthode de type Newton qui entraîne la résolution d'une séquence de systèmes linéaires. De plus, la prise en compte de relations linéaires imposées à l'aide de multiplicateurs de Lagrange induit souvent une structure de point-selle aux systèmes, symétriques mais non-définis. Dans l'optique d'améliorer l'efficacité des méthodes itératives de Krylov dans ce cadre (cf. [7]), nous avons travaillé sur l'enrichissement de préconditionneur en généralisant la famille de "Limited Memory Preconditioners" (LMP) développée dans le cas symétrique défini positif dans [3] et [8].

## 2 "Limited Memory Preconditioners" (LMP) pour systèmes symétriques indéfinis

Afin d'alléger les notations, nous considérons dans cette section un seul système linéaire à résoudre :

$$Ax = b \quad \text{avec } A \in \mathbb{R}^{n \times n}, \quad x \text{ et } b \in \mathbb{R}^n.$$

Supposons que l'on bénéficie déjà d'un préconditioneur *M*, dit de *premier niveau*, permettant d'accélérer la convergence de la méthode de Krylov choisie. Afin d'améliorer l'effet de celui-ci, un préconditioneur de *second niveau* appelé LMP a été étudié dans [6] puis dans [3] et [8] dans le cas symétrique défini positif (SDP). L'idée principale est de récupérer de l'information sur certaines directions spectrales qui ralentissent la convergence afin de réduire leur influence.

**Définition 2.1** Soit *A* et *M* un SDP d'ordre *n*. Supposons que *S* est une matrice *n* par *k* de rang *k*, avec *k* ≤ *n*. La matrice symétrique suivante est appelée LMP :

$$H = [I_n - S(S^T A S)^{-1} S^T A] M^{-1} [I_n - A S(S^T A S)^{-1} S^T] + S(S^T A S)^{-1} S^T.$$

Ce préconditioneur *H* permet notamment de transformer les vecteurs colonnes de *S* en vecteurs propres de *HA* associés à 1. De plus, le reste du spectre de *HA* est contracté par rapport à celui de $M^{-1}A$ (cf. théorème 3.4 dans [3]). L'objectif de cette étude est de généraliser cet effet au cas symétrique indéfini.

<footnote n="1">¹ EDF R&D, 1 av. du Général De Gaulle, 92141 Clamart, France, {sylvain-s.mercier,nicolas.tardieu}@edf.fr</footnote>

<footnote n="2">² CERFACS, 42 av Gaspard Coriolis, 31057 Toulouse Cedex 1, France, xavier.vasseur@cerfacs.fr</footnote>

<footnote n="3">³ INPT-IRIT, Université de Toulouse et ENSEEIHT, 2 rue Camichel 31071 Toulouse Cedex 7, France, serge.gratton@enseeiht.fr</footnote>
## 2.1 Généralisation au cas symétrique indéfini

Supposons maintenant que A soit symétrique indéfinie. L'idée est de travailler séparément dans les invariants positif et négatif de A, notés $I_+(A)$ et $I_-(A)$. Ces derniers sont définis comme étant les espaces engendrés par les vecteurs propres associés respectivement aux valeurs propres positives et négatives. Sous certaines hypothèses, nous avons obtenu une généralisation du LMP et démontré un résultat de contraction spectrale.

**Définition 2.2** Soit A symétrique indéfinie et S une matrice n par k de rang k (avec $k \le n$). Supposons que M est SDP avec $M = LL^T$ sa décomposition de Cholesky et notons $K = L^{-1}AL^{-T}$. Si $S = [S_+, S_-]$ où les vecteurs colonnes de $S_+$ (resp. $S_-$) appartiennent à $I_+(K)$ (resp. $I_-(K)$), on définit la matrice symétrique suivante comme une généralisation du LMP :

$$H = [I_n - S(S^T KS)^{-1}S^T K][I_n - KS(S^T KS)^{-1}S^T] + S(S^T KS)^{-1}S^T. \quad (1)$$

**Théorème 2.1** Sous les hypothèses de 2.2, les k vecteurs colonnes de S sont des vecteurs propres de HK associés à 1. De plus, les autres valeurs propres sont contractées par rapport à celles de K de chaque côté de l'axe imaginaire.

Le point faible de cette méthode vient de l'hypothèse sur S. Pour avoir les vecteurs colonnes dans les invariants positifs ou négatifs de K, il faut utiliser la notion de matrice signe de K (décrite dans [4]) qui nécessite un coût de calcul très important. Cependant, une méthode de Krylov comme GMRES permet de récupérer des vecteurs propres approchés (appelés vecteurs de Ritz). On peut montrer que ces vecteurs sont proches de leurs projections dans $I_+(K)$ ou $I_-(K)$.

## 2.2 Application du LMP

Soit une séquence de systèmes linéaires à résoudre, issue par exemple d'une méthode de Newton :

$$A_i x_i = b_i \quad i = 1, \dots, I \quad \text{avec } A_i^T = A_i \in \mathbb{R}^{n \times n}, \quad x_i \text{ et } b_i \in \mathbb{R}^n.$$

Supposons de plus que $A_1 \approx A_2 \approx \dots \approx A_I$ et que l'on dispose d'un préconditionneur de premier niveau $M = LL^T$ SDP pour l'ensemble de la séquence. Nous détaillons ci-dessous les différentes étapes de la résolution.

- Résoudre $A_1x_1 = b_1$ avec GMRES(30) :

$$A_1 x_1 = b_1 \iff \left\{ \begin{array}{l} \overbrace{L^{-1} A_1 L^{-T}}^{K_1} \tilde{x}_1 = L^{-1} b_1 \\ x_1 = L^{-T} \tilde{x}_1 \end{array} \right.$$

- Calculer et choisir k vecteurs de Ritz issus de l'étape précédente. On constitue ainsi S puis H à partir de (1) en prenant $K = K_1$.

- Résoudre la suite de la séquence avec GMRES(30) en appliquant le LMP constitué auparavant :

$$A_i x_i = b_i \iff \left\{ \begin{array}{l} \overbrace{L^{-1} A_i L^{-T}}^{K_i} H \tilde{x}_i = L^{-1} b_i, \quad i = 2, \dots, I. \\ x_i = L^{-T} H \tilde{x}_i \end{array} \right.$$

# 3 Résultats numériques

Nous avons implémenté cette méthode dans le logiciel *Code_Aster* à l'aide de la bibliothèque PETSc, spécialisée dans l'utilisation de méthodes itératives. Les calculs présentés ont été effectués dans un environnement séquentiel.

On peut retrouver en mécanique des structures des problèmes de type point-selle lorsque l'on dualise par exemple les conditions de Dirichlet. Il faut ainsi résoudre des séquences sous la forme suivante :

$$A_i x_i = b_i \iff \begin{pmatrix} G_i & B^T \\ B & 0 \end{pmatrix} \begin{pmatrix} u_i \\ v_i \end{pmatrix} = \begin{pmatrix} f_i \\ g_i \end{pmatrix}, \quad i = 1, \dots, I.$$
Le préconditionnement de ce type de matrice est largement discuté dans la littérature, comme dans [2] ou plus généralement dans [1]. Cependant, d'après les hypothèses énoncées plus haut, il est nécessaire d'appliquer un préconditionneur de *premier niveau* SDP. Nous avons choisi d'utiliser une factorisation de Cholesky incomplète d'un préconditionneur bloc diagonal défini étudié dans [5] :

$$ M = \begin{pmatrix} G_1 + \frac{1}{\gamma} B^T B & 0 \\ 0 & \gamma I \end{pmatrix}, \gamma \in \mathbb{R}^+. $$

## 3.1 Enceinte de confinement d'un réacteur nucléaire

La première étude industrielle concerne l'enceinte de confinement d'un réacteur nucléaire, soumise à des efforts extérieurs et prenant en compte la modélisation des câbles de précontrainte (cf. figures 1 et 2). Ce calcul est linéaire mais le chargement est appliqué graduellement en 4 pas, ce qui donne une séquence à seconds membres multiples. La taille des systèmes est de l'ordre de 400000 inconnues, dont un tiers de multiplicateurs de Lagrange. Le tableau 1 illustre les gains cumulés obtenus grâce au LMP, pour différentes valeurs de $k$.

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Maillage de l'enceinte</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Câbles de précontrainte</figcaption></figure>

<table><caption>TABLE 1 – Statistiques de la résolution pour les 4 systèmes cumulés</caption><thead><tr><th></th><th>Sans LMP</th><th>LMP, k = 5</th><th>LMP, k = 20</th><th>LMP, k = 30</th></tr></thead><tbody><tr><td>Nombre total d'itérations</td><td>675</td><td>555</td><td>509</td><td>438</td></tr><tr><td>Gain en itérations (%)</td><td>×</td><td>18</td><td>25</td><td>35</td></tr><tr><td>Flops<sup>2</sup>(×10<sup>12</sup>)</td><td>6.138</td><td>5.25</td><td>4.901</td><td>4.266</td></tr><tr><td>Gain en Flops (%)</td><td>×</td><td>14</td><td>20</td><td>30</td></tr><tr><td>Mémoire (Mo)</td><td>6436</td><td>6474</td><td>6575</td><td>6643</td></tr><tr><td>Perte en mémoire (%)</td><td>×</td><td>0.5</td><td>2</td><td>3</td></tr></tbody></table>

## 3.2 Agrégats polycristallins

Cette étude est utilisée notamment pour retrouver des lois macroscopiques à partir de considérations microscopiques. La représentation se fait ainsi à une échelle intermédiaire sur un volume élémentaire représentatif (VER). On peut voir sur les figures 3 et 4 le maillage complet ainsi que quelques grains. Le calcul est cette fois non-linéaire (premiers et seconds membres multiples) et de taille de l'ordre de 3 millions d'inconnues, dont 1% de Lagrange. Le tableau 2 présente les performances de la résolution des 4 systèmes linéaires correspondants aux itérations de Newton.

<footnote n="1">1. Pour le choix de $\gamma$, on utilise un coefficient calculé dans *Code_Aster* qui permet de mettre à l'échelle la partie lagrangienne par rapport à la partie physique. Dans ce cas, $G_1$ et $\frac{1}{\gamma}B^T B$ ont le même ordre de grandeur.</footnote>

<footnote n="2">2. 1 flop = 1 opération élémentaire de type (multiplication/division/addition/soustraction), valeur obtenue par PETSc</footnote>
<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Maillage complet</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Vue éclatée de quelques grains</figcaption></figure>

<table><caption>TABLE 2 – Statistiques de la résolution pour les 4 systèmes cumulés</caption><thead><tr><th></th><th>Sans LMP</th><th>LMP, k = 5</th><th>LMP, k = 20</th><th>LMP, k = 30</th></tr></thead><tbody><tr><td>Nombre total d'itérations</td><td>1865</td><td>1605</td><td>1558</td><td>1547</td></tr><tr><td>Gain en itérations (%)</td><td>×</td><td>14</td><td>17</td><td>17</td></tr><tr><td>Flops (×10<sup>11</sup>)</td><td>7.687</td><td>6.734</td><td>6.753</td><td>6.915</td></tr><tr><td>Gain en Flops (%)</td><td>×</td><td>12</td><td>12</td><td>10</td></tr><tr><td>Mémoire (Mo)</td><td>74594</td><td>74871</td><td>75626</td><td>76129</td></tr><tr><td>Perte en mémoire (%)</td><td>×</td><td>0.3</td><td>1</td><td>2</td></tr></tbody></table>

# 4 Conclusions

Un préconditionneur de second-niveau a été développé dans le cadre d'une séquence de systèmes linéaires symétriques indéfinis. L'idée générale est de récupérer de l'information spectrale à l'issue de la première résolution et de l'utiliser pour accélérer la convergence du reste de la séquence. Les résultats obtenus numériquement sur des études industrielles présentent des gains allant jusqu'à 30% en termes d'opérations flottantes sans pour autant augmenter significativement la mémoire nécessaire. Cependant, l'hypothèse de symétrie définie positivité du préconditionneur de premier niveau paraît assez restrictive. Des travaux sont en cours pour généraliser cette méthode au cas non-symétrique, permettant notamment de s'affranchir de cette hypothèse et autorisant ainsi l'utilisation de nombreux préconditionnements de premier niveau.

# Références

[1] M. Benzi, G.H. Golub, J. Liesen. *Numerical solution of saddle-point problems*, Acta Numerica 14, no. 1, 1-137, 2005.

[2] P. Krzyżanowski. *On block preconditioners for saddle point problems with singular or indefinite (1, 1) block*, Numerical Linear Algebra with Applications 18, no.1, 123-140, 2011.

[3] S. Gratton, A. Sartenaer, J. Tshimanga. *On a class of limited memory preconditioners for large scale linear systems with multiple right-hand sides*, SIAM Journal on Optimization 21, no.3, 912-935, 2011.

[4] N.J. Higham. *Functions of matrices : theory and computation*, SIAM, 2008.

[5] TZ. Huang, GH. Cheng, L. Li. *New block triangular preconditioners for saddle point linear systems with highly singular (1, 1) blocks*, Mathematical Problems in Engineering, 2009, 2009.

[6] J.L. Morales, J. Nocedal. *Automatic preconditioning by limited memory quasi-Newton updating*, SIAM Journal on Optimization 10, no.4, 1079-1096, 2000.

[7] Y. Saad. *Iterative methods for sparse linear systems*, SIAM, 2003.

[8] J. Tshimanga. *On a class of limited memory preconditioners for large-scale nonlinear least-squares problems*, PhD thesis, University of Namur, Belgium, 2007.