Inversion matrice bande (et non pas "résolution")

Bonjour à toutes et à tous,

Je souhaite inverser une matrice bande de très grande dimension (du style (m*m)**2 où m peut valoir 10^2) ce qui prend avec un algorithme classique du style décomposition LU un temps phénoménal, ce qui est problématique car inclus dans une boucle itérative. La matrice est quelconque (donc non symétrique)

Avant tout, je tiens à préciser, puisque c'est visiblement la source de beaucoup de confusions sur le forum, que je souhaite bien INVERSER la matrice et non résoudre directement un système (l'inversion fait partie d'un problème dont les détails sont je pense inutiles)

L'ennuie est que les libraires que j'ai pu trouver dans LAPACK par exemple ne proposent que des RESOLUTIONS de systèmes et non des inversions matricielles. Je cherche donc un programme capable d'inverser une une telle matrice bande.

Merci d'avance pour vos propositions.

Mathieu

Salut!

Citation:

---

L'ennuie est que les libraires que j'ai pu trouver dans LAPACK par exemple ne proposent que des RESOLUTIONS de systèmes et non des inversions matricielles. Je cherche donc un programme capable d'inverser une une telle matrice bande.

---

Ce n'est pas un oubli de la part des auteurs de LAPACK, mais une décision inévitable: l'inverse d'une matrice bande est une matrice pleine, d'où le besoin d'une place monstrueuse, sans commune mesure avec celle nécessitée par la matrice originale.

Citation:

---

l'inversion fait partie d'un problème dont les détails sont je pense inutiles

---

Là, j'ai des doutes!

Jean-Marc Blanc

Merci pour votre réponse. Voici donc plus de détails :

J'ai besoin de résoudre un système A*x=b où :
- x et b sont de dimensions (2*(m-2)*(n-2),1)
- A une matrice de dimensions (2*(m-2)*(n-2),2*(m-2)*(n-2))

J'arrive à résoudre le problème à l'aide d'une décomposition LU mais la durée n'est raisonnable que pour un nombre de points de maillage (m,n) très faible.
Or la matrice A se compose de 4 blocs : 2 pentadiagonaux et 2 diagonaux (c'est une jacobienne que j'utilise dans un algorithme de Newton-Raphson pour résoudre un problème non linéaire).
Je souhaitais donc inverser la matrice A de départ par blocs, et c'est là qu'intervient l'inversion des matrices pentadiagonales.

PS: Bien entendu si vous avez d'autres idées pour résoudre mon système A*x=b, je suis preneur puisque l'inversion que je souhaitais faire n'est pas une fin en soit. J'étais juste un peu à court d'idées ! ;)

Salut!

Citation:

---

J'ai besoin de résoudre un système A*x=b

---

Citation:

---

c'est une jacobienne que j'utilise dans un algorithme de Newton-Raphson pour résoudre un problème non linéaire

---

Le calcul de la matrice inverse n'est pas une bonne méthode pour résoudre un système linéaire. Pour plus de détails, regarde dans le tutoriel que j'ai écrit: "Résolution des systèmes linéaires". Tu trouveras aussi d'excellents codes dans la bibliothèque LAPack, disponible sur www.netlib.org

Jean-Marc Blanc

je me suis peut-être mal exprimé :

En fait, j'arrive déjà à résoudre mon système général A*x=b, justement en utilisant une librairie LAPack qui est la suivante (puisque que ma matrice est vraiment quelconque) :

file dgesv.f dgesv.f plus dependencies
prec double
for Solves a general system of linear equations AX=B.
gams d2a1

Cependant, la taille du système étant trop importante pour avoir une convergence en un temps raisonnable, je cherche un moyen plus rapide. D'où mon essai d'inversion par bloc, peu fructueux je vous l'accorde.
Je suis donc toujours dans l'impasse... :/

PS: En tout cas merci pour les tutoriels (que j'ai déjà lu) et surtout "le language fortran" qui ne quitte pas mon bureau ! ;)

Salut!

Citation:

---

file dgesv.f dgesv.f plus dependencies

---

Si tu as une matrice "bande", pourquoi n'utilises-tu pas plutôt DGBSV ?
Jean-Marc Blanc

Bonjour,

Citation:

---

la matrice A se compose de 4 blocs : 2 pentadiagonaux et 2 diagonaux

---

Etant donné la structure très creuse de ta matrice, il est clair qu'il faut tirer parti de cette structure et éviter de la manipuler ou l'inverser d'un coup.
Puisqu'au final il s'agit bien de résoudre des systèmes linéaires, je pense qu'une approche du type factorisation de matrice, à la Faddeev & Faddeva serait payante: c'est en gros du LU par bloc, mais sur des blocs particuliers (ceux contenant les bandes, pour lesquels la décomposition est bien moins couteuse que pour une matrice pleine, puisque la décompistion LU d'une matrice bande donne des matrices L et U de même largeur de bande).
L'algo. est décrit (entre autres) dans le bouquin de J.P. Boyd (appendix B5) que tu peux télécharger (gratuitement!) sur son site.

Bon courage.

Ehouarn

Jean-Marc, en fait c'est un des blocs composant la matrice qui est "bande" et non pas la matrice complète, qui est elle quelconque.

Ehouarn , je vais regarder cet algorithme qui, d'après la description que tu m'as faite, semble correspondre à l'idée que j'avais en tête.

Merci à tous les deux.

Salut!
La difficulté à te répondre provient de ce que tu n'as pas posé ta question de manière judicieuse, c'est-à-dire en commençant par le début. Si tu es étudiant, c'est un prof. qui t'a posé un problème, si tu es dans l'industrie, c'est ton patron ou un client. Dans un cas comme dans l'autre, le problème initial n'est pas de résoudre un système non linéaire mais quelque-chose d'autre que tu ne nous as pas révélé. Alors dis-nous d'abord quel est le problème qui t'a été posé, puis ce que tu as essayé de faire. Ainsi, on pourra te dépanner.
Jean-Marc Blanc

Bonjour,

Pour une matrice creuse de grande taille j'essaierais plutôt une méthode itérative ne faisant que des produits matrice x vecteur.
Comme en plus il s'agit en fait de résoudre un système non-linéaire par la méthode de Newton-Raphson, il n'est peut-être même pas nécessaire de résoudre parfaitement chaque système linéaire et donc de faire beaucoup d'itérations.
En plus comme il y a des blocs bien identifiés, il faut voir si une parallélisation des produits matrice x vecteur avec OpenMP ne donnerait pas de bons résultats.

JBF

Salut à tous,

- Ehouarn, j'ai mis en place la méthode que tu m'as conseillée et qui est, dans mon cas, vraiment avantageuse, car je pense avoir divisé par 10 le temps de calcul. Même si ce dernier n'est pas encore suffisamment réduit, il y a déjà un grand pas en avant.

- Jean-Marc, je n'ai en effet pas décrit le problème dans son ensemble : Il s'agit de résoudre simultanément deux équations non-linéaires : l'équation de Reynolds pour les écoulements en film mince et une loi d'état de comportement du fluide en question, avec pour inconnu la pression et la masse volumique.

- JB, merci pour ton conseil, je vais regarder plus en détails ce que cela pourrait donner sur mon problème.

Merci à tous!

Mathieu