Discussion :

[Aleph69]

Bonjour,

Nan mais de toute façon les interactions à N corps nécessite d'être traiter en parallèle si possible, car le temps de calcul est proportionnel en O(N^2).

C'est la complexité d'un bête produit scalaire...

la décomposition LU me donnera juste le résultat qui m'intéresse, et en difficulté O(N^2) également (au lieu de O(N^3) du pivot de Gauss), je suis vraiment débile de ne pas avoir vu ça avant

La décomposition LU se fait par la méthode du pivot de Gauss et nécessite N^3/3 opérations élémentaires.
Ensuite, la résolution des deux systèmes triangulaires pour calculer la Nième colonne de ton inverse se fait en O(N^2).
Ta matrice étant générale, il faut appliquer une stratégie de pivotage partiel pour assurer l'existence de la factorisation.
Si tu veux des résultats précis, il est conseillé d'appliquer a priori une méthode de préconditionnement diagonal et a posteriori un procédé de raffinement itératif sur la solution calculée.

La rapidité de calcul est cruciale dans mon cas.

Si, par hasard, les coefficients d'une matrice NxN ne sont pas modifiés dans une matrice N'xN' lorsqu'on passe de N à N' (N'>N), alors tu peux gagner énormément de temps en adoptant un format de stockage par bloc pour tes matrices et en calculant tes factorisations LU par bloc.
Note également que si tu n'as réellement besoin que de calculer la dernière colonne des inverses de tes matrices, il n'est aucunement nécessaire de stocker tous les coefficients du facteur L mais seulement le dernier.

[Asohan]

La complexité d'un bête produit scalaire ?? Tu plaisantes ? Imagine que tu aies 1000 FLOP (par couple de particules) de calcul et que ta simulation nécessite 2500 particules, pour environ 10000 pas de temps ? ça fait 2500^2 / 2 * 1000 * 10000 = 31250 GFLOP nécéssaires pour achever la simulation. Sur un ordi pas trop mal (j'ai un Core i5 655k) on peut monter à environ 10 GFLOPS en peak pour du monothread en double (c'est mon cas pour l'instant).

