# Localisation robuste et dénombrement de valeurs propres

## Robust localization and counting of eigenvalues

L. B. NGUENANG * — E. KAMGNIA * — B. PHILIPPE **

* Département d'Informatique, Université de Yaoundé I, Yaoundé, Cameroun.
Emails: lguenang@yahoo.fr, erkamgnia@yahoo.fr.

** INRIA Rennes Bretagne Atlantique, Campus de Beaulieu, 35042 Rennes cedex, France.
Email: Bernard.Philippe@inria.fr

Ce travail a été réalisé dans l'équipe MOMAPPLI du LIRIMA.
This work was pursued within the team MOMAPPLI of the LIRIMA.

LIRIMA=Laboratoire International de Recherche en Informatique et Mathématiques en Afrique,
http://lirima.org/

**RÉSUMÉ.** L'article se consacre à la localisation de valeurs propres pour une grande matrice creuse, a priori non symétrique, dans un domaine du plan complexe. Il combine deux notions déjà étudiées. La première précise l'effet de perturbations sur la matrice par la définition de $\varepsilon$-spectre ou pseudospectre. La deuxième consiste à dénombrer les valeurs propres entourées par une courbe a priori donnée dans le plan complexe. A partir de travaux antérieurs, on combine ici les deux approches avec l'objectif de mettre en commun les factorisations LU de la résolvante nécessaire aux deux approches et d'en diminuer le nombre. Les codes obtenus sont parallélisés.

**ABSTRACT.** This article deals with the localization of eigenvalues of a large sparse and not neces-
sarily symmetric matrix in a domain of the complex plane. It combines two studies carried out earlier.
The first work deals with the effect of applying small perturbations on a matrix, and referred to as
$\varepsilon$-spectrum or pseudospectrum. The second study describes a procedure for counting the number
of eigenvalues of a matrix in a region of the complex plain surrounded by a closed curve. The two
methods are combined in order to share the LU factorization of the resolvent, that intervenes in the
two methods, so as to reduce the cost. The codes obtained are parallelized.

**MOTS-CLÉS** : Pseudospectre, dénombrement de valeurs propres, déterminant, valeur singulière minimum, orbite, trace, client-serveur, Maître, Travailleurs, Accélération, Efficacité

**KEYWORDS** : Pseudospectrum,eigenvalue counting, determinant, minimum singular value, orbit, trace, Client-server,Master, Workers, speedup, efficiency
# 1. Sensibilité aux perturbations des valeurs propres d'une matrice

La localisation de valeurs propres pour une grande matrice creuse dans un domaine du plan complexe est une tâche délicate, spécialement lorsque la matrice n'est pas symétrique. En effet, pour une matrice symétrique A et pour toute matrice symétrique de perturbation $\Delta$, on est assuré que pour toute valeur propre $\lambda$ de A, il existe une valeur propre $\tilde{\lambda}$ de $A + \Delta$ telle que

$$|\lambda - \tilde{\lambda}| \le ||\Delta||.$$

où la norme matricielle considérée est la norme $L_2$, comme elle l'est dans tout l'article. Par contre si la matrice A n'est plus supposée symétrique, la variation des valeurs propres dépend de la proximité d'une matrice non diagonalisable dans le voisinage de A puisque la perturbation d'un bloc de Jordan de dimension $k$ éclate les valeurs propres en fonction de la racine $k$-ième de la perturbation.

En effet si $A = XJX^{-1}$ est la forme de Jordan et $m$ la dimension du plus grand bloc de Jordan pour $\lambda \in \Lambda(A)$ spectre de A, alors pour $\Delta$ il existe $\tilde{\lambda} \in \Lambda(A + \Delta)$ tel que [10] pp.174 :

$$\frac{|\lambda - \tilde{\lambda}|}{(1 + |\lambda - \tilde{\lambda}| + \dots + |\lambda - \tilde{\lambda}|^{m-1})^{1/m}} \le (||X^{-1}|| ||X|| ||\Delta||)^{1/m}.$$

Afin de prendre en compte l'effet de perturbations sur la matrice, la notion de $\varepsilon$-spectre ou pseudo-spectre d'une matrice $A \in \mathbb{C}^{n \times n}$ a été définie indépendamment mais simultanément par Godunov [3] et Trefethen [11].

Soit $A \in \mathbb{C}^{n \times n}$, et $\varepsilon \ge 0$. L' $\varepsilon$-spectre d'une matrice A est défini par :

$$
\begin{align*}
\Lambda_{\varepsilon}(A) &= \{z \in \mathbb{C} : \exists \Delta \in \mathbb{C}^{n \times n} \text{ avec } ||\Delta|| \le \varepsilon \text{ et z est une valeur propre de } A + \Delta\} \text{[1]} \\
&= \{z \in \mathbb{C} : ||(A - zI)^{-1}|| \ge \frac{1}{\varepsilon}\} \tag*{[2]} \\
&= \{z \in \mathbb{C} : \sigma_{\min}(A - zI) \le \varepsilon\}. \tag*{[3]}
\end{align*}
$$

où $\sigma_{\min}(A - zI)$ est la plus petite valeur singulière de $(A - zI)$. Construire un $\varepsilon$-spectre d'une matrice A consiste donc à déterminer une courbe de niveau $(\Gamma)$ de la fonction

$$z \in \Omega \subset \mathbb{C} \rightarrow \sigma_{\min}(A - zI),$$

où $\Omega$ est la région en considération dans $\mathbb{C}$.

On peut aussi développer une approche alternative : étant donnée une courbe $(\Gamma)$ dans le plan complexe, il s'agit de compter le nombre de valeurs propres de A entourées par $(\Gamma)$. Cela revient à évaluer l'intégrale

$$\frac{1}{2i\pi} \int_{\Gamma} \frac{d}{dz} \log \det(zI - A) dz.$$

Ce problème a été étudié dans [1], [2], [9] où plusieurs méthodes d'intégration étaient proposées. Plus récemment, dans [5] le contrôle du pas d'intégration a fait l'objet d'une étude fine. L'objectif est maintenant de combiner les deux approches pour la localisation et le dénombrement de valeurs propres d'une matrice A.

Il existe deux types d'approche pour la construction d'un pseudo-spectre : la méthode de maillage uniforme et les méthodes de suivi de contour. Dans cet article nous utiliserons une méthode de suivi de contour appelée méthode PAT.
## 2. La méthode de suivi de contour à l'aide de triangles : Méthode PAT

Cette méthode est présentée dans [6, 7]. On construit une séquence de triangles équilatéraux adjacents qui recouvrent le contour. On détermine un point du contour par bisssection de chaque coté du triangle qui intersecte le contour (voir la figure 1).

<figure><img src="image_1.png" /><figcaption>Figure 1. PAT : Exemple de suivi de contour par des triangles équilatéraux.</figcaption></figure>

### 2.1. Définition du maillage

