résolution de systéme linéaire [Résolu]

[odessos]

Bonjour à tous,

je souhaiterai résoudre un systéme linéaire de type Ax = b pour toute matrice de taille
N² x N² ayant une structure semblable à celle de l'image mise en lien (dans laquelle N=36).

Cependant, la résolution ne pose pas de problème jusqu'à N=35 inclus, mais donne des résultats aberrants pour les valeurs plus grandes

Pour quelles raisons pensez-vous qu'il existe-t-il un tel seuil?
Est-ce plutôt dû à un problème de mémoire ou à un problème de conditionnement de A?

Merci d'avance

p.s. je fais la résolution avec pardiso (décomposition LU)

[Aleph69]

Salut,

impossible de te répondre sans les valeurs de ta matrice. Peux-tu les envoyer?

[odessos]

J'ai essayé d'envoyer la matrice N=36 à partir de laquelle le résultat est faux mais elle est trop volumineuse.
Je vous transmet la matrice N=10 et sa structure, les ordres de grandeurs sont les même quelquesoit la taille de la matrice.

p.s. j'ai mis le fichier en .cpp pour que le site accepte le fichier
la premiére ligne indique le nombre d'éléments non-nuls suivi de la dimension

[Aleph69]

Salut,

la matrice que tu m'as envoyée est inversible et parfaitement bien conditionnée. Ses valeurs propres sont strictement positives et contenues dans l'intervalle [9,10]. En revanche, elle admet des valeurs propres multiples : il n'y a que 31 valeurs propres distinctes sur 100. Cela dit, je ne vois pas ce qui empêcherait pardiso de factoriser ta matrice avec un pivotage partiel.

La structure ne correspond pas à l'image que tu as envoyée. Chez moi, on voit grossièrement un bloc diagonal alors que dans ton image il est antédiagonal. Est-ce normal?

Il faut envoyer la matrice pour N=36 sinon je ne peux rien faire. Ta matrice étant creuse, écris ta matrice avec l'algo suivant :

Code :

1
2
3
4
5
6
7
for i=1:N*N
    for j=1:N*N
        if mat(i,j)~=0
             écrire sur une même ligne "i j mat(i,j)" avec des espaces pour séparer
        end
    end
end

[odessos]

Non, même comme ça il ne passe pas...

[Aleph69]

alors envoie en plusieurs parties

[odessos]

Voilà la matrice obtenue pour N=36.

[odessos]

la matrice que tu m'as envoyée est inversible et parfaitement bien conditionnée

Qu'utilises-tu comme bibliothèque/logiciel pour connaitre ce genre de chose?

La structure ne correspond pas à l'image que tu as envoyée. Chez moi, on voit grossièrement un bloc diagonal alors que dans ton image il est antédiagonal. Est-ce normal?

Effectivement... Je vais revoir mon algo d'affichage.

[Aleph69]

Bonjour,

Envoyé par odessos

Qu'utilises-tu comme bibliothèque/logiciel pour connaitre ce genre de chose?

Matlab mais ce n'est pas le seul choix possible.

Ta dernière matrice est également inversible et parfaitement bien conditionnée (4.9819). J'ai résolu Ax=b pour x=(1,1,1,......,1) par factorisation LU et j'obtiens une erreur relative sur la solution exacte de l'ordre de la précision machine (9.8524e-016). En revanche, il y a une chose bizarre : les valeurs propres de ta matrice sont maintenant dans le plan complexe, ce qui n'était pas le cas pour la première matrice que tu as envoyé.

Conclusion : la matrice étant inversible, il en existe une permutation admettant une factorisation LU. Pardiso doit pouvoir factoriser ta matrice.

Deux explications possibles :

1. tu n'utilises pas Pardiso avec les bonnes options;
2. tu considères aberrantes des solutions qui ne le sont pas.

Pour le 1., envoie le bout de code d'appel à Pardiso en précisant toutes les options passées en paramètres.

Pour le 2., envoie le second membre associé pour N=36, la solution que calcule Pardiso et éventuellement la solution exacte ou que tu penses être juste si disponible.

[odessos]

Deux explications possibles :
1. tu n'utilises pas Pardiso avec les bonnes options;
2. tu considères aberrantes des solutions qui ne le sont pas.