Donc au final il me faudrait : 31250 / 10 soit environ 1 heure. C'est le cas ideal. Sauf qu'en realité : j'ai énormement d'appel a des tableaux de données avec interpolation linéaire, donc ce serait plutôt 150000 FLOPS par couple de particule et plutôt 100000 pas de temps (sachant que mon soucis principal est la précision de mon système et la conservation systématique de l'énergie).

Tout ça est assez technique et il faudrait que j'explique la physique que j'utilise pour vous montrer la complexité du problème (les intéractions en physique subatomique et dans Crysis c'est pas vraiment le même niveau de complexité ).

Sinon pour la décomposition LU, je viens de revérifier ça ....effectivement j'ai un doute sur la procédure pour ne calculer que l'inverse de la dernière colonne. Je vais sans doute garder mon pivot de Gauss. Bien sur ce n'est qu'un calcul local et temporaire, je ne stocke que les valeurs propres calculées (qui change radicalement à chaque pas de temps ).

Les coefficient de la matrice sont crées de manière à avoir la certitude que l'inversion de matrice est possible, mais j'ai tout de même des sécurités.

J'ai déjà pensé à de nombreuse optimisations (stocker uniquement ce qui est nécessaire, réduire les formules à leur minimum d'opération, redimensionner ou calculer uniquement si nécessaire, ...) mais malgré tout, ce sont de longues simulations.

Par exemple, par rapport à de "simples" programmes de calcul de structure moléculaires, mes calculs d'interactions sont 4 à 5 fois plus lourds. Et eux ils utilisent déjà le parallélisme des GPU alors dans mon cas ... j'ai plutôt intérêt trouver vite une solution. D'autant plus qu'il me faut de la statistique, et donc : des milliers de simulations

Merci à toi Aleph69 pour tes conseils et ton intéressement au sujet. Je continue de creuser la piste de l'inversion de matrice pour voir un peu les solutions. Mais le pivot de gauss sera très facilement parallélisable sur GPU, et je pourrait utiliser à coté le CPU pour les tris de valeurs donc : vive Gauss ! (pour le moment )

[Aleph69]

Bonjour,

La complexité d'un bête produit scalaire ?? Tu plaisantes ?

J'aurais bien aimé mais non... j'ai juste écrit une grosse c.......
C'est la complexité d'un produit matrice-vecteur.
Désolé, je met ca sur le compte de la fatigue.

Donc au final il me faudrait : 31250 / 10 soit environ 1 heure. C'est le cas ideal. Sauf qu'en realité : j'ai énormement d'appel a des tableaux de données avec interpolation linéaire, donc ce serait plutôt 150000 FLOPS par couple de particule et plutôt 100000 pas de temps (sachant que mon soucis principal est la précision de mon système et la conservation systématique de l'énergie).

Oui, je connais bien le problème, je me retrouve souvent dans la même situation avec mes simulations par éléments finis.
Pour réduire le nombre de pas de temps, tu peux tenter un schéma pas à pas adaptatif.

Sinon pour la décomposition LU, je viens de revérifier ça ....effectivement j'ai un doute sur la procédure pour ne calculer que l'inverse de la dernière colonne. Je vais sans doute garder mon pivot de Gauss. Bien sur ce n'est qu'un calcul local et temporaire, je ne stocke que les valeurs propres calculées (qui change radicalement à chaque pas de temps ).

Si tu veux uniquement la dernière colonne de l'inverse, il faut vraiment passer par la résolution du système linéaire que je t'ai donnée.
Il n'y a pas de doute à avoir sur la procédure.
Cela découle de l'identité A*inv(A)=I, avec A une matrice, inv(A) sont inverse et I la matrice identité.
C'est mieux que d'inverser A : tu gagnes en performances, en stockage et surtout en précision de calcul.

Les coefficient de la matrice sont crées de manière à avoir la certitude que l'inversion de matrice est possible, mais j'ai tout de même des sécurités.

Sur le papier, je veux bien croire que ta matrice est inversible.
Le problème survient en précision finie.
Ta matrice peut être très mal conditionnée.
Or, le conditionnement correspond géométriquement à l'inverse de la distance qui sépare ta matrice de l'espace des matrices de même taille qui sont singulières.
La conséquence de cela est que ta matrice peut être singulière relativement à la précision de tes calculs.

D'autant plus qu'il me faut de la statistique, et donc : des milliers de simulations

Peut-être qu'il est possible de savoir le nombre d'échantillons qu'il te faut si tu sais à peu près l'intervalle de confiance que tu veux sur tes résultats.
Tu fais du Monte-Carlo?

[Asohan]

Mon pas de temps est adaptatif, bien sur (heureusement sinon bonjour les dégats ).

C'est bien ce qui me semblait pour l'inversion de matrice par LU. Je vais creuser ça à tout hasard, et voir si d'office c'est facilement parallélisable pour la suite.

Sinon je sais très bien que c'est la précision qui joue lors de l'inversion, mais dans mon cas le problème n'apparait jamais (par expérience). Je n'ai que des matrice du style : plein de (1 + epsilon[i][j]) partout, mais avec des epsilon qui permettent de rester loin des limites de précision en double float (10^-15 décimale je crois).

Quand je disais des millier, je me suis mal exprimé : pour 1 condition initiale donnée, seule une dizaine de simulations suffisent (en effet avec 2500 particules, les effets globaux sont gommés). Le problème c'est que j'ai 8 paramètres initiaux, et que je dois tous les tester (différentes combinaisons).
Ça fait un paquet de simulations totalement différentes au final.

J'utilise Monte-Carlo, mais pas ici. Ici il s'agit d'étudier les effets microscopiques sur la dynamique des formules appliquées, et au final de comparer aussi les résultats globaux avec ceux des autres méthodes de simulations (qui utilisent en général du tirage Monte-Carlo) qui elles n'ont pas une approche microscopique.

[Aleph69]

Bonsoir,

C'est bien ce qui me semblait pour l'inversion de matrice par LU. Je vais creuser ça à tout hasard, et voir si d'office c'est facilement parallélisable pour la suite.

La parallélisation des factorisations LU est un domaine de recherche très actif.
Ceci dit, cela concerne plutôt les matrices creuses.
Pour les matrices denses, ça ne me semble pas particulièrement compliqué dans la mesure où tu imbriques correctement tes trois boucles.
La manière classique (k-i-j avec k l'étape, i la ligne, j la colonne) est utilisée dans les approches multifrontales : c'est utile pour les matrices creuses ou assemblées par éléments finis.
La version left-looking (i-k-j ou i-j-k je ne sais plus) te permet de calculer tes facteurs par ligne.
La version right-looking (j-k-i ou j-i-k parce que je ne sais toujours pas plus) calcule par colonne.

Il faut choisir la bonne selon tes besoins.
Tu regarder ce petit document qui fait un petit résumé des techniques mais il faut chercher d'autres sources.

Le problème c'est que j'ai 8 paramètres initiaux, et que je dois tous les tester (différentes combinaisons).
Ça fait un paquet de simulations totalement différentes au final.

Il est sûrement possible de faire un plan d'expérience pour réduire le nombre de simulations.

J'utilise Monte-Carlo, mais pas ici. Ici il s'agit d'étudier les effets microscopiques sur la dynamique des formules appliquées, et au final de comparer aussi les résultats globaux avec ceux des autres méthodes de simulations (qui utilisent en général du tirage Monte-Carlo) qui elles n'ont pas une approche microscopique.

Ton projet a l'air complet et vraiment très intéressant, je t'envie!
Bon courage et bonne continuation!

[Asohan]

Merci infiniment pour ton aide, c'est vrai qu'a force de faire de la physique, on oublie vite les maths de base et les langages informatiques.

J'espère bien que mon projet est intéressant A vrai dire je suis en doctorat de physique et c'est mon sujet de thèse . Mon sujet de base est l'étude d'un modèle physique + des calculs nouveaux. Finalement, j'ai aussi eu besoin de rédiger en LaTeX, de programmer un code de simulation complet (c'est moi qui ai écrit et optimisé chaque routine depuis deux ans) et un autre programme pour le calcul intensif d'intégrales multiples avec minimisation. Bref, des compétence vraiment hors de mon domaine. Mais c'est très intéressant aussi .

Tous mes problèmes à la base sont dus au fait que les formules que j'applique donne des variation très brutales des trajectoires avec en plus : des collisions, des fusions (N--), des fissions (N++), ...(j'ai des particules élémentaires et aussi des composées).

Ah oui et pour finir j'ai fait cette dynamique d'une manière "relativiste" (c'est à dire des particules dont la vitesse est proche de celle de la lumière). Pour faire simple : on utilise de nouvelles équations de mouvement des particules, basées sur la résolution d'un système d'équations linéaires (dues à des contraintes imposées sur le système).

Bref : que du bonheur tout ça sauf que je suis le seul de mon labo en tant que physicien à pouvoir également implémenter ses idée de manière appliquée. Et c'est ça qui est bien : montrer ce que l'on fait en vrai (des courbes de maths, c'est bien, mais une petite vidéo 3D OpenGL de son travail, c'est encore mieux pour vulgariser son travail ). A oui j'ai fait aussi du Qt et (Qt)OpenGL afin de faire un logiciel interactif et facile d'utilisation (et bientôt OpenCL !)

Merci encore de m'avoir aidé. C'est vrai que ce qu'il manque pendant la thèse : c'est des collègues thésards informaticiens et matheux afin de s'entraider .

Si jamais j'ai l'autorisation du directeur, je mettrais mon code complet sur le site de manière libre. Ça pourra servir de base à d'autres travaux. Et puis mon code a très peu (si ce n'est aucune) dépendance à d'autre librairies (hormis Qt dans la version interactive du code).

[Aleph69]

Joli CV!
C'est ton futur recruteur qui va être content!

Si jamais j'ai l'autorisation du directeur, je mettrais mon code complet sur le site de manière libre. Ça pourra servir de base à d'autres travaux. Et puis mon code a très peu (si ce n'est aucune) dépendance à d'autre librairies (hormis Qt dans la version interactive du code).