La définition d'un maillage adapté à la ligne de niveau repose sur la recherche de deux points; l'un $z_i$, intérieur (bord compris) et l'autre $z_e$, extérieur au pseudo spectre $\Lambda_e(A) = \{z \in \mathbb{C} | \sigma_{\min}(zI - A) \le \varepsilon\}$ tels que $|z_e - z_i| = \tau$, où $\tau$ désigne la résolution du maillage.

Ces deux points sont dits $\varepsilon$-séparés et permettent de générer un maillage triangulaire $S(z_i, z_e)$ où $T_0 = \{S_{0,0}, S_{1,0}, S_{0,1}\}$ constitue le premier triangle. Un moyen pour trouver $z_i$ et $z_e$ consiste à partir d'une estimation $z_0$ d'une valeur propre proche d'un point de référence $z_{ref}$ et à construire la suite

$$z_{k+1} = z_0 + 2^k \tau e^{i\theta} \text{ pour un angle } \theta \text{ donné.}$$

Étant donnés un point intérieur $z_i$ et un point extérieur $z_e$, le réseau uniforme $\mathcal{S}$ de triangles équilatéraux

$$S(z_i, z_e) = \{S_{kl} = z_i + k (z_e - z_i) + l (z_e - z_i)e^{i\frac{\pi}{3}}, (k,l) \in \mathbb{Z}^2\}$$

définit un maillage uniforme de points vérifiant : $|S_{k,l+1} - S_{k,l}| = |S_{k+1,l} - S_{k,l}| = \tau$, où $\tau = |z_i - z_e|$, pour tout $(k,l) \in \mathbb{Z}^2$. Cette grille permet de définir deux classes de triangles $T_{kl} = \{S_{k,l}, S_{k+1,l}, S_{k,l+1}\}$ et $\tilde{T}_{kl} = \{S_{k,l}, S_{k+1,l}, S_{k+1,l-1}\}$ qui sont liées par les rotations $R(S_{k,l}, \frac{\pi}{3})$ et $R'(S_{k,l}, -\frac{\pi}{3})$ (voir la figure 2).
<figure><img src="image_2.png" /><figcaption>**Figure 2.** PAT : Le réseau des triangles équilatéraux.</figcaption></figure>

## 2.2. Construction d'une orbite

Soit $\Omega(z_i, z_e)$ l'ensemble des triangles équilatéraux défini à partir d'un point intérieur $z_i$ et d'un point extérieur $z_e$ ; soit $\mathcal{O}_L$ la partie de $\Omega(z_i, z_e)$ constituée des triangles qui ont au moins deux sommets $\varepsilon$-séparés ; $\mathcal{O}_L$ est un ensemble fini. Pour chaque triangle $T \in \mathcal{O}_L$, il n'y a qu'un sommet qui soit seul de sa catégorie(point intérieur ou point extérieur). Ce point est appelé le pivot de T et est noté $p(T)$. On définit alors la bijection $F$ de l'ensemble $\mathcal{O}_L$ dans lui-même par