Alors je dois être dans le premier cas, car b-Ax me donne certaines valeurs à 10⁺¹²...

Voici les paramétres de Pardiso que j'utilise (la majorité ont été laissé par défaut) :

iparm[0] = 1; /* No solver default */
iparm[1] = 2; /* Fill-in reordering from METIS
iparm[2] = 8; /* Number of threads*/
iparm[3] = 0; /* CGS */
iparm[4] = 0; /* No user fill-in reducing permutation */
iparm[5] = 0; /* Write solution into x */
iparm[6] = 0; /* Not in use */
iparm[7] = 0; /* Max numbers of iterative refinement steps */
iparm[8] = 0; /* Not in use */
iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
iparm[11] = 0; /* Not in use */
iparm[12] = 0; /* Not in use */
iparm[13] = 0; /* Output: Number of perturbed pivots */
iparm[14] = 0; /* Not in use */
iparm[15] = 0; /* Not in use */
iparm[16] = 0; /* Not in use */
iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
iparm[18] = -1; /* Output: Mflops for LU factorization */
iparm[19] = 0; /* Output: Numbers of CG Iterations */
iparm[27] = 1; /* check the data structure */
iparm[31] = 1; /* iterative solver*/
maxfct = 1; /* Maximum number of numerical factorizations. */
mnum = 1; /* Which factorization to use. */

Je poste mon problème sur la mailing list de Pardiso, peut-être sauront-ils quel paramétre devrait être changé.

Je t'envoie quand même le systéme linéaire complet. La matrice A est un peu différente de celle que je t'ai envoyé, mais conserve la structure.

Merci beaucoup en tout cas

[Aleph69]

Envoyé par odessos

Alors je dois être dans le premier cas, car b-Ax me donne certaines valeurs à 10⁺¹²...

Cela ne remet pas en cause la pertinence de ton calcul. Si x=10^{13} et y =10^{13}+10^{12}, alors y-x=10^{12} mais |y-x|/|x|=10^{-1}. Ce qu'il faut mesurer c'est ||b-Ax||/||b||.

Je te tiens au courant pour le dernier système linéaire.

[Aleph69]

Salut,

le système linéaire que tu as envoyé ne pose aucun problème de résolution. J'obtiens une approximation à la précision machine. Je ne sais pas ce qu'est censé représenter le vecteur x, j'imagine la solution calculée par pardiso, mais il est très loin de la solution exacte du problème.

Bref, le problème ne vient pas de mauvaises propriétés du système linéaire mais bien d'une erreur d'algorithme.

[Aleph69]

Je te suggères une modification dans tes options : ajouter quelques itérations de raffinement itératif (<10) :
iparm[7] = 5; /* Max numbers of iterative refinement steps */
Cela permet d'améliorer grandement la précision de la solution calculée à un coût dérisoire (O(n²)) devant le coût de la factorisation (O(n^3)).

Cependant, je ne pense pas que le problème vienne de là.

EDIT : à tout hasard, essaye de résoudre sur un seul thread pour voir si le problème persiste. L'autre source potentiel d'erreur qui me vient à l'esprit est la mise en entrée d'une matrice erronée, au sens différente de ce que tu souhaites en fait passer en paramètre, à cause d'un changement de format de stockage creux mal implémenté par exemple (bogue sur la distribution des indices). Une mesure judicieuse de ce qui se passe serait de tracer l'évolution de l'erreur rétrograde ||b-Ax||/||b|| pour des valeurs croissantes de N.

[odessos]

le vecteur x est très loin de la solution exacte du problème

Tu veux dire de la solution que tu obtiens avec matlab?
Pourrais-tu m'envoyer cette matrice pour que j'ai une référence pour la suite?

L'autre source potentiel d'erreur qui me vient à l'esprit est la mise en entrée d'une matrice erronée

La matrice est calculée à l'intérieur du logiciel directement au format requis pour utiliser Pardiso... il ne devrait donc pas y avoir de problème

Je pense que l'erreur a de grandes chances de venir des paramètres de Pardiso. Je ne maîtrise pas suffisemment certaines notions pour être sur de ceux que j'ai choisi.

[Aleph69]

Envoyé par odessos

Tu veux dire de la solution que tu obtiens avec matlab?

