r6.02.01 A propos des méthodes de décomposition de type GAUSS#
Résumé :
Ce document présente quelques aspects liés à la méthode de décomposition de GAUSS.
Après une rapide présentation, nous rappelons les principaux avantages et inconvénients liés à cette méthode directe. Puis, nous détaillons l’implémentation de l’algorithme LDLT mis en oeuvre dans le Code Aster .
GAUSS (1777 - 1855) est à l’origine de toutes les méthodes directes de résolution numérique de systèmes linéaires. Qu’il en soit ici remercié.
Table des matières
Inconvénients des méthodes de type GAUSS#
Les inconvénients des méthodes de type GAUSS sont essentiellement de trois types :
un nombre élevé d’opérations
un remplissage de la matrice
une perte de précision en cours de calcul
Les deux premiers points sont souvent qualifiés de défauts majeurs alors que le troisième est considéré comme un défaut mineur.
Le nombre d’opérations#
Pour un système plein de taille n, à la p-ième étape, nous devons effectuer pour calculer les nouveaux coefficients de la matrice et le second membre :
(n-p) divisions
(n-p+1)(n-p) additions et multiplications
Le nombre d’opérations est donc :
\(\sum_{p=1}^{n-1}(n-p)=\frac{n(n-1)}{2}\) divisions
\(\sum_{p=1}^{n-1}(n-\text{p+1})(n-p)=\sum_{q=1}^{n-1}{q}^{2}+\sum_{q=1}^{n-1}q=(\frac{n(n-1)(n-2)}{6}+\frac{n(n-1)}{2})\) additions et autant de multiplications.
Soit \(\frac{1}{3}\) \(n(n-1)\) \((n+\frac{1}{2})\) opérations auxquelles il convient d’ajouter les \({n}^{2}\) opérations de la résolution du système triangulaire.
En résumé : |
Pour un grand système plein, l’algorithme de GAUSS nécessite de l’ordre de
|
Remarque :
Dans le cas d’une matrice stockée bande, le nombre d’opérations est \({\mathrm{n.b}}^{2}\) où \(b\) est la largeur de la bande.
Le remplissage de la matrice#
Commençons par un exemple classique de matrice dite « flèche » que l’on rencontre par exemple en chimie [bib9].
Soit \(A\) la matrice telle que \(a(1,i)\ne 0,a(i,1)\ne 0,a(i,i)\ne 0\) et tous ses autres termes sont nuls, la matrice a alors l’allure suivante :
A =
Après la 1ère étape de factorisation (par l’algorithme de GAUSS), la matrice est pleine au sens où il n’y a plus de termes théoriquement nuls.
Regardons plus formellement le phénomène de remplissage du à l’algorithme ; pour cela récrivons l’algorithme :
Pour k variant de 1 à n - 1 faire % boucle sur les étapes
Pour i variant de k + 1 à n faire % boucle sur les lignes
Pour j variant de k + 1 à n faire % boucle sur les colonnes
\(a(i,j)=a(i,j)-a(i,k).a(k,j)/a(k,k)\)
fin faire
fin faire
fin faire
Ce que l’on peut schématiser graphiquement, à la k-ième étape, pour le calcul du terme a(i,j) par :
Le terme \(a(i,j)\) est non nul à la fin de la k-ième étape :
s’il était non nul au début de cette k-ième étape,
ou si les termes \(a(i,k)\) et \(a(k,j)\) sont tous les deux non nuls au début de la k-ième étape, et ceci indépendamment de la valeur initiale du terme \(a(i,j)\) .
De plus, on voit que la méthode de GAUSS remplit le profil enveloppe de la matrice au cours des étapes de factorisation.
Dans l’exemple de la matrice « flèche » : le profil enveloppe est la matrice pleine, d’où le résultat constaté.
Cet exemple met en évidence l’importance de la numérotation des inconnues de la matrice puisque la matrice peut être récrite, après permutation des inconnues, sous la forme :
dont le profil enveloppe est la matrice elle même (il n’y a donc pas de remplissage).
Nous venons de voir l’importance de la numérotation des inconnues.
Nous ne saurions trop insister sur le fait que des algorithmes de renumérotation « optimale » doivent être utilisés pour minimiser le remplissage de la matrice.
Ces algorithmes reposent sur des heuristiques et sont spécialisés.
Parmi les algorithmes les plus couramment utilisés citons
Algorithmes |
Objectifs |
CUTHILL - Mc KEE |
minimiser la largeur de la bande |
Reverse CUTHILL - Mc KEE |
minimiser le profil |
Minimum Degré |
minimiser les multiplications par 0 |
Formulation intrinsèque du remplissage
On peut donner une formulation intrinsèque du remplissage lors de l’élimination d’inconnue en termes de graphe [bib5].
Soit une matrice A à laquelle nous associons le graphe \(G(X,E)\) , où \(X\) est l’ensemble des noeuds et \(E\) l’ensemble des arêtes non orientées.
Le problème d’élimination d’une inconnue de la matrice est alors équivalent à éliminer un noeud du graphe.
Définition : Soient \(x,y\in X\) , on dira que \(x\) et \(y\) sont adjacents si et seulement si \(x,y\in E\)
Si \(Y\) est un sous-ensemble de nœuds de \(X(X\supset Y)\) nous pouvons définir les ensembles suivants :
l’ensemble des noeuds adjacents à Y |
\(\mathrm{Adj}(Y)=x/x\in X-Y\mathrm{et}y\in Y\mathrm{tel}\mathrm{que}x,y\in E\) |
l’ensemble des cotés incidents à Y |
\(\mathrm{Inc}(Y)=x,y/y\in Y,x\in \mathrm{Adj}(Y)\) |
l’ensemble de définition de X |
\(\text{Def}(X)=y,z/y,z\in \mathrm{Adj}(X),y\ne z\mathrm{et}z\notin \mathrm{Adj}(Y)\) |
Opération d’élimination de y
L’élimination de y consiste alors à considérer le sous-ensemble
\(\mathrm{Elim}(y,G)=X-y,(E-\mathrm{Inc}(y))U\stackrel{\scriptscriptstyle\mathrm{def}}{=}(y)\)
Il faut donc bien considérer le remplissage qui est lié à \(\text{Def}(y).\)
Pour minimiser le remplissage on peut utiliser une heuristique consistant à éliminer le y de degré minimal, le degré de {y} étant le cardinal de l’ensemble Adj(y). C’est l’idée maîtresse de l’utilisation de l’algorithme du minimum degré avant une factorisation de type GAUSS.
Cette approche est utilisée dans la méthode multi-frontale mise en oeuvre dans Aster [bib15].
La perte de précision en cours de calcul#
Le problème provient du fait qu’en cours d’algorithme les pivots décroissent et qu’ils sont utilisés comme dénominateur pour les étapes suivantes [bib13].
Etude sommaire de la perte de précision#
Notons \({A}^{(k)}\) la matrice à l’étape k (c’est-à-dire après l’élimination de la k-ième variable) ; avec par convention \({A}^{(o)}=\mathrm{A.}\)
Nous pouvons alors écrire que A(1) vérifie \({L}^{(1)}.{A}^{(1)}={A}^{(0)}\)
en prenant L(1) =
par identification
Puis après (n-1) factorisations de ce type nous obtenons
\(({L}^{(n-1)}•\dots •{L}^{(1)})U=\mathrm{L.U}={A}^{(0)}\)
Du fait des erreurs E sur la décomposition, il convient d’écrire : \(\mathrm{L.U}={A}^{(0)}+E\)
Les coefficients de la matrice d’erreur E peuvent être évalués [bib12] en tenant compte de l’erreur commise sur l’opération flottante que nous noterons fl :
\(\mid {e}_{ij}\mid \le 3.\varepsilon .{m}_{ij}.\underset{\text{k,i,j}}{\max}\mid {a}_{ij}^{(k)}\mid {\forall}_{\text{i,j}}\) \(\forall i,j\)
avec \({m}_{ij}:\) |
nombre de termes tel que \({a}_{\text{ik}}^{(k)}.{a}_{\text{kj}}^{(k)}\ne 0\) |
\(\varepsilon\) : |
erreur relative des opérations machine |
En effet, en se plaçant dans le cas où l’élimination de GAUSS « réussit » (par exemple : la matrice est définie positive ou bien on utilise une technique de pivotage).
Notons \({\eta}_{\text{ik}}=\beta ()\equiv ().(1+{\varepsilon}_{1})\) pour \(i>k,k=1,\dots ,n-1\)
avec \(\mid {\varepsilon}_{1}\mid <\varepsilon\)
Le terme
est alors évalué par :
\({a}_{ij}^{(k)}=\beta ({a}_{ij}^{(k-1)}-{\eta}_{\text{ik}}.{a}_{\text{kj}}^{(k-1)})\) pour \(i,j>k,k=1,\dots ,n-1\)
Soit encore \({a}_{ij}^{(k)}=\beta ({a}_{ij}^{(k-1)}-{\eta}_{\text{ik}}.{a}_{\text{kj}}^{(k-1)})\) avec \(\mid {\varepsilon}_{2}\mid ,\mid {\varepsilon}_{3}\mid <\varepsilon\)
La « perturbation » \({e}_{ij}^{k}\) subie par \({a}_{ij}^{(k)}\) peut alors être évaluée ; à partir de la définition de \({m}_{ij}\) nous en déduisons la relation :
\(\mid {e}_{ij}^{k}\mid \equiv \mid {\varepsilon}_{1}.{a}_{ij}^{(k-1)}\mid \le \varepsilon .\mid {a}_{ij}\mid\) pour \(i>k,k=1,\dots ,n-1\)
et de l’évaluation première de \({a}_{ij}^{(k)}\) nous déduisons :
\(\mid {\mu}_{\text{ik}}.{a}_{\text{kj}}^{(k-1)}\mid \equiv ({a}_{ij}^{(k-1)}-{a}_{ij}^{(k)}/(1+{\varepsilon}_{3}))/(1+{\varepsilon}_{2})\)
et finalement
\(\mid {e}_{\text{kj}}^{(k)}\mid \equiv \mid {a}_{ij}^{(k)}(1-\frac{1}{(1+{e}_{2})(1+{e}_{3})})-{a}_{ij}^{(k-1)}(1-\frac{1}{(1+{e}_{2})})\mid\)
\(\mid {e}_{\text{kj}}^{(k)}\mid \le 3.\mid {a}_{ij}\mid .\varepsilon\)
or la décomposition \(\mathrm{L.U}=A(0)+E\) nous indique que \({e}_{ij}\) est la somme des erreurs.
Estimation de l’erreur sur la solution#
Nous avons résolu le système
\(\mathrm{L.}Ux=A(0).x+\mathrm{E.x}\)
où |
\(\mathrm{E.x}\) est le terme d’erreur induit par les erreurs d’arrondi/troncature dans les opérations de la factorisation. |
\(A(0)x\) est le second membre (en fait \(b\) ). |
La solution approchée trouvée \(\tilde{x}\) qui approxime la vraie solution x est :
\(\tilde{x}\) \(={(A+E)}^{\text{-1}}b\)
On montre alors que l’évaluation de l’erreur sur x est liée au conditionnement de A.
Posons le problème sous la forme : \((A+\delta A)(x+\delta x)=b\)
En supposant \(\delta \mathrm{A.}\delta x\) petit on a : \(\delta x=A-1.\delta \mathrm{A.x}\)
d’où en normant \(\frac{\parallel \mathrm{\delta \; x}\parallel }{\parallel x\parallel }\le \text{Cond}(A).\frac{\parallel \mathrm{\delta \; A}\parallel }{\parallel A\parallel }\)
où \(\mathrm{Cond}(A)\) = \((\parallel A\parallel .{\parallel A\parallel }^{-1})\) est le conditionnement de la matrice A.
Remarque:
On constate que l’erreur induite sur le second membre est faible et ne perturbe la solution qu’à travers un mauvais conditionnement de la matrice A.
En effet, si l’on considère le système \(A(x+\delta x)=b+\delta b,\)
du fait des égalités : \(\mathrm{\delta x}=A-1.\mathrm{\delta b}\mathrm{et}\mathrm{Ax}=b\)
\(\begin{array}{}\text{on a}\text{}\parallel \mathrm{\delta \; x}\parallel \le \parallel {A}^{\text{-1}}\parallel .\parallel \mathrm{\delta \; b}\parallel \text{et}\parallel b\parallel <\parallel A\parallel .\parallel x\parallel \\ \text{et donc}\text{}\parallel \frac{\mathrm{\delta \; x}}{x}\parallel \le (\parallel A\parallel .\parallel {A}^{\text{-1}}\parallel )\parallel \frac{\mathrm{\delta \; b}}{b}\parallel \end{array}\)
Remarque:
Si l’on considère les variations sur \(A,x,\) et b simultanément, on a l’estimation suivante [bib14]:
\(\frac{\parallel \mathrm{\delta \; x}\parallel }{\parallel x\parallel }\le \frac{\text{Cond}(A)}{1-\parallel {A}^{-1}\parallel .\parallel \mathrm{\delta \; A}\parallel }.(\frac{\parallel \mathrm{\delta \; A}\parallel }{\parallel A\parallel }+\frac{\parallel \mathrm{\delta \; b}\parallel }{\parallel b\parallel })\)
Quelques remarques sur le conditionnement
le conditionnement n’est défini que pour une matrice régulière,
le conditionnement dépend de la norme choisie sur \({ℝ}^{n}\) ,
quelque soit la norme choisie, nous avons \(1\le \mathrm{Cond}(A)\) et une matrice est d’autant mieux conditionnée que son conditionnement est proche de 1.
Si la norme euclidienne est choisie pour norme, alors
le conditionnement d’une matrice A quelconque est
\(\mathrm{Cond}(A)=\) \(\mid \frac{{\mu}_{n}}{{\mu}_{1}}\mid\)
où \({\mu}_{1}\) et \({\mu}_{n}\) sont les valeurs singulières extrêmes de A (c’est-à-dire la plus petite et la plus grande des valeurs propres de \({A}^{\text{*}}.A\) ).
Si la matrice A est symétrique (ou hermitienne) alors
\(\mathrm{Cond}(A)\) = \(\mid \frac{{\lambda}_{n}}{{\lambda}_{1}}\mid\)
où \({\lambda}_{1}\) et \({\lambda}_{n}\) sont les valeurs propres de module minimal et maximal de \(A\) .
Estimation du nombre de chiffres significatifs de la solution#
Si l’on possède une précision de p chiffres (décimaux) significatifs, on a alors :
\(\frac{\parallel \mathrm{dA}\parallel }{\parallel A\parallel }\simeq {10}^{-p}\)
Si l’on désire une précision de s chiffres (décimaux) significatifs sur la solution
\(\frac{\parallel \deltax \parallel }{\parallel x\parallel }\le {10}^{-s}\)
d’où l’estimation du nombre de chiffres significatifs décimaux exacts de la solution
\(s\ge p-{\log}_{10}(\mathrm{Cond}(A))\)
Méthode pour réduire le conditionnement#
La méthode la plus simple est celle de la mise à l’échelle :
On « passe » de A à \({\varphi}_{1}\text{.A.}{\varphi}_{2}\) avec \({\varphi}_{i}\) matrice diagonale telle que Cond(\({\varphi}_{1}\text{.A.}{\varphi}_{2}\) ) soit meilleur que \(\mathrm{Cond}(A)\) .
Ceci est très théorique et il n’existe pas de méthode universelle pour déterminer \({\varphi}_{1}\text{et}{\varphi}_{2}\) .
Notons que si la matrice A est symétrique et que l’on désire conserver cette propriété, il faut alors prendre \({\varphi}_{1}={\varphi}_{2}\) .
Exemple de matrice mal conditionnée.#
Cet exemple, très significatif et instructif est du à R.S. WILSON.
Soit le système \(\mathrm{Ax}=b\) avec :
\(A\) = \((\begin{array}{cccc}10& 7& 8& 7\\ 7& 5& 6& 5\\ 8& 6& 10& 9\\ 7& 5& 9& 10\end{array})\) et \(b\) = \((\begin{array}{c}32\\ 23\\ 33\\ 31\end{array})\) et dont la solution est \(x\) = \((\begin{array}{c}1\\ 1\\ 1\\ 1\end{array})\)
Si l’on perturbe le second membre de l’ordre de 0.5 % en prenant :
\(\tilde{b}=(\begin{array}{cccc}32.1,& 22.9,& 33.1,& 30.9\end{array})\)
alors la solution est : \(\tilde{x}=(\begin{array}{cccc}9.2,& \text{-12}.6,& 4.5,& \text{-1}.1\end{array})\) .
Si la matrice est perturbée de l’ordre de 1 % :
\(Ã\) = \((\begin{array}{cccc}10.& 7& 8.1& 7.2\\ 7.08& 5.04& 6& 5\\ 8& 5.98& 9.89& 9\\ 6.99& 4.99& 9& 9.98\end{array})\)
alors la solution est \(\tilde{x}=(\begin{array}{cccc}\text{-81,}& 137,& \text{-34,}& 22\end{array})\)
Remarques sur les propriétés de la matrice A:
Elle est symétrique, définie positive, de déterminant 1et d’inverse « sympathique ».
\({A}^{\text{-1}}\) = \((\begin{array}{cccc}25& \text{-41}& 10& \text{-6}\\ \text{-41}& 68& \text{-17}& 10\\ 10& \text{-17}& 5& \text{-3}\\ \text{-6}& 10& \text{-3}& 2\end{array})\)
Son conditionnement au sens de la norme euclidienne est :
\(\text{Cond}(A)=\frac{{l}_{4}}{{l}_{1}}=\frac{30.2887}{0.01015}\text{=2984}.11\)
Une Interprétation géométrique du mauvais conditionnement#
On peut donner une interprétation très simple du mauvais conditionnement d’un système linéaire \(\mathrm{Ax}=b\) dans le cas particulier où A est normale (i.e. \({A}^{\text{*}}.A={\mathrm{A.A}}^{\text{*}}\) ).
Soient \({\lambda}_{1}\) la plus petite valeur propre de la matrice \(A\) et \({\lambda}_{n}\) sa plus grande valeur propre et soient \({v}_{1}\) et v_n les vecteurs propres associés.
pour \(b={v}_{n}\text{et}\deltab ,\) on a \(\parallel \frac{\deltax }{x}\parallel =\text{cond}(A)\parallel \frac{\deltab }{b}\parallel\)
d’où si \(\mathrm{cond}(A)\) = \(\mid \frac{{\lambda}_{n}}{{\lambda}_{1}}\mid\) est grand, une petite perturbation \(\deltab\) de \(b\) entraîne une grande variation sur la solution \(v\) .
Critères de détermination d’un pivot nul#
Définition : |
Un système numériquement dégénéré est un système pour lequel un pivot est nul ou n’a pas de chiffre significatif exact. |
Notons qu’un système peut être dégénéré numériquement sans l’être mathématiquement.
Dans ces deux cas, il convient de ne pas résoudre le système d’où la nécessité de déterminer un critère d’arrêt dès que l’un des pivots n’a plus de chiffres significatifs exacts.
Soient A la matrice à factoriser et F la matrice de diagonale libre issue de la factorisation.
Critère 1 : |
Le critère le plus simple est de considérer que le pivot est nul dès qu’il est inférieur, en valeur absolue à un seuil donné. |
\(\mid {f}_{ii}\mid <{\varepsilon}_{1}\)
\({\varepsilon}_{1}\) est un « nombre petit » en dessous duquel on considère que les valeurs sont arbitrairement nulles.
Critère 2 : |
Ce critère s’applique au nombre de chiffres significatifs encore disponibles. En constatant que l’on ne peut avoir plus de p chiffres significatifs sur une machine donnée, on vérifiera que la décroissance du pivot ne s’effectue pas dans un rapport supérieur à 10-p. |
\(\mid \frac{{f}_{ii}}{{a}_{ii}}\mid \le {\varepsilon}_{2}={10}^{\text{-p}}\)
Notons que le rapport \(\mid \frac{{f}_{ii}}{{a}_{ii}}\mid\) est toujours inférieur ou égal à 1 du fait de la décroissance des pivots.
La cause de base d’un mauvais conditionnement numérique est l’erreur d’arrondi provoquée par l’introduction de grands nombres sans signification physique.
Méthode LDLT par blocs mise en oeuvre dans Aster#
Ce paragraphe détaille la mise en oeuvre dans Aster de la résolution du système linéaire A.x = b par la méthode de factorisation LDLT de la matrice symétrique A.
La matrice A est stockée en profil (ou ligne de ciel) par bloc.
Principe de base : |
Une colonne de la matrice est contenue toute entière dans un seul bloc : nous ne segmentons pas les colonnes. |
Les tableaux permettant la description de la matrice stockée profil par bloc sont :
HCOLhauteur de colonne de la matrice
HCOL(i)hauteur de la i-eme colonne
ADIAadresse du terme diagonal dans son bloc
ADIA(i)renvoie l’adresse du i-eme terme diagonal dans son bloc
ABLOpointeur de bloc
ABLO(i+1) renvoie le numéro de la dernière équation dans la numérotation globale contenu dans le i-eme bloc
Par convention ABLO(1) = 0 et le nombre d’équations dans le i-eme bloc est donné par la relation ABLO(i+1) - ABLO(i) + 1
Le nombre total d’équation se déduit comme étant ABLO(nombre_de_bloc + 1)
Il est également nécessaire de mémoriser le nombre total de blocs utilisés pour contenir les coefficients de la matrice.
Remarque:
Formellement le tableau HCOLest inutile car il se déduit des tableaux ADIA et ABLO, mais il permet d’effectuer les calculs plus rapidement.
Mise en oeuvre de la factorisation#
Les principales particularités de la mise en oeuvre de la factorisation de GAUSS par la variante de CROUT sous forme LDLT d’une matrice symétrique stockée profil par bloc sont :
la factorisation est effectuée en place, c’est à dire en écrasant la matrice initiale,
la factorisation peut-être partielle,
les critères de détections de pivot nuls peuvent être adaptés à la factorisation de matrices quasi-singulières,
en cas de détection de pivot nul, ce pivot est remplacé par une très grande valeur (1040) ce qui revient à introduire une condition de blocage par pénalisation.
Remarque:
On crée deux tableaux de travail :
un tableau qui contiendra la diagonale de la matrice factorisée (minimisation du nombre d’accès au bloc),
un tableau pour la colonne courante (minimisation du nombre de calculs effectués).
DEBUT Algorithme ;
création d’un tableau intermédiaire pour la colonne courante
création d’un tableau intermédiaire pour la colonne courante
POUR ibloc VARIANT_DE 1 A nombre_de_bloc FAIRE
Détermination des colonnes de départ et de fin pour le bloc courant.
Recherche de la plus petite équation en relation avec une équation contenue dans le bloc courant. Cette recherche se fait en exploitant le tableau HCOL
Recherche du bloc d’appartenance de l’équation trouvé précédemment. Cette recherche se fait en exploitant le tableau ABLO.
Requête en mode écriture du i-eme bloc.
POUR jbloc VARIANT_DE plus_petit_concerne A ibloc-1 FAIRE
Requête en mode lecture du j-eme bloc
POUR iequa CONTENUE DANS LE i-eme bloc FAIRE
calcul de l’adresse de départ de la colonne dans le bloc
calcul de la hauteur de la colonne
POUR jequa CONTENUE DANS LE j-eme bloc FAIRE
calcul de l’adresse de départ de la colonne dans le bloc
calcul de la longueur de la colonne
A(ibloc,jequa) = A(ibloc,jequa)- < A(ibloc,*), A(jbloc,*) >
FIN_POUR
FIN_POUR
libération du j-eme bloc qui n’a pas été modifié
FIN_POUR
POUR iequa CONTENUE DANS LE i-eme bloc FAIRE
calcul de l’adresse de départ de la colonne dans le bloc
calcul de la longueur de la colonne
POUR jequa CONTENUE DANS LE i-eme bloc et < iequa FAIRE
calcul de l’adresse de départ de la colonne dans le bloc
calcul de la hauteur de la colonne
A(ibloc,lm) = A(ibloc,lm) - < A(ibloc,*) , A(jbloc,*) >
FIN_FAIRE
% utilisation de la colonne iequa (calcul du pivot)
calcul de l’adresse de départ de la colonne dans le bloc
calcul de la hauteur de la colonne
sauvegarde de la colonne: trav(i) \(\neg\) A(ibloc,i)
normalisation de la colonne en utilisant le tableau diagonal :
A(ibloc,*)
A(ibloc,*) / diag(*)
calcul du terme diagonal et actualisation du tableau de travail:
tabr8(iadia) = tabr8( iadia ) - < A(ibloc,*) , trav(*) >
test du pivot par rapport à \(\varepsilon\)
test du pivot par rapport au nombre de chiffres significatifs
FIN_POUR
libération du bloc courant que l’on vient de modifier
FIN_POUR
libération des tableaux de travail
FIN Algorithme ;
Détermination de la colonne de départ et de fin pour le bloc courant.
Cet phase est due à la notion de factorisation partielle.
SI (dernière colonne du bloc < début de factorisation) ALORS
requête en mode lecture du i-eme bloc remplir le tableau de travail contenant la diagonale. libération du i-eme bloc ALLER_AU bloc suivant
SINON_SI (première colonne du bloc > fin de factorisation) ALORS
SORTIR
SINON
SI (première colonne du bloc < début de factorisation) ALORS % compléter le tableau « diagonal » requête en mode lecture du i-eme bloc remplir le tableau de travail contenant la diagonale. FIN_SI SI (dernière colonne du bloc > fin de factorisation) ALORS modification de la dernière équation à prendre en compte FIN_SI
FIN_SI
Remarque sur l’encombrement :
Il est obligatoire que l’on puisse avoir au moins simultanément en mémoire :
deux blocs de la matrice,
deux vecteurs de travail de taille : le nombre d’équations du système à résoudre.
Mise en oeuvre de la résolution#
La mise en oeuvre de la résolution simultanée de n seconds membres du système A.x = b, où la matrice A est symétrique et a été factorisée sous forme L D LT (la résolution est en place)
Remarque :
On crée un tableau de travail qui contiendra la diagonale de la matrice factorisée, en vue de minimiser le nombre accès au bloc de la matrice et donc de limiter les lectures pour les grosses matrices ne pouvant résider totalement en mémoire.
DEBUT Algorithme ;
Création d’un tableau pour stocker la diagonale pour éviter des lectures lors de l’étape de résolution diagonale.
POUR ibloc VARIANT_DE 1 AU nombre_de_bloc |
% résolution descendante et |
% remplissage du tableau diagonal |
requête en mode lecture du i-eme bloc
POUR iequa contenue DANS LE BLOC
Calcul de la hauteur de la colonne et
calcul de l’adresse de départ de la colonne dans son bloc
POUR chaque second membre FAIRE
xsol(isol) = xsol(isol) - < x(isol) , U) > FIN_POUR
sauvegarde du terme diagonal dans le tableau de travail
FIN_POUR
libération du i-eme bloc
FIN_POUR
POUR chaque second membre FAIRE |
% résolution diagonale |
POUR toute les équations FAIRE xsol(iequa,isol) = xsol(iequa,isol) / diag(iequa-1) FIN_POUR
POUR ibloc VARIANT_DE nombre_de_bloc à 1 PAR_PAS_DE -1% résolution remontante
requête en mode lecture du i-eme bloc POUR iequa contenue DANS LE BLOC
Calcul de la hauteur de la colonne et
calcul de l’adresse de départ de la colonne dans son bloc
POUR chaque second membre FAIRE xsol(ixx+i,isol)= xsol(ixx+i,isol)- xsol(isol)*L(ide+i) FIN_POUR FIN_POUR
libération du i-eme bloc
FIN_POUR
Libération de la zone de travail (c’est à dire du tableau diagonal)
FIN Algorithme ;
Remarque sur l’encombrement:
Il est nécessaire que l’on puisse avoir simultanément en mémoire :
un bloc de la matrice,
un vecteur de travail de taille : le nombre d’équations du système à résoudre,
les n seconds membres.
Mise à l’échelle#
Il est possible de faire une mise à l’échelle de la matrice à factoriser ; cette mise à l’échelle simple se fait de façon à obtenir une matrice dont les termes diagonaux valent 1.
La matrice diagonale est telle que :
\({\varphi}_{i}=\lbrace \begin{array}{cc}\frac{1}{\sqrt{\mid {a}_{ii}\mid }}& {\text{si a}}_{ii}\ne 0\\ 1& {\text{si a}}_{ii}=0\end{array}\)
Il est à noter que lors de la résolution, on n’obtient la solution du système qu’après déconditionnement.
En effet : le système initial est \(\mathrm{Ax}=b\)
après multiplication à gauche par \(\varphi\) , on a :
\(\varphi \mathrm{Ax}=\varphi .b\)
Or le système résolu est :
\(\varphia \varphix =\varphib\)
d’où la solution \(\varphix\) obtenue qu’il faut « déconditionner ».
Tests sur le pivot#
Deux critères de détections de pivot nuls sont implémentés :
le test en valeur absolue \(\mid {a}_{ii}\mid <\varepsilon\) , avec \(\varepsilon\) donné,
le test en valeur relative sur le nombre de chiffres significatifs exacts.
Notons que ces tests peuvent être réduits à leur plus simple expression en fournissant un \(\varepsilon =0\) . et en donnant, par convention, un nombre de chiffres significatifs exacts nul.
Cette option est rendue nécessaire par le fait que des algorithmes tels que les algorithmes de recherche de valeurs propres [R5.01.01] [R5.01.02] cherchent à factoriser des matrices quasi‑singulière.
Factorisation de matrices complexes#
L’algorithme implémenté dans Aster permet également de traiter les matrices symétriques à coefficients complexes.
L’algorithme implémenté ne traite pas les matrices hermitiennes, bien que ce soit théoriquement possible.
Méthodes de stockage classiques
A1.1 Matrice pleine
Une matrice pleine non symétrique de taille n possède n2 coefficients.
Si la matrice est symétrique, on peut ne stocker que sa triangulaire inférieure ou supérieure soit \(\frac{n(n+1)}{2}\) valeurs.
Aucun tableau de description de la matrice n’est nécessaire.
A1.2 Matrice bande
b + 1 = largeur de bande
Dans ce cas on stocke la bande (appelée parfois matrice redressée) dans un tableau rectangulaire n (2b + 1) ; on inclut alors les b (b + 1) valeurs nulles correspondant aux compléments des pointes.
Dans le cas d’une matrice symétrique, on peut ne stocker que n (b + 1) valeurs dont \(\frac{b(b+1)}{2}\) valeurs nulles (inutiles).
Cette méthode nécessite uniquement de connaître la largeur de la bande.
A1.3 Matrice profil ou matrice à ligne de ciel
Cette technique consiste à stocker les termes de la matrice par colonnes et lignes de longueurs variables. Les termes extérieurs à la « ligne de ciel », qui est l’enveloppe des sommets des colonnes étant supposés n’avoir aucune contribution dans les calculs, ne sont pas stockés.
Le profil de la i-ème ligne (resp. colonne) est déterminé par :
min {j tel que 1 j n aij:math:ne`0} (resp min {j tel que 1 j n aji:math:ne`0})
Si le profil est symétrique, on parle de matrice à profil symétrique.
Cette méthode de stockage nécessite des tableaux de stockage que nous allons détailler dans le cas d’une matrice à profil symétrique.
Classiquement, dans cette option de stockage, la matrice est rangée sous la forme d’un tableau mono dimensionné nécessitant un tableau de pointeur d’entrée de colonne ADIA pour explorer la matrice : l’entrée se fait par les termes diagonaux, le nombre de termes de la colonne est obtenu par différences de deux termes successifs : ADIA(i+1) - ADIA(i).
Si la matrice est non symétrique, mais à profil symétrique, il est nécessaire de stocker la i-ème ligne et la i-ème colonne. Classiquement, on les mets « bouts à bouts » et le nombre de termes de la colonne ou de la ligne est (ADIA(i+1) - ADIA(i))/2.
A1.4 Stockage par bloc
Les méthodes de stockage vues précédemment supposent implicitement que la matrice puisse résider en mémoire centrale, ce qui n’est pas toujours le cas.
D’où la notion de matrice stockée par bloc (ou segmentée sur disque).
Tous les stockages précédents peuvent être segmentés, mais nous n’exposerons que le cas de la matrice symétrique stockée profil.
Matrice profil stockée par bloc
Nous ne considérons ici que le cas des matrices symétriques, ce qui n’ôte rien à la généralité du propos.
Figure A1.4-a: Taille maximale d’un bloc : 20 éléments
Dans cet exemple, nous supposons les blocs de même taille, pour utiliser des blocs de taille variable, il faut introduire un tableau supplémentaire contenant la taille de chaque bloc.
Nous considérons également qu’une colonne ne peut appartenir qu’à un seul bloc : « nous ne coupons pas les colonnes ».
Pour gérer la matrice, il est toujours nécessaire de connaître l’adresse des termes diagonaux ; mais maintenant, cette adresse est relative au bloc d’appartenance de la colonne.
A ce tableau, il est nécessaire de joindre un tableau donnant les équations contenues dans un bloc.
Ce tableau dimensionné au nombre de bloc plus 1 contient la dernière équation du bloc.
Variations sur l’algorithme de GAUSS
Comme nous l’avons vu avec la variante de CROUT, il existe plusieurs implémentations de l’algorithme de GAUSS : elle consiste à effectuer les calculs des coefficients dans un ordre différent.
Schématiquement, on considère qu’il y a trois boucles imbriquées :
boucle i sur les lignes,
boucle j sur les colonnes,
boucle k sur les étapes.
L’algorithme standard se caractérise par la séquence kij, mais il existe 5 autres permutations des indices qui donnent lieu à autant de variantes (ou d’algorithmes).
l’algorithme de CROUT se caractérise par la séquence jki,
l’algorithme correspondant à la séquence ikj, qui travaille par ligne, est connu sous le nom d’algorithme de « Doolitle ».
Donnons ici une représentation graphique tirée de [bib10].
Bibliographie#
P.D. CROUT A short method for evaluating determinants and solving systems of linear equations with real or complex coefficients.-AIEE Trans Vol 60, 1941, pp 1235 - 1240
CUTHILL & J. Mc KEE « Reducing the band with of sparse symmetric matrices » Proc 24th Nat Conf Assoc Comput Mach ACM Publ (1969)
CUTHILL « Several strategies for reducing the bandwith of the matrices » dans Sparse Matrices and their applications - D.J. ROSE & R.A. WILLOUGHBY Editeurs, Plenum Press, New York (1972) pp 157 - 166
Von FUCHS, J.R. ROY, E. SCHREM Hypermatrix solution of large sets of symmetric positive definite linear equations - Computer methods in Applied Mechanics and Engineering (1972)
GEORGE, D.R. Mc INTYRE « On the application of the minimum degree algorithm to finite element systems » SIAM J. Num. An. Vol 15, (1978) pp90 - 112
G.H. GOLUB et C.F. Van LOAN Matrix computations. Johns Hopkins University Press - Baltimore (1983)
IRONS Roundoff criteria in direct stiffness solutions - AIAA Journal 6 n° 7 pp 1308 -1312 (1968)
IRONS A frontal solution program for finite elements analysis - Int Journal Num. Meth. Eng.,2, 1970
O’LEARY & STEWART, Computing the eigenvalues and eigenvectors of symmetric arrowhead matrices, J. of comp physics 90, 497-505 (1990)
J.M. ORTEGA « Introduction to parallel and vector solution of linear systems Plenum Press (1988)
RADICATI di BROZOLO, M. VITALETTI Sparse matrix-vector product and storage representations on the IBM 3090 with vector facility - IBM-ECSEC Report G513 - 4098 (Rome) July 1986
J.K. REID A note on the stability of gaussian elimination. J. Int Maths Applies, 1971, 8 pp 374 - 375
J.H. WILKINSON Rounding Errors in Algebraic. Processes Her majesty’s stationery office (1963)
J.H. WILKINSON The algebraic eigenvalue problem Clarendon Press Oxford (1965)
ROSE « Une méthode multifrontale pour la résolution directe de systèmes linéaires » Note EDF - DER HI-76/93/008 (1993)
SELIGMANN « Algorithmes de résolution pour le problème généralisé aux valeurs propres » [R5.01.01] - Note EDF - DER HI-75/7815 (1992)
SELIGMANN, R. MICHEL « Algorithmes de résolution pour le problème quadratique aux valeurs propres » [R5.01.02] - Note EDF - DER HI-75/7816 (1992)
Description des versions du document#
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
2.3 |
D.SELIGMANNEDF-R&D/AMA |
Texte initial |