$$ 
\left\{ 
\begin{array}{l}
F(T) = \mathcal{R}(p(T), \theta) T, \\
\text{où } \theta = \frac{\pi}{3} \text{ si } p(T) \text{ est intérieur, ou } \theta = -\frac{\pi}{3}, \text{ sinon} \\
\text{et où } \mathcal{R}(z, \theta) \text{ est la rotation de centre } z \in \mathbb{C} \text{ et d'angle } \theta.
\end{array}
\right. 
$$

La transformation $F$ est illustrée par la figure 3.

En partant d'un triangle quelconque $T_0 \in \mathcal{O}_L$, on peut construire la $F$-orbite de $T_0$ :

$$ O(T_0) = \{T_k \equiv F^k(T_0), k \in \mathbb{Z}\} \subset \mathcal{O}_L. $$

Puisque tous les triangles construits intersectent la ligne de niveau, l'ensemble est constitué dans une partie bornée du plan complexe. Il s'en suit qu'il existe un entier $n$ tel que $F^n(T) = T$ et par conséquent la suite $(T_k)_{k \in \mathbb{Z}}$, est périodique ; comme nous utilisons des coordonnées entières pour identifier les points $S_{kl}$, le test $F^n(T) = T$ pourra être implémenté de manière sûre. Comme la somme des angles des rotations construisant une orbite doit être nulle à $2\pi$ près, le nombre de triangles distincts d'une orbite est toujours pair.

## 2.3. L'algorithme PAT

L'algorithme est décrit dans l'algorithme 1 et la méthode PAT est définie en trois étapes :

**ETAPE I : construction du triangle initial.** Choisir $z_0$ point intérieur. Cela est réalisé en calculant une valeur propre de la matrice $A$ dans le voisinage d'une valeur de référence $z_{ref}$ fournie par l'utilisateur. A partir des paramètres $\tau > 0$ et $\theta$, on construit
<figure><img src="image_3.png" /><figcaption>**Figure 3.** PAT : La transformation F.</figcaption></figure>

**Algorithm 1** PAT : suivi d'une ligne de niveau $\partial\Lambda_{\varepsilon}(A)$ à l'aide de triangles.

**Require:** $A \in \mathbb{C}^{m \times n}, \varepsilon > 0, \tau > 0$ la taille des mailles du réseau, $\theta$ l'inclinaison du réseau, $z_0$ un point intérieur ($\sigma_{\min}(A - z_0I) \le \varepsilon$).

**Ensure:** Un triangle initial $T_0 = (z_i, z_e, z') \in \mathcal{O}_L$, son orbite $O(T_0)$, une ligne polygonale de $N$ sommets ($N = \text{card}(O(T_0))$) inscrite dans le bord $\partial\Lambda_{\varepsilon}(A)$ du pseudo-spectre $\Lambda_{\varepsilon}(A)$.

1: $h = \tau e^{i\theta}$; $Z = z_0 + h$; $I_t = 0$; {Recherche d'un point extérieur}

2: **while** $\sigma_{\min}(A - z_0I) \le \varepsilon$,

3: $I_t = I_t + 1$; $h = 2h$; $Z = z_0 + h$;

4: **end while**

5: $z_i = z_0$; $z_e = Z$; {ETAPE I : Construction du triangle initial}

6: **do** $k = 1 : I_t$,

7: $z = (z_e + z_i)/2$;

8: **if** $\sigma_{\min}(A - zI) \le \varepsilon$, **then**

9: $z_i = z$;

10: **else**

11: $z_e = z$;

12: **end if**

13: **end**

14: $z' = z_i + \tau e^{i(\theta+\pi/3)}$; $T_0 = (z_i, z_e, z')$;

15: $T_1 = F(T_0)$; $i = 1$; $\mathcal{I} = \{[z_i, z_e], [z', p(T_1)]\}$; {ETAPE II : Construction d'une F-orbite}

16: **while** $T_i \ne T_0$,

17: $T_{i+1} = F(T_i)$; $i = i + 1$;

18: Choisir le nouveau côté de $T_i$ qui intersecte $\partial\Lambda_{\varepsilon}(A)$ et l'insérer dans la liste $\mathcal{I}$;

19: **end while**
{ETAPE III : Construction de la ligne polygonale inscrite dans $\partial\Lambda_{\varepsilon}(A)$}

20: **doall** $I \in \mathcal{I}$,

21: Extraire $z_i \in \partial\Lambda_{\varepsilon}(A) \cap I$ par bissections itérées;

22: **end**
alors un point Z du réseau, extérieur à $\Lambda_{\varepsilon}(A)$ sur la droite $z = z_0 + k\tau e^{i\theta}$, point qui existe car $\lim_{|z|\to\infty} \sigma_{\min}(zI - A) = \infty$. Pour cela, on recherche pour $k = 1, 2, \dots$, le premier point $z_k = z_0 + 2^k \tau e^{i\theta}$ qui soit extérieur. Puisque ce point est sommet dans le réseau $\mathcal{S}$, par bissection successive de l'intervalle $[z_{k-1}, z_k]$, on en déduit, par ajustement, un intervalle de longueur de côté $\tau$ qui intersecte le bord $\partial\Lambda_{\varepsilon}(A)$ de $\Lambda_{\varepsilon}(A)$ et on construit un triangle équilatéral sur ce segment.

**ETAPE II : construction de l'orbite.** La fonction $F$ est successivement appliquée jusqu'à la fermeture de l'orbite. Le test d'égalité de triangles est sûr puisqu'il est fondé sur un test d'entiers.

**ETAPE III : construction de points de la ligne de niveau.** Extraction par bissection. Cette étape peut être omise si on ne recherche qu'une ligne polygonale contenant le pseudo-spectre (voir section 4).

Le nombre minimum de triangles dans une orbite est de 6. Ceci correspond à la situation où le pivot est le même pour tous les triangles de l'orbite et par conséquent les angles de toutes les rotations sont égaux à $\pi/3$. Une telle orbite correspond souvent à une valeur de $\tau$ qui est grande par rapport à la longueur de $\partial\Lambda_{\varepsilon}(A)$.

**Proposition 1** *Supposons que l'algorithme 1 détermine une orbite $O(T_0)$ qui contienne $N$ triangles et que $N > 6$. Soit $\ell$ la longueur de la courbe $\partial\Lambda_{\varepsilon}(A)$ et $\eta$ la précision pour la détermination des points du contour. La complexité du calcul est définie par le nombre $\mathcal{N}_{\varepsilon}$ de points où $s(z) = \sigma_{\min}(A - zI)$ est calculé. Ce nombre est tel que*

$$ \mathcal{N}_{\varepsilon} \le (N+1)[\log \frac{\tau}{\eta}] \quad \text{et} \quad N = O\left(\frac{\ell}{\tau}\right). $$

Pour calculer la valeur singulière minimale de $B = zI - A$, on choisit la méthode de Lanczos pour calculer la valeur propre maximale de la matrice augmentée

$$ \begin{pmatrix} 0 & B^{-1} \\ B^{-*} & 0 \end{pmatrix} $$

grâce à une factorisation LU de la matrice $B$.

Quand la matrice $A$ et le point $z_0$ sont réels avec une inclinaison $\theta$ nulle du réseau, alors l'orbite admet comme axe de symétrie l'axe des réels. Dans ce cas, on peut arrêter la détermination de l'orbite dès que le nouvel intervalle intégrant la liste $\mathcal{I}$ est réel. Ce cas réduit de moitié le travail nécessaire à la construction de l'orbite et de la ligne polygonale recherchée.

### 3. La méthode EIGENCNT

Soit une matrice $A \in \mathbb{C}^{n \times n}$ et une ligne polygonale $\Gamma$ ne contenant pas de valeur propre de $A$. Alors le nombre de valeurs propres de $A$ qui appartiennent au domaine de frontière $\Gamma$ est donné par la formule de Cauchy

$$ N_{\Gamma} = \frac{1}{2i\pi} \int_{\Gamma} \frac{f'(z)}{f(z)} dz $$

où $f(z) = \det(zI - A)$ est le polynôme caractéristique de $A$. Pour calculer effectivement cette intégrale par morceaux, on considère une discrétisation $\Gamma = \bigcup_{i=0}^{N-1} [z_i, z_{i+1}]$ de la ligne
polygonale où $[z_i, z_{i+1}]$ dénote un segment de la ligne de niveau avec pour extrémité $z_i$ et $z_{i+1}$. Soit $z$ et $z+h$ deux points de $(\Gamma)$. Puisque

$$ (z+h)I - A = (zI - A) + hI \\ = (zI - A)(I + hR(z)), $$

où $R(z) = (zI - A)^{-1}$, on peut écrire

$$ f(z+h) = f(z) \det(I + hR(z)). \quad [4] $$

En posant $\Phi_z(h) = \det(I + hR(z))$, on a

$$ \int_z^{z+h} \frac{f'(z)}{f(z)} dz = \log(f(z+h)) - \log(f(z)) \\ = \log\left(\frac{f(z+h)}{f(z)}\right) \\ = \log(\Phi_z(h)) \\ = \log|\Phi_z(h)| + i\arg(\Phi_z(h)). $$

Le déterminant est obtenu à partir de la factorisation LU de la matrice $(zI - A)$ et seule la partie imaginaire est utile. Une méthode robuste et parallèle de calcul de déterminant de grandes matrices creuses est proposée dans [4].

## 3.1. Contrôle du pas d'intégration

Afin d'intégrer de manière sûre, (c'est-à-dire de ne pas perdre de valeurs propres pendant l'intégration), il est nécessaire d'introduire un contrôle de pas. Dans [5], on montre que la condition

$$ |\operatorname{Arg}(\Phi_z(s))| < \pi, \quad \forall s \in [0, h], \quad [5] $$

où Arg représente l'argument principal, est appelée **Condition (A)** et qu'elle est nécessaire et suffisante ; elle est équivalente à :

$$ \Phi_z(s) \notin (-\infty, 0], \quad \forall s \in [0, h]. $$

Une condition plus restrictive que la **Condition (A)** mais suffisante est définie par :

$$ |\Phi_z(s) - 1| < 1, \quad \forall s \in [0, h]; \quad [6] $$

elle sera appelée **Condition (B)**. Si elle n'est satisfaite qu'en $z+h$, elle se réduit à :

$$ |\Phi_z(h) - 1| < 1, \quad [7] $$

et sera appelée **Condition (B')** ; cette condition est celle utilisée dans [1]. S'il est clair que **Condition (B)** implique **Condition (A)** il n'en est pas le cas pour **Condition (B')**. Comme il serait difficile de vérifier la **Condition (B)**, on utilise une approximation linéaire donnée par :

$$ \Phi_z(u) \approx 1 + u \Phi'_z(0) $$

d'où la **Condition (C)** suivante :

$$ |\delta| |\Phi'_z(0)| < 1. \quad [8] $$
où

$$ \Phi'_{z}(0) = \left[ \frac{d}{du} (\det(I + u R(z))) \right]_{u=0}, $$

avec $R(z) = (A - zI)^{-1}$. Lorsque $\Phi'_{z}(0) \approx 0$, l'approximation linéaire n'est pas valable et

$$ \Phi_{z}(\delta) \not\approx 1 + \delta \Phi'_{z}(0). $$

Pour garantir un bon pas d'intégration, on impose que la **Condition (C)** soit satisfaite aux deux extrémités de l'intervalle :

$$ |\delta| \le \min (|\Phi'_{z}(0)|, |\Phi'_{z+\delta}(0)|). $$

Les trois conditions sont logiquement reliées comme suit :

(A) Condition nécessaire et suffisante

↑

(B) Condition suffisante

≈

(C) Approximation de la condition suffisante

Pour que la **Condition C** soit satisfaite à chaque morceau, il sera nécessaire d'ajouter des points intermédiaires supplémentaires. Plus précisément : si **Condition (C)** est violée, on insère M points uniformément répartis entre z and z + h où

$$ M = \min (\lceil|h| |\Phi'_{z}(0)| \rceil, M_{\max}), \quad [9] $$

et $M_{\max}$ étant un paramètre pré-défini. Lorsque la **Condition (B')** est violée et la **Condition (C)** satisfaite, on insère le point $z + h/2$.

## 3.2. Calcul de la dérivée

Pour le calcul de dérivée, à partir d'un simple calcul polynomial, on peut montrer que :

$$ \Phi'_{z}(0) = \left[ \frac{d}{du} (\det(I + u R(z))) \right]_{u=0} \quad [10] $$

$$ = \operatorname{trace}(R(z)). \quad [11] $$

La factorisation $P(A - zI)Q = LU$ est déjà faite pour le calcul du $\det(I + \delta R(z))$. Cependant le coût du calcul de la trace($R(z)$) est très élevé puisqu'il nécessite de résoudre n systèmes linéaires ; en effet, puisque $\operatorname{trace}(R(z)) = \sum_{i=1}^{n} R(z)_{ii}$ où

$$ R(z)_{ii} = e_i^T (zI - A)^{-1} e_i, \quad [12] $$

il faut résoudre les systèmes linéaires suivants :

$$ (zI - A)x_i = e_i, \text{ pour } i = 1, \dots, n. \quad [13] $$

Nous proposons deux méthodes pour estimer la dérivée à moindre coût : par taux d'accroissement et par estimation stochastique.
### 3.2.1. Estimation de la dérivée par taux d'accroissement

Dans [5] on envisage d'approcher la dérivée par taux d'accroissement. Nous choisissons ici

$$ \left[ \frac{d}{d\delta} (\det(I + \delta R(z))) \right]_{\delta=0} \approx \frac{\Phi_z(s) - \Phi_z(0)}{s}, $$

où $s = O(\sqrt{\epsilon}|z|)$; on exploitera la structure creuse de la matrice lors de sa factorisation LU. Étant donnés $z$ et $z+h$, la dérivée de $\Phi_z$ au point 0 est estimée par :

$$ \Phi'_{z}(0) \approx \frac{\Phi_{z}(s) - 1}{s}, $$

où $s = \alpha h$ avec $\alpha = \min(10^{-6}\mu/|h|, 1)$, and $\mu = \max_{z \in \Gamma} |z|$. Ainsi, une deuxième factorisation LU est nécessaire pour le calcul de $\Phi_z(s)$.

### 3.2.2. Estimation de la dérivée par méthode stochastique

Cette approche consiste à trouver un estimateur de la trace en diminuant le nombre de systèmes linéaires à résoudre. On choisit de façon aléatoire, $N$ indices entre 1 et $n$, la moyenne arithmétique de la trace obtenue avec les $N$ indices est supposée approcher $\frac{\text{trace}(R(z))}{n}$, ce qui nous amène à écrire :

$$ \sum_{i=1}^{n} R(z)_{ii} \approx \frac{n}{N} \sum_{i=1..N; j_i \in 1..N} R(z)_{ji j_i}. \quad [14] $$

#### Dimensionnement de l'échantillon

Pour estimer le nombre de systèmes linéaires à résoudre afin d'avoir une approximation de la trace de la résolvante avec au moins un chiffre significatif, nous choisissons un point $z_k$ de la ligne de niveau et nous approchons la trace pour différentes valeurs de $N$. Pour chaque valeur de $N$, nous effectuons 10 tirages et nous calculons la moyenne des erreurs, l'erreur maximale est également mentionnée pour illustrer le pire des cas.

La figure 4 présente les résultats obtenus pour la matrice ADD32 de taille $n = 4960$ et la matrice CRY10000 de taille $n = 10000$. Sur cette figure, nous pouvons remarquer qu'avec $N = 400$ nous avons une bonne approximation de la trace pour la matrice ADD32 de taille $n = 4960$ alors que pour la matrice CRY10000 de taille $n = 10000$, il faut résoudre $N = 1000$ pour être sûr de toujours avoir une bonne approximation de la trace de la résolvante.

### 3.3. L'algorithme EIGENCNT

La méthode EIGENCNT est décrite dans l'algorithme 2. A partir d'une liste de points Z, l'algorithme inclut de nouveaux points selon que la **Condition (C)** ou la **Condition (B')** est satisfaite ou pas. La complexité de l'algorithme dépend du nombre de déterminants à calculer et du nombre de systèmes linéaires à résoudre. Pour chaque point $z \in Z$, on calcule $\det(zI - A)$ et $\Phi'_z(0)$ qui nécessite soit un nouveau calcul de déterminant lorsque l'on estime la dérivée par le taux d'accroissement soit la résolution de $N$ systèmes linaires avec $N \ll n$, lorsque l'on utilise la méthode stochastique.

Lorsque l'estimation de la dérivée se fait avec le taux d'accroissement, si en sortie, la liste Z contient $N_z$ points, alors la complexité de l'algorithme peut s'exprimer sous la forme :

$$ \mathcal{C}_{\text{eig}} = 2\mathcal{C}_{\text{LU}}N_z, $$

où $\mathcal{C}_{\text{LU}}$ est le nombre d'opérations pour la factorisation LU complexe de $(zI - A)$.

Lorsque l'estimation de la dérivée est stochastique, pour chaque point $z \in Z$, on calcule
<figure><img src="image_4.png" /><figcaption>Figure 4. Estimation du nombre de systèmes nécessaires pour approximer la trace.</figcaption></figure>

$det(zI - A)$, avec la même factorisation LU nous résolvons N systèmes linéaires pour estimer $\Phi'_z(0)$ ; avec $N_z$ points en sortie, la complexité de l’algorithme peut alors s’exprimer sous la forme :

$$ \mathcal{C}_{eig} = \mathcal{C}_{LU}N_z + \mathcal{C}_{solve}N, $$

où $\mathcal{C}_{solve}$ est le nombre d’opérations pour la résolution du système $(zI - A)x_i = e_i$ (c’est à dire une descente-remontée des systèmes triangulaires), pour un indice $i$ donné.

# 4. La combinaison des deux méthodes et leur parallélisation

## 4.1. Combinaison des méthodes PAT et EIGENCNT

On combine maintenant les deux procédures PAT et EIGENCNT : deux approches sont possibles.

### 4.1.1. Première approche

L’algorithme PAT détermine une liste de points Z de la courbe de niveau ($\Gamma_\epsilon$). A chaque point $z \in Z$, $det(A - zI) = sign(P)sign(Q)det(U), \sigma(z) = \sigma_{min}(A - zI)$ sont calculés, où $P(A - zI)Q = LU$ est la factorisation LU. Le calcule de $\sigma_{min}(A - zI)$ se fait par Lanczos sur

$$ M_z = \begin{pmatrix} 0 & U^{-1}L^{-1} \\ L^{-*}U^{-*} & 0 \end{pmatrix}, $$

puisque $\sigma_{min}(A - zI) = \sigma_{min}(LU)$. On fait appel à EIGENCNT en lui envoyant la liste des points.

### 4.1.2. Seconde approche

L’algorithme PAT construit plutôt une liste de triangles qui recouvrent la courbe $\Gamma$, sans interpolation pour obtenir des points la courbe $\Gamma$. On peut se servir de l’algorithme 3. A partir des noeuds des triangles de la liste, on déduit deux listes de points : $Z_{int}$ de points intérieurs et $Z_{ext}$ de points extérieurs. On fait appel à EIGENCNT en lui envoyant la liste $Z_{ext}$. Nous utilisons ici la seconde approche.

### 4.1.3. Méthode combinée

Les grandes étapes du calcul sont :
**Algorithm 2 EIGENCNT**: Algorithme de dénombrement de valeurs propres dans un domaine du plan entouré par $(\Gamma)$.

**Require:** ENTREES : $Z=\{\text{noeuds de } (\Gamma)\}$, $M_{\text{pts}}$ = nombre maximum de points autorisés, $M_{\text{max}}$ = nombre maximum de points à insérer simultanément, SORTIES : $\text{neg}$ = nombre de valeurs propres entourées par $(\Gamma)$;

1: Statut($Z$)=-1;

2: **while** Statut($Z$)≠ 0 et longueur($Z$) < $M_{\text{pts}}$,

3: **do** $z \in Z$ tel que Statut($z$)===-1,

4: calcul $det(zI - A)$ et $\Phi'_z(0)$;

5: Statut($z$) = 1;

6: **end**

7: **do** $z \in Z$ tel que Statut($z$)=1,

8: **if** Condition (C) non satisfaite au point $z$, **then**

9: Générer $M$ points $\tilde{Z}$ comme dans (9);

10: $Z=Z \cup \tilde{Z}$; Statut($\tilde{Z}$)=-1;

11: **else if** Condition (B') non satisfaite au point $z+h$; **then**

12: $Z=Z \cup \{z+h/2\}$; Statut($z+h/2$)=-1;

13: **else**

14: Statut($z$)=0;

15: **end if**

16: **end**

17: **if** pas de nouveaux points à insérer dans $Z$; **then**

18: **do** $z \in Z$,

19: **if** Condition (C) est non satisfaite anormalement; **then**

20: $Z=Z \cup \{z-h/2\}$; Statut($z-h/2$)=-1;

21: **end if**

22: **end**

23: **end if**

24: **end while**

25: Integral = $\sum_{z \in Z} \text{Arg}(\Phi_z(0))$; neg = round(Integral/2$\pi$);

**Algorithm 3 TRIANGLESET**: Algorithme de construction de liste de triangles recouvrant $\Gamma$.

**Require:** ENTREES : $\varepsilon, \tau$ et $\tilde{z}_0$ tel que $s(\tilde{z}_0) \le \varepsilon$

**Ensure:** construire $T_0 \in T_L$ triangle initial;

1: **do** $k \ge 0$,

2: $T_{k+1} = F(T_k)$;

3: **if** $T_{k+1} = T_0$, **then**

4: $L = k$;

5: **break**;

6: **end if**

7: **end**
– Choisir la résolution $\tau > 0$ et l'inclinaison $\theta$ du maillage triangulaire et définir un point de référence $z_{ref} \in \mathbb{C}$ autour duquel on recherche des valeurs propres. Calculer une valeur propre $z_0$ voisine de $z_{ref}$ par itération inverse ou autre méthode plus rapide.

– Par PAT, construire une liste de triangles à partir de $z_0$ et $\tau e^{i\theta}$ (étapes I et II de PAT). En déduire la ligne polygonale $Z_{ext}$ des sommets extérieurs de la liste des triangles de l'orbite. Pour tout sommet $z$, les nombres $\sigma_{min}(zI - A)$ et $\det(zI - A)$ sont donc calculés. Ces deux quantités complexes sont obtenues à partir de la même factorisation LU de la matrice $(zI - A)$.

– A l'aide de EIGENCNT, calculer le nombre de valeurs propres entourées par la ligne polygonale trouvée précédemment. Cette étape rajoutera des points sur les segments de la ligne polygonale si la longueur est trop grande pour une intégration sûre. Chaque nouveau point entraîne une ou deux factorisations LU de la résolvante.

### 4.1.4. Résultats Numériques

Les tests sont conduits sur des matrices issues soit de Matrix Market [14] soit de la collection de Tim Davis[15]. Les résultats sont rassemblés dans le tableau 1. Pour chaque exécution en plus de la matrice, on doit spécifier un point de référence $z_{ref}$, une taille du maillage $\tau$, un niveau de la ligne de niveau à détecter $\eta$ et le nombre $N_s$ de systèmes linéaires à résoudre pour calculer la dérivée $\Phi'_z(0)$ en un point $z$. Le programme nous retourne le nombre $nvp$ de valeurs propres entourées par la ligne polygonale formée par les points extérieurs de chaque triangle de la chaîne des triangles construite, le nombre $nT$ de triangles constituant la chaîne et le nombre $nintv$ d'intervalles obtenu à la fin de l'intégration. Les tests numériques ont été menés sur une machine équipée de deux processeurs,

<table><caption>**Tableau 1.** Nombre de valeurs propres localisées pour différentes matrices en différentes régions</caption><thead><tr><th>Name</th><th>n</th><th>Zref</th><th>τ</th><th>η</th><th>nvp</th><th>Ns</th><th>nT</th><th>nintv</th></tr></thead><tbody><tr><td>GRCAR100</td><td>100</td><td>1.7+1.1i</td><td>10<sup>-1</sup></td><td>10<sup>-6</sup></td><td>100</td><td>10</td><td>344</td><td>2109</td></tr><tr><td>FS-680-1</td><td>680</td><td>7x10<sup>9</sup></td><td>10<sup>5</sup></td><td>10<sup>5</sup></td><td>6</td><td>10</td><td>646</td><td>2763</td></tr><tr><td>ADD20</td><td>2395</td><td>0.1+0.1i</td><td>10<sup>-5</sup></td><td>10<sup>-4</sup></td><td>43</td><td>100</td><td>496</td><td>251</td></tr><tr><td>CRY2500</td><td>2500</td><td>1.5+1.7i</td><td>10<sup>-1</sup></td><td>10<sup>-3</sup></td><td>569</td><td>100</td><td>170</td><td>13199</td></tr><tr><td>UTM3060</td><td>3060</td><td>-1.5+1.2i</td><td>10<sup>-1</sup></td><td>10<sup>-3</sup></td><td>4</td><td>100</td><td>6</td><td>3083</td></tr><tr><td>MHD4800B</td><td>4800</td><td>i</td><td>10<sup>-4</sup></td><td>10<sup>-6</sup></td><td>1230</td><td>100</td><td>1370</td><td>29630</td></tr><tr><td>ADD32</td><td>4960</td><td>0.1+0.1i</td><td>10<sup>-5</sup></td><td>10<sup>-4</sup></td><td>89</td><td>100</td><td>408</td><td>1240</td></tr><tr><td>UTM5940</td><td>5940</td><td>1.5+1.2i</td><td>10<sup>-2</sup></td><td>10<sup>-6</sup></td><td>1</td><td>100</td><td>6</td><td>926</td></tr><tr><td>CRY10000</td><td>10000</td><td>1.5+1.7i</td><td>10<sup>-1</sup></td><td>10<sup>-3</sup></td><td>652</td><td>1000</td><td>312</td><td>31552</td></tr><tr><td>BCSSTK17</td><td>10974</td><td>1+i</td><td>10<sup>-1</sup></td><td>10<sup>-6</sup></td><td>518</td><td>1000</td><td>6</td><td>5656</td></tr><tr><td>BCSSTM25</td><td>15439</td><td>0.1+0.1i</td><td>10<sup>-3</sup></td><td>10<sup>-2</sup></td><td>6</td><td>1000</td><td>150</td><td>706</td></tr><tr><td>BCSSTM25</td><td>15439</td><td>0.1+0.1i</td><td>10<sup>-3</sup></td><td>10<sup>-1</sup></td><td>30</td><td>1000</td><td>1950</td><td>3477</td></tr><tr><td>E40R5000</td><td>17281</td><td>1+2i</td><td>10<sup>-1</sup></td><td>10<sup>-6</sup></td><td>2</td><td>1000</td><td>6</td><td>2090</td></tr><tr><td>AF23560</td><td>23560</td><td>i</td><td>10<sup>-1</sup></td><td>10<sup>-6</sup></td><td>1</td><td>1000</td><td>6</td><td>731</td></tr><tr><td>BCSSTK31</td><td>35588</td><td>0.1+0.1i</td><td>10<sup>-1</sup></td><td>10<sup>-6</sup></td><td>2318</td><td>1000</td><td>6</td><td>15606</td></tr><tr><td>FINAN512</td><td>74752</td><td>2+25i</td><td>10<sup>-1</sup></td><td>10<sup>-3</sup></td><td>611</td><td>1000</td><td>238</td><td>26884</td></tr></tbody></table>

chacun ayant 6 cœurs et 24 GB de mémoire; chaque cœur a un cache L1 de 32 KB et un cache L2 de 256 KB; les 6 cœurs d'un processeur partagent un cache L3 de 12 MB; chaque cœur est un Intel(R) Xeon(R) de 3.47 GHZ de vitesse d'horloge; la machine est équipée d'un système d'exploitation Red had 4.6.2-1; du compilateur gcc 4.6.2; le code est écrit en C et compilé avec pour flag -g -w -O3; le solveur UMFPACK est utilisé pour la factorisation LU de matrice creuses et la résolution des systèmes linéaires.
### 4.2.2. Le Maître

Pour maximiser l'efficacité, le processus *Maître* doit générer le maximum possible de tâches pour occuper les *Travailleurs*.

Au démarrage, le *Maître* génère six tâches `compute_chain()`. Ces tâches sont envoyées aux *Travailleurs* et le *Maître* passe dans un état d'attente de messages.

A la réception d'un message, le *Maître* analyse les données reçues, si elles proviennent d'une tâche `compute_chain()`, le *Maître* teste les conditions, génère de nouveaux points $z_k$ puis pour chaque point $z_k$ généré, crée une tâche `compute Det$(z_k)$`. Ces tâches sont rangées dans une file d'attente puis servies aux *Travailleurs* au fur et à mesure qu'ils se libèrent. Si les données proviennent d'une tâche `compute Det$(z_k)$`, le *Maître* les rangent et attribue une autre tâche au *Travailleur*.

### 4.2.3. Le Travailleur

Le *Travailleur* est un processus simple qui répond aux requêtes du processus *Maître*. Il offre aux processus *Maître* les services suivants :

`compute\_chain()` : étant donné $\tau$ une taille du maillage, $\eta$ un niveau de la ligne à détecter et $d$, une direction, construire l'orbite dans la direction $d$ en retournant chaque point extérieur au *Maître*. Cette tâche nécessite pour chaque sommet du triangle, le calcul d'un $\sigma_{min}$. Lorsqu'il s'agit d'un point extérieur $z_k$ (dont à retourner au *Maître*), il faut en plus du calcul de $\sigma_{min}$, calculer $\Phi'_{z_k}(0)$ en résolvant N systèmes linéaires lorsque l'on utilise l'estimation stochastique pour la dérivée ; cette résolution utilise la même factorisation LU effectué pour le calcul de $\sigma_{min}$.

`compute\_det(z_k)` : étant donné un point $z_k$, calculer $det(z_kI - A)$ et $\Phi'_{z_k}(0)$. Cette tâche nécessite une factorisation LU pour le calcul du déterminant de $(z_kI - A)$. Lorsque l'estimation de la dérivée est stochastique, il faut résoudre N systèmes linéaires pour le calcul de $\Phi'_{z_k}(0)$ alors que lorsque l'estimation de la dérivée utilise le taux d'accroissement, il faut effectuer une nouvelle factorisation LU en un point voisin de $z_k$.

Chacun de ces services nécessite un accès à la matrice. Pour cela cette matrice est diffusée aux différents *Travailleurs* au démarrage de l'application. Dans le cas des matrices de très grande taille, un processeur physique risque de manquer de ressources mémoire pour assurer le service. En plus la complexité de ces services augmente considérablement car le calcul de valeurs singulières, de déterminants et de dérivées devient très long. Un second niveau de parallélisme peut être introduit si les nœuds de calcul utilisés ont plusieurs cœurs, ceci en exploitant soit la méthode décrite dans [4] lorsque la dérivée est approchée par taux d'accroissement soit la résolution parallèle des systèmes linéaires.

## 4.3. Résultats Numériques

Nous avons implémenté notre algorithme sur la grille de calcul igrida du centre de recherche INRIA de Rennes et particulièrement sur le cluster lambda qui est équipé de 11 nœuds de calculs. Chaque nœud possède 2 x 6 CPU de modèle Westmere-EP et de référence Intel(R) Xeon(R) CPU E5645 @ 2.40GHz, ainsi que 48GB de RAM, les nœuds sont reliés entre eux grâce à un réseau infiniband de 1GB/s. Les nœuds de calculs sont équipés du système Debian version 6.0.8, de openmpi 1.6. Le code est écrit en C et compilé avec pour flag -g -w -O3 ; le solveur UMFPACK [12, 13] est utilisé pour la factorisation LU de matrices creuses et la résolution des systèmes linéaires. Les calculs sont effectués en arithmétique complexe double précision et la bibliothèque MPI[16, 17] assure la connexion entre le *Maître* et les *Travailleurs*.

L'algorithme parallèle est exécuté sur plusieurs matrices ; dans un premier exemple nous illustrons sur la figure 6 la construction de la chaîne de triangles en démarrant dans 6 directions ; dans l'exemple 2, nous illustrons le fait que les processus lancés dans ces 6 directions peuvent découvrir plusieurs composantes du pseudo spectre ; ceci est illustré par la figure 7. Les figures 8 et 9 montrent l'accélération et l'efficacité de l'algorithme
parallèle sur trois matrices.

### 4.3.1. Premier cas test du programme parallèle

Le premier cas test concerne la matrice BCSSTM25 qui est issue d'un problème gé- néralisé de valeurs propres. Cette matrice est tirée de la collection de Harwell-Boeing de Matrix Market[14] sa taille est $n = 15439$, sa norme 2 vaut $7 \times 10^8$ et son conditionnement $6.1 \times 10^9$, pour un point de référence $z_{ref} = 0.1 + 0.1i$ une taille de maillage $\tau = 10^{-3}$ et un niveau de ligne $\eta = 10^{-2}$. La figure 6 présente une orbite de 150 triangles, la ligne polygonale obtenue en considérant les points extérieurs des triangles contient 6 valeurs propres et le nombre d'intervalle à la fin de l'intégration est 494.

<figure><img src="image_5.png" /><figcaption>**Figure 6.** : Ligne de niveau construite pour la matrice BCSSTM25 (les étoiles rouge représentent les points nécessaires à l'intégration)</figcaption></figure>

### 4.3.2. Deuxième cas test du programme parallèle

Nous illustrons dans ce cas test la découverte de plusieurs composantes du pseudo-spectre lorsque l'on démarre la construction dans 6 directions la matrice considérée est FS-680-1, une matrice issue de l'étude d'un problème de radiations chimiques elle est également tirée de la collection Harwell-Boeing da Matrix Market[14], sa taille est $n = 680$, sa norme 2 vaut $7.2 \times 10^{13}$ et son conditionnement $2.1 \times 10^4$, pour un point de référence $z_{ref} = 7 \times 10^9$ une taille de maillage $\tau = 10^5$ et un niveau de ligne $\eta = 10^6$. Deux composantes du pseudo spectre sont construites :

- une première orbite de 646 triangles dont l'intégration sur les points externes donne 6 valeurs propres pour un nombre total de points insérés égale à 2677;

- une deuxième de 1250 triangles dont l'intégration sur les points extérieurs donne 488 valeurs propres avec un nombre de points égale à 12893.

La figure 7 présente ces résultats.
<figure><img src="image_6.png" /><figcaption>Figure 7. : Ligne de niveau construite pour la matrice FS-680-1</figcaption></figure>

### 4.3.3. Autres cas test

Pour illustrer le gain obtenu en parallélisant l'algorithme combiné, nous avons effectué des tests avec des matrices dont les caractéristiques sont décrites dans le tableau 2. Les

<table><caption>Tableau 2. Caractéristiques de quelques matrices tests.</caption><thead><tr><td>Matrice</td><td>Origine</td><td>Type</td><td>n</td><td>nZ</td></tr></thead><tbody><tr><td>E40R5000</td><td>modélisation d'écoulement</td><td>réelle non symétrique</td><td>17281</td><td>553956</td></tr><tr><td>CHEVRON1</td><td>modélisation sismique</td><td>complexe non hermitienne</td><td>37365</td><td>330633</td></tr><tr><td>CIRCUIT4</td><td>modélisation électrique</td><td>réelle non symétrique</td><td>80209</td><td>307604</td></tr></tbody></table>

zones d'intérêt ($z_r, e, f$), la résolution du maillage ($\tau$), le niveau de la ligne ($\eta$), le nombre de triangles, le nombre de sommets de la ligne polygonale et le nombre de valeurs propres encerclées sont présentés dans le tableau 3. Le tableau 4 présente les temps de calcul, l'accélération et l'efficacité pour les matrices test du tableau 2. De ces données, il ressort que, au dessus de 7 processeurs (84 cœurs), on commence à perdre en accélération et efficacité ; ceci est dû non seulement aux tranches incomplètes mais également au temps de communication qui augmente avec le nombre de processeurs(cœurs).

Les courbes de la figure 8 et 9 illustrent l'accélération et l'efficacité obtenue pour les matrices du tableau 2 en fonction du nombre de processeurs employés de 1 à 10. Chaque processeur a 12 cœurs et les calculs sont effectués sur des cœurs.
<figure><img src="image_7.png" /><figcaption>Figure 8. Accélération première série de tests par rapport à 120 cœurs</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 9. Efficacité première série de tests</figcaption></figure>
<table><caption>Tableau 3. Nombre de valeurs propres par zone d'intérêt.</caption><thead><tr><th>Matrice</th><th>z<sub>ref</sub></th><th>τ</th><th>η</th><th>Lignes</th><th>Triangles</th><th>Sommets</th><th>Valeurs P</th></tr></thead><tbody><tr><td>e40r5000</td><td>1.7 + 1.1i</td><td>0.1</td><td>10<sup>-6</sup></td><td>1</td><td>6</td><td>1516</td><td>2</td></tr><tr><td>Chevron1</td><td>-0.039</td><td>10<sup>-4</sup></td><td>10<sup>-4</sup></td><td>1</td><td>48</td><td>1377</td><td>2</td></tr><tr><td>Circuit4</td><td>0</td><td>10<sup>-4</sup></td><td>10<sup>-4</sup></td><td>2</td><td>122(278)</td><td>11243(78463)</td><td>128(11507)</td></tr></tbody></table>

## 5. Conclusion

Dans cet article, nous avons développé un algorithme parallèle en combinant deux approches : PAT qui construit une orbite de triangles équilatéraux en suivant une ligne de niveau $η$ et la méthode EIGENCNT, qui compte le nombre de valeurs propres entourées par une ligne polygonale a priori donnée dans le plan complexe. Cet algorithme réduit le nombre de factorisations LU en exploitant les factorisations LU obtenues lors du calcul de $σ_{min}$ pour calculer la dérivée. L'utilisation d'un modèle *Client-Serveur* pour l'implémentation de cet algorithme permet d'optimiser l'équilibrage de charges entre les différents processeurs d'un réseau de machines, ainsi que de réduire la complexité des données échangées. Les communications sont réduites à un échange de messages courts, constitués de quelques nombres complexes entre le *Maître* et les *Travailleurs* permettant d'utiliser des machines géographiquement éloignées.

Les bonnes performances de cet algorithme nous permettent de dénombrer les valeurs propres de grandes matrices dans une région du plan complexe en construisant des lignes polygonales faisant intervenir un nombre de factorisation LU très important pour plus de valeurs propres à dénombrer.

Malgré le fait que cet algorithme permette une localisation robuste de valeurs propres de matrice dans le plan complexe, un deuxième niveau de parallélisme doit être considéré pour paralléliser le calcul de déterminant et/ou d'une dérivée et même de $σ_{min}$. Cela serait possible en considérant des matrices bloc-diagonales comme dans [4].

## 6. Bibliographie

[1] O. Bertrand et B. Philippe. Counting the eigenvalues surrounded by a closed curve. *Siberian Journal of Industrial Mathematics*, 4:73-94, 2001.

[2] D. Bindel, Bounds and error estimates for nonlinear eigenvalue problems. Berkeley Applied Mathematical Seminar, October 2008. D. Bindel, Bounds and error estimates for nonlinear eigenvalue problems. Berkeley Applied Mathematical Seminar, October 2008.

[3] S. Godunov. Spectral portraits of matrices and criteria of spectrum dichotomy. In L. Atanasova and J. Hezberger, editors, *Third Int'l. IMACS-CAMM Symposium on Computer Arithmetic and Enclosure Methods*, Amsterdam, 1992.

[4] E. Kamgnia et L. B. Nguenang. Some efficient methods for computing the determinant of large matrices. *ARIMA Journal*, Special issue CARI'12, 17:73-92, 2014.

[5] E. Kamgnia et B. Philippe. Counting eigenvalues in domains of the complex field. *ETNA*, 40:1-16, 2013.

[6] D. Mezher et B. Philippe. Parallel computation of pseudospectra of large sparse matrices. *Parallel Comput.*, 28(2):199-221, 2002.

[7] D. Mezher et B. Philippe. PAT : A reliable path following algorithm. *Numer. Algorithms*, 29:131-152, 2002.

[8] E. Polizzi et A. Sameh. A parallel hybrid banded systems solver : the SPIKE algorithm. *Parallel Comput.*, 32:177-194, 2006.
<table><caption>**Tableau 4.** *Temps(en s), accélération et efficacité pour e40r5000, Chevron1 et Circuit4.*</caption><thead><tr><th rowspan="2">Proc</th><th colspan="3">E40R5000</th><th colspan="3">CHEVRON1</th><th colspan="3">CIRCUIT4</th></tr><tr><th>Temps</th><th>Acc.</th><th>Eff.</th><th>Temps</th><th>Acc.</th><th>Eff.</th><th>Temps</th><th>Acc.</th><th>Eff.</th></tr></thead><tbody><tr><td>1(12)</td><td>1475.86</td><td>1.0</td><td>1.0</td><td>4831.73</td><td>1.0</td><td>1.0</td><td>56095.59</td><td>1.0</td><td>1.0</td></tr><tr><td>2(24)</td><td>729.11</td><td>2.02</td><td>0.96</td><td>2379.44</td><td>2.03</td><td>0.97</td><td>26828.33</td><td>2.09</td><td>1.0</td></tr><tr><td>3(36)</td><td>501.09</td><td>2.94</td><td>0.92</td><td>1625.50</td><td>2.97</td><td>0.93</td><td>18636.98</td><td>3.00</td><td>0.94</td></tr><tr><td>4(48)</td><td>389.13</td><td>3.79</td><td>0.88</td><td>1250.23</td><td>3.86</td><td>0.90</td><td>14655.26</td><td>3.82</td><td>0.89</td></tr><tr><td>5(60)</td><td>327.19</td><td>4.51</td><td>0.84</td><td>1025.75</td><td>4.71</td><td>0.87</td><td>12269.34</td><td>4.57</td><td>0.85</td></tr><tr><td>6(72)</td><td>280.76</td><td>5.25</td><td>0.81</td><td>877.10</td><td>5.50</td><td>0.85</td><td>10827.51</td><td>5.18</td><td>0.80</td></tr><tr><td>7(84)</td><td>241.97</td><td>6.09</td><td>0.80</td><td>765.03</td><td>6.31</td><td>0.83</td><td>9678.83</td><td>5.79</td><td>0.76</td></tr><tr><td>8(96)</td><td>226.44</td><td>6.51</td><td>0.75</td><td>690.19</td><td>7.00</td><td>0.81</td><td>8810.29</td><td>6.36</td><td>0.73</td></tr><tr><td>9(108)</td><td>196.84</td><td>7.49</td><td>0.77</td><td>686.09</td><td>7.04</td><td>0.72</td><td>8219.26</td><td>6.82</td><td>0.70</td></tr><tr><td>10(120)</td><td>187.64</td><td>7.86</td><td>0.72</td><td>683.78</td><td>7.06</td><td>0.65</td><td>7713.35</td><td>7.27</td><td>0.67</td></tr></tbody></table>

[9] Y. Maeda, Y. Futamura, AND T. Sakurai, Stochastic estimation method of eigenvalue density for nonlinear eigenvalue problem on the complex plane, JSIAM Letters, 3:61-64, 2011.

[10] G. W. Stewart and Ji-guang Sun. Matrix Perturbation Theory Academic Press Inc. page 174, 1990.

[11] L. Trefethen. Pseudospectra of Matrices. In D. F. Griffiths and G. A. Watson, editors, *Numerical Analysis*, pages 234–266. Longman, 1992.

[12] UMFPACK; Unsymmetric Multifrontal Sparse LU Factorization Package; http://www.cise.ufl.edu/research/sparse/umfpack/

[13] T. A. Davis Algorithm 832: UMFPACK, an unsymmetric-pattern multifrontal method *ACM Transactions on Mathematical Software*, 30:2, 199-200, 2004.

[14], Matrix Market, *Service of the Mathematical and Computational Sciences Division / Information Technology Laboratory / National Institute of Standards and Technology*, http://math.nist.gov/MatrixMarket/.

[15], T. A. Davis and Y. Hu The University of Florida Sparse Matrix Collection, *ACM Transactions on Mathematical Software*, Vol 38, Issue 1, 2011, pp 1:1 - 1:25, http://www.cise.ufl.edu/research/sparse/matrices.

[16], Naoyuki Shida and Shinji Sumimoto and Atsuya Uno MPI Library and Low-Level Communication on the K computer, *FUJITSU Scientific & Technical Journal*, Vol 48, 2012, pp 3:324 - 3:330,

[17], MPI Forum MPI: A Message-Passing Interface Standard, *Message Passing Interface Forum*, http://www.mpi-forum.org/