et de la solution exacte puisque l'ordre d'approximation est de la précision machine.

Envoyé par odessos

Pourrais-tu m'envoyer cette matrice pour que j'ai une référence pour la suite?

Tu veux dire la solution calculée? Je te l'envoie dans la soirée.

Envoyé par odessos

Je pense que l'erreur a de grandes chances de venir des paramètres de Pardiso. Je ne maîtrise pas suffisemment certaines notions pour être sur de ceux que j'ai choisi.

Il y a un moyen simple de le vérifier. Tu construis le vecteur (1,1,.....,1) de la taille de ta matrice A, puis tu calcules b=A*x, puis tu résous Ay=b avec pardiso et enfin tu calcules l'erreur relative ||x-y||/||x||. Si tu as un problème avec l'utilisation de pardiso, la solution calculée y devrait être très éloignée de la solution x. A toi ensuite de trouver le bon jeu de paramètres pour que l'erreur relative soit petite. Si tu as des questions sur le sens de tel ou tel paramètre, je peux sûrement te renseigner.

[Aleph69]

Comme convenu, voici une approximation numérique de la solution exacte. L'erreur rétrograde associée, ||b-Ax||/||b||, pour la norme euclidienne, est de l'ordre de la précision machine : 4.9096e-016. A noter que le conditionnement estimé de A dans cette même norme est 1.0006; on peut difficilement espérer mieux.

EDIT : je ne sais pas d'où provient ta matrice mais elle est très bizarre. Elle est pratiquement semblable à la matrice identité. Précisément, ses valeurs propres sont toutes égales à 1,0002.

[odessos]

Je pense que je peux mettre cette conversation en "résolue". Je suis en train de voir avec l'équipe de Pardiso à quel niveau je commet l'erreur.

Pour info, la matrice A en question est de type (1 + B/c). Voyant que la résolution ne fonctionnait pas, j'ai simplifié le problème en mettant c >> Bij pour faire tendre A vers 1.
La résolution de ce systéme servira à terme à simuler la dispersion d'un polluant dans un paysage

Encore merci pour ton aide

[Aleph69]

Merci d'avoir pensé à passer en résolu.
J'espère que tu vas réussir à résoudre ton problème rapidement.
Bon courage pour la suite.

[odessos]

C'est bon, j'ai enfin résolu mon problème.

C'était bien un des paramétres de Pardiso qu'il fallait activer. Il faut mettre iparm[12] = 1 pour ceux qui aurait un problème semblable

iparm[12]=0 veut dire "Maximum weighted matching algorithm is switched-off". Je doit avouer ne pas trés bien comprendre ce que ça fait.

J'aurais jamais réussi à le trouver par moi-même : dans la doc que j'utilise, il est écris :
"iparm(12) This parameter is reserved for future use. Its value must be set to 0."

Content que ce soit résolu en tout cas.

Encore merci,
odessos

[Aleph69]

Salut,

je suis content d'apprendre que tu as réussi à résoudre ton problème. En revanche, je suis un peu surpris sur la solution proposée. L'algorithme "Maximum weighted matching" est une méthode de prétraitement permettant de permuter les lignes et colonnes de ta matrice de manière à ce que les coefficients de plus grande valeur absolue soient positionnés sur la diagonale. L'idée est de renforcer la (pseudo-)dominance diagonale. Bien que cela ait un effet sur les techniques de pivotage utilisées dans les méthodes directes, c'est surtout utile pour accélérer la convergence des méthodes itératives. Cela ne devrait pas être nécessaire de faire cette opération pour les matrices envoyées : elles sont beaucoup trop faciles à factoriser.

S'il s'agit de la seule manière de factoriser tes matrices avec pardiso, c'est que c'est un (très) mauvais solveur direct. Il n'est même pas capable de rivaliser avec matlab (solveur umfpack me semble-t-il). Je te conseille vivement de te tourner vers des solutions plus sûres telles que superlu ou mumps.

Je t'avoue que je suis vraiment très surpris car pardiso jouit d'une bonne réputation dans la communauté. Pour moi, le problème que tu rencontres relève plus du bogue qu'autre chose. Je n'ai aucun doute sur le fait que les gens qui t'ont répondu le savent mais ils se sont bien abstenus de te le dire.