Je suis preneur!
Quand tu auras rédigé et soutenu ta thèse, je serai ravi d'en avoir une copie.

[Asohan]

Ah non, pas de recruteur dans mon domaine. Seul 1 docteur sur 3 est recruté dans la recherche, 1/3 dans le privé (c-à-d qu'il faut accepter une "reconversion") et 1/3 ne trouve rien et recommence d'autre étude ("Game Over ... try again ? Insert coin !"). Vive la France !

D'autant plus que les salaires des chercheurs en France sont les seuls au monde à être inférieurs (en moyenne sur 40 ans de travail) à ceux des ingénieurs ! Bah oui dans les autres pays : un docteur est un docteur. Que ce soit en médecine, en informatique, en chimie ou autre. En France, un médecin touche quasiment le double de ses homologues scientifiques et un ingénieur environ 10 à 50% de plus.

Alors pourquoi se casser la tête à résoudre ces problèmes de C++ quand on est pas informaticien ... on va dire que c'est parce qu'on est passionnés . En plus dans la recherche, on il n'y a pas vraiment de limite à notre activité (tant qu'on fait ce pour quoi on est payé, on peut en plus faire autre chose à coté) et on a pas le droit (en physique en tout cas) de déposer de brevet en tant que chercheur sur "ce que l'on fait"

