Publicité
+ Répondre à la discussion
Page 1 sur 2 12 DernièreDernière
Affichage des résultats 1 à 20 sur 21
  1. #1
    Invité de passage
    Homme Profil pro
    Ingénieur développement logiciels scientifiques
    Inscrit en
    décembre 2011
    Messages
    13
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France

    Informations professionnelles :
    Activité : Ingénieur développement logiciels scientifiques

    Informations forums :
    Inscription : décembre 2011
    Messages : 13
    Points : 2
    Points
    2

    Par défaut résolution de systéme linéaire

    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)
    Images attachées Images attachées

  2. #2
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 191
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 191
    Points : 1 546
    Points
    1 546

    Par défaut

    Salut,

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

  3. #3
    Invité de passage
    Homme Profil pro
    Ingénieur développement logiciels scientifiques
    Inscrit en
    décembre 2011
    Messages
    13
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France

    Informations professionnelles :
    Activité : Ingénieur développement logiciels scientifiques

    Informations forums :
    Inscription : décembre 2011
    Messages : 13
    Points : 2
    Points
    2

    Par défaut

    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
    Images attachées Images attachées
    Fichiers attachés Fichiers attachés

  4. #4
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 191
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 191
    Points : 1 546
    Points
    1 546

    Par défaut

    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

  5. #5
    Invité de passage
    Homme Profil pro
    Ingénieur développement logiciels scientifiques
    Inscrit en
    décembre 2011
    Messages
    13
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France

    Informations professionnelles :
    Activité : Ingénieur développement logiciels scientifiques

    Informations forums :
    Inscription : décembre 2011
    Messages : 13
    Points : 2
    Points
    2

    Par défaut

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

  6. #6
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 191
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 191
    Points : 1 546
    Points
    1 546

    Par défaut

    alors envoie en plusieurs parties

  7. #7
    Invité de passage
    Homme Profil pro
    Ingénieur développement logiciels scientifiques
    Inscrit en
    décembre 2011
    Messages
    13
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France

    Informations professionnelles :
    Activité : Ingénieur développement logiciels scientifiques

    Informations forums :
    Inscription : décembre 2011
    Messages : 13
    Points : 2
    Points
    2

    Par défaut

    Voilà la matrice obtenue pour N=36.
    Fichiers attachés Fichiers attachés

  8. #8
    Invité de passage
    Homme Profil pro
    Ingénieur développement logiciels scientifiques
    Inscrit en
    décembre 2011
    Messages
    13
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France

    Informations professionnelles :
    Activité : Ingénieur développement logiciels scientifiques

    Informations forums :
    Inscription : décembre 2011
    Messages : 13
    Points : 2
    Points
    2

    Par défaut

    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.

  9. #9
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 191
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 191
    Points : 1 546
    Points
    1 546

    Par défaut

    Bonjour,

    Citation Envoyé par odessos Voir le message
    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.

  10. #10
    Invité de passage
    Homme Profil pro
    Ingénieur développement logiciels scientifiques
    Inscrit en
    décembre 2011
    Messages
    13
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France

    Informations professionnelles :
    Activité : Ingénieur développement logiciels scientifiques

    Informations forums :
    Inscription : décembre 2011
    Messages : 13
    Points : 2
    Points
    2

    Par défaut

    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
    Fichiers attachés Fichiers attachés

  11. #11
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 191
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 191
    Points : 1 546
    Points
    1 546

    Par défaut

    Citation Envoyé par odessos Voir le message
    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.

  12. #12
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 191
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 191
    Points : 1 546
    Points
    1 546

    Par défaut

    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.

  13. #13
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 191
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 191
    Points : 1 546
    Points
    1 546

    Par défaut

    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.

  14. #14
    Invité de passage
    Homme Profil pro
    Ingénieur développement logiciels scientifiques
    Inscrit en
    décembre 2011
    Messages
    13
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France

    Informations professionnelles :
    Activité : Ingénieur développement logiciels scientifiques

    Informations forums :
    Inscription : décembre 2011
    Messages : 13
    Points : 2
    Points
    2

    Par défaut

    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.

  15. #15
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 191
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 191
    Points : 1 546
    Points
    1 546

    Par défaut

    Citation Envoyé par odessos Voir le message
    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.

    Citation Envoyé par odessos Voir le message
    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.


    Citation Envoyé par odessos Voir le message
    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.

  16. #16
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 191
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 191
    Points : 1 546
    Points
    1 546

    Par défaut

    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.
    Fichiers attachés Fichiers attachés

  17. #17
    Invité de passage
    Homme Profil pro
    Ingénieur développement logiciels scientifiques
    Inscrit en
    décembre 2011
    Messages
    13
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France

    Informations professionnelles :
    Activité : Ingénieur développement logiciels scientifiques

    Informations forums :
    Inscription : décembre 2011
    Messages : 13
    Points : 2
    Points
    2

    Par défaut

    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

  18. #18
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 191
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 191
    Points : 1 546
    Points
    1 546

    Par défaut

    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.

  19. #19
    Invité de passage
    Homme Profil pro
    Ingénieur développement logiciels scientifiques
    Inscrit en
    décembre 2011
    Messages
    13
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France

    Informations professionnelles :
    Activité : Ingénieur développement logiciels scientifiques

    Informations forums :
    Inscription : décembre 2011
    Messages : 13
    Points : 2
    Points
    2

    Par défaut

    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

  20. #20
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 191
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 191
    Points : 1 546
    Points
    1 546

    Par défaut

    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.

+ Répondre à la discussion
Cette discussion est résolue.
Page 1 sur 2 12 DernièreDernière

Liens sociaux

Règles de messages

  • Vous ne pouvez pas créer de nouvelles discussions
  • Vous ne pouvez pas envoyer des réponses
  • Vous ne pouvez pas envoyer des pièces jointes
  • Vous ne pouvez pas modifier vos messages
  •