Alors : salaire moyen + recrutement moyen + aucun moyen de rémunérer ses compétences supplémentaires = je vous laisse deviner ma pensée

En tout cas je suis content que vous m'ayez aidé. Le savoir ne doit pas avoir de limite, et il ne doit pas se payer. Je pense que même si je n'ai pas l'accord du directeur, je vais poster mon code quand il sera fini. Si ça peut éviter au autre de refaire les même erreurs que moi ... (je mets un maximum de commentaires dans le code pour expliquer la méthodo)

Par contre ce code c'est : 11 fichiers pour environ 6500 lignes de code (tout compris : commentaires, espaces et fichiers .h) J'espère que ça passera sur ce site

[Aleph69]

Bonsoir,

Ah non, pas de recruteur dans mon domaine. Seul 1 docteur sur 3 est recruté dans la recherche, 1/3 dans le privé (c-à-d qu'il faut accepter une "reconversion") et 1/3 ne trouve rien et recommence d'autre étude ("Game Over ... try again ? Insert coin !"). Vive la France !

C'est pas mal, ça fait déjà 2/3!

D'autant plus que les salaires des chercheurs en France sont les seuls au monde à être inférieurs (en moyenne sur 40 ans de travail) à ceux des ingénieurs ! Bah oui dans les autres pays : un docteur est un docteur. Que ce soit en médecine, en informatique, en chimie ou autre. En France, un médecin touche quasiment le double de ses homologues scientifiques et un ingénieur environ 10 à 50% de plus.

En fait ça dépend de ce que tu appelles un chercheur.
Les situations d'un MCF, d'un prof, d'un CR ou DR CNRS, etc, sont complètement différentes.

Et puis, si tu veux voir le verre à moitié plein, tu regardes du côté de l'Italie (peu de rectrutement dans la recherche publique) ou de l'angleterre (nécessité de faire vivre soi-même son labo avec des contrats).

Si tu considères la Belgique, le Luxembourg ou la Suisse, le verre commence à se vider plus sérieusement...

Aux états-unis, ton diplôme de docteur sera mieux reconnu mais il faut vouloir s'expatrier.
Et là encore, tu vas devoir travailler avec des financements privés.
Ce n'est pas nécessairement une mauvaise chose sauf lorsque tu souhaites faire des sciences non appliquées (mathématiques pures, physique théorique, etc).

Il ne faut pas être aussi pessimiste.
Les profils de "numériciens" comme le tien sont très recherchés.
Tu trouveras des travaux similaires aux tiens au CEA, à l'INRIA, à EdF, au Cermics, ... et bien sûr dans le privé!

Je pense que même si je n'ai pas l'accord du directeur, je vais poster mon code quand il sera fini.

J'espère ne pas te l'apprendre mais en fait ton code ne t'appartient pas : ton labo en est propriétaire.

[Asohan]

C'est si recherché que ça les "numériciens" ? Je croyais que tous les informaticiens étaient déjà, de base, des "numériciens".

Pour ce qui est du boulot, c'est vrai que je vais sans doute me débrouiller, comme on dit : wait and see

Pour ce qui est du code, c'est un peu plus compliqué je crois En fait je suis hébergé par le labo, donc mon rapport de thèse et mes résultats appartiennent au labo, c'est clair. En revanche, les moyens mis en œuvre pour arriver a ces résultats (code source, templates LaTeX, ...) sont a mon entière discrétion si je me trompe pas.

Et puis de toute façon on doit publier son code entre labo pour que les autres vérifie les résultats alors ... dans tout ça je suppose qu'il y a des fuites