réflexion sur le calcul matriciel

bonjour,

j'ai l'ambition de créer une bibliothèque en c++ pour faire du calcul matriciel. Dans l'idéal cette bibliothèque devra s'utiliser avec une syntaxe proche de celle de matlab.

Cette bibliothèque utilise les template.

Dans une première version, je me suis inspiré du conteneur vector pour la possibilité d'écrire les surcharges d'opérateur indépendamment de la dimension de la matrice créé.

Ainsi les vecteur<float> utilisent les mêmes opérateurs (en théorie) que les vecteurs<vecteurs<float>> ... peu importe la taille.

Un gros défaut de ma méthode actuelle : la mémoire.

En effet j'ai fait l'expérience de créer la matrice de transformé de fourrier. Ma matrice crée est de la taille 8192 au delà, le programme plante. cette matrice occupe 1Go de mémoire et met beaucoup plus de temps à être calculée que sous matlab.

ma question: quel est selon vous la meilleur façon de bien partir pour créer la bibliothèque?
quelle classe? s'appuyer sur des classes standard existante?
bref je vous propose d'y réfléchir avec moi.:mouarf:

Spa trivial le calcul matriciel en C++.

Deja vector<vector<T>> est tres tres mauvais comme conteneur de base car il alloue des bloc discontinue et donc empeche le cache et les pages d'etre accédé proprmeent. COmme d'hab, c'est dans les vieux pots etc... donc regarde du coté de l'allcoation façon NRC.

Ensuite, expression tempaltes ou alors tu va crouler sous le point des temporaires.

je cherche pas quelque chose de trivial :mouarf:.

Le coté mathématique de l'affaire ne me fait absolument pas peur. C'est plus la partie c++ j'ai découvert pas mal de chose coté c++ depuis que j'ai commencé. Et je réalise que je suis pas suffisamment au point pour faire quelque chose de propre seul.


NRC => ok petite recherche rapide dessus et je tombe sur pleins de doc. Avant de tout lire, ce qui me prendra pas mal de temps, est-ce que c'est allouer un grand tableau et simuler une matrice2D en calculant les offset des différentes lignes?

"Ensuite, expression tempaltes ou alors tu va crouler sous le point des temporaires." => j'ai découvert le principe des Template que récement, je suis pas sur d'en connaitre assez pour comprendre ce que tu veux me dire.


lien en rapport NRC (à lire demain :lol:)

http://www.nrbook.com/a/bookcpdf.php

quand je dis pas trivial c'ets au niveau implantation pas math. Pour bosser depuis 5 ans sur un projet simialire ... je t'assure que c'ets pas le truc qu'on design sur un coin de nappe le dimanche au pique-nique.

Les expressions tempaltes sont une technique qui permet de travailler à un grain plus gros que l'operateur et de proposer des optimizations global.

Citation:

---

NRC => ok petite recherche rapide dessus et je tombe sur pleins de doc. Avant de tout lire, ce qui me prendra pas mal de temps, est-ce que c'est allouer un grand tableau et simuler une matrice2D en calculant les offset des différentes lignes?

---

Oui, et ca a l'avantage de s'écrire récursivement

Dasn ton lien, le chapitre débutant page 15 doit contenir ce que tu cherches

Pour les template, Il y a un cours d'introduction à la méta-prog sur ce site où il est question à un moment de matrice (écrit par Laurent il me semble).

L'idée est d'utilisé la méta-prog pour que le compilateur génère le plus de calcul possible, ca allègera le traville de l'ordinateur au moment d la compilation (en gros il n'aura plus que des addition et multiplication à faire pour calcul chaque élément de la matrice et pas à faire les calcul un à un.

Bonjour,
La majorité des lib d'algèbre linéaire se basent sur blas (level 1,2 et 3) et lapack. Ce sont des librairies plutôt bas niveau écrites en fortran. Ces librairies sont des dinosaures du calcul matriciel, elles sont terriblement rapides et débuggées.
Ne pas les utiliser serait - à mon avis - une faute de goût :lol:. Il y a des passerelles avec le C++ via les numeric bindings et boost uBlas.
Bon courage,
Damien.

Bonsoir,

les solutions à adopter dépendent de tes propres besoins : que souhaites-tu faire exactement? Tu veux manipuler des matrices denses, des matrices structurées, des matrices creuses? Souhaitent faire des factorisations (LU,QR,SVD), résoudre des systèmes linéaires, calculer des valeurs propres?

Il faut décrire un peu ton projet avant d'avoir une réponse adaptée à ta situation.

Mon but n'est pas de recréer une bibliothèque qui surpassera de loin (ou pas ;)) les bibliothèques existantes.

C'est principalement à titre pédagogique. Je suis curieux de comprendre profondément comment marchent les outils que j'utilise tous les jours et tous les jours je travail avec matlab :lol:.

Bref c'est un prétexte pour apprendre. Mais je veux que cette bibliothèque soit exploitable sans avoir à rougir (trop) de celles déjà existante.

Je tiens à pouvoir utiliser cette bibliothèque sur des jeux de données de taille importante, pouvoir faire des transformations (Fourier, ondelette, hough,...) mais aussi des résolutions d'optimisation (LMS,...) et pour tout cela j'ai fermement l'intention d'utiliser les factorisations LU,QR,SVD.

Je ne vais pas plus loin dans les projets, il y à déjà la du travail pour ...

Citation:

---

Envoyé par omc24

Bonjour,
La majorité des lib d'algèbre linéaire se basent sur blas (level 1,2 et 3) et lapack. Ce sont des librairies plutôt bas niveau écrites en fortran. Ces librairies sont des dinosaures du calcul matriciel, elles sont terriblement rapides et débuggées.
Damien.

---

Je ne connais absolument pas le fortran mais j'en ai beaucoup entendu parlé en bien et en mal. Son plus gros défaut: de toutes mes études en école je n'en ai jamais entendu parlé de la part de mes professeurs.

Une bibliothèque basé sur ce langage ne risque-t-elle pas de ne pas bien évoluer avec les avancées matérielles futur? Notamment avec l'utilisation du calcul GPU?

à part qu'il faut voir à arreter de se palucher sur les GPU, les gens de LAPACK n'ont pas attendu pour ocntinuer à le faire évoluer (cf PLASMA, MAGMA et autre SCALAPACK)

Citation:

---

Envoyé par Joel F

à part qu'il faut voir à arreter de se palucher sur les GPU, les gens de LAPACK n'ont pas attendu pour ocntinuer à le faire évoluer (cf PLASMA, MAGMA et autre SCALAPACK)

---

:lol:, youpi ! enfin quelqu'un de raisonnable ! Je partage l'opinion...

Citation:

---

Envoyé par wurst

Je ne connais absolument pas le fortran mais j'en ai beaucoup entendu parlé en bien et en mal. Son plus gros défaut: de toutes mes études en école je n'en ai jamais entendu parlé de la part de mes professeurs.

---

Ce n'est pas un problème, avec les numeric bindings, tu peux manipuler tout ce beau monde en C++. Sinon comme je l'ai dit, ces librairies sont vraiment des "bestOf". Je dirais même qu'il n'y a pas mieux.
La majorité des applications/libs de calcul les utilise (y compris Matlab).

Bonjour,

Citation:

---

Envoyé par wurst

Mon but n'est pas de recréer une bibliothèque qui surpassera de loin (ou pas ;)) les bibliothèques existantes.

C'est principalement à titre pédagogique. Je suis curieux de comprendre profondément comment marchent les outils que j'utilise tous les jours et tous les jours je travail avec matlab :lol:.

---

Tu peux déjà commencer à regarder les algorithmes de matlab en utilisant la commande 'open'. Des choix d'implémentation faits dans matlab sont particulièrement discutables donc essaye aussi de voir du côté des livres.

Citation:

---

Envoyé par wurst

Bref c'est un prétexte pour apprendre. Mais je veux que cette bibliothèque soit exploitable sans avoir à rougir (trop) de celles déjà existante.

---

Si tu te limites aux matrices denses, tu ne devrais pas avoir trop de problèmes pour faire une bibliothèque honnête.

Citation:

---

Envoyé par wurst

Je tiens à pouvoir utiliser cette bibliothèque sur des jeux de données de taille importante, pouvoir faire des transformations (Fourier, ondelette, hough,...) mais aussi des résolutions d'optimisation (LMS,...) et pour tout cela j'ai fermement l'intention d'utiliser les factorisations LU,QR,SVD.

---

Pour les transformations de Fourier etc, il faut que tu cherches à faire un noyau d'algèbre linéaire dense très performant. Dans un premier temps, je te conseille de faire un noyau non optimisé pour mettre en place tes structures de données. Il faut que tu codes des algorithmes en O(n^2), typiquement des produits entre matrices et vecteurs, et en O(n^3), typiquement des produits entre matrices. Il ne s'agit que d'imbrications de boucles mais il faut faire attention à l'ordre d'imbrication pour ne pas détériorer la localité des données.

La lecture de la thèse de David Parello me semble être pertinente pour bien programmer ton noyau :
http://perso.univ-perp.fr/david.pare...ello_these.pdf
Tu y trouveras plusieurs techniques de transformations de programmes (à faire manuellement dans ton cas) pour améliorer les performances de ton code de façon simple et pédagogique.

Pour les factorisations matricielles, c'est un peu plus compliqué car il faut les implémenter de façon stable. Tout cela va dépendre des propriétés numériques de tes matrices (symétriques ou non, hermitiennes ou non, définies ou non, positives ou non, de plein rang ou non, nombre de conditionnement). De façon grossière, voici les conditions à remplir pour tes factorisations :

Cholesky : matrice hermitienne définie positive à un signe près
Cholesky+pivotage : matrice hermitienne indéfinie
Gauss : matrice n'admettant que des sous-matrices inversibles
Gauss+pivotage : matrice inversible
QR : matrice de plein rang colonne
SVD : aucune condition

Le pivotage consiste à permuter des lignes ou des colonnes pour éviter les divisions par zéro. Dans le cas des factorisations LU (Cholesky et Gauss), l'approche moderne consiste à pivoter suivant des blocs 2*2.

Si tu veux faire cela très proprement, je te conseille de lire les notes de travail de Lapack pour les algos qui t'intéressent :
http://www.netlib.org/lapack/lawns/downloads/
Les personnes qui les écrivent sont des spécialistes mondialement reconnus de l'algèbre linéaire numérique (Dongarra et Demmel en particulier).

Pour les factorisations LU, si tu es très motivé, tu peux jeter un oeil aux différentes versions du processus d'élimination de Gauss et Cholesky (left-looking, right-looking, frontal, multifrontal, supernodal pour matrices creuses).

Pour QR, il ne faut surtout pas implémenter l'algorithme de Gram-Schmidt mais lui préférer l'algorithme de Householder (cas général) ou celui de Givens (matrices structurées) qui sont beaucoup plus stables.

Citation:

---

Envoyé par omc24

La majorité des lib d'algèbre linéaire se basent sur blas (level 1,2 et 3) et lapack. Ce sont des librairies plutôt bas niveau écrites en fortran. Ces librairies sont des dinosaures du calcul matriciel, elles sont terriblement rapides et débuggées.

---

C'est effectivement une bonne idée que tu ailles voir comment sont faites ces bibliothèques. BLAS n'évolue plus beaucoup et est un peu dépassée en terme de conception et de performances mais soulève des problématiques intéressantes, notamment au niveau de l'orientation des matrices (par ligne ou par colonne), de la vectorisation des calculs et de la nomenclature des routines. Mais, puisque tu es en C++, autant que tu essayes de bénéficier des expression templates qui généralisent la raison d'être initiale de BLAS (au-delà de la volonté d'avoir voulu faire un standard).

LAPACK évolue encore mais n'est pas exempt de bogues.

BLAS et LAPACK ne sont pas particulièrement rapides et ce n'est pas très difficile de les surpasser de ce point de vue. La bibliothèque ATLAS te permet par exemple de générer des versions de ces bibliothèques optimisées pour ta machine. Elle fournit également une version optimisée pour la transformation de Fourier. Pour le C++, tu peux aussi regarder du côté des instructions SIMD et la bibliothèque Eigen qui semble très bien mais que je n'ai pas eu l'occasion de tester.

Citation:

---

Envoyé par wurst

Je ne connais absolument pas le fortran mais j'en ai beaucoup entendu parlé en bien et en mal. Son plus gros défaut: de toutes mes études en école je n'en ai jamais entendu parlé de la part de mes professeurs.

---

C'est parce que tu n'as pas fait d'études en calcul scientifique. C'est un langage plutôt simple d'accès donc tu peux t'y mettre tout seul sans peiner.

Citation:

---

Envoyé par wurst

Une bibliothèque basé sur ce langage ne risque-t-elle pas de ne pas bien évoluer avec les avancées matérielles futur? Notamment avec l'utilisation du calcul GPU?

---

Le fortran étant le langage par excellence pour le calcul scientifique, tu n'as pas de soucis à te faire pour l'évolution de ta bibliothèque vis-à-vis des avancées matérielles : des personnes se chargeront de fournir les outils adéquats avant même que tu n'ais réalisé ton besoin. Les personnes qui développent BLAS, LAPACK, ATLAS, MAGMA, PLASMA, MPI, etc, ont accès aux nouvelles technologies bien avant toi parce qu'ils travaillent sur des supercalculateurs et ont les fonds pour investir (et ils travaillent tous ensemble de surcroît).

Bonjour,

Je n'ai pas lu tous les messages, mais je vais essayer de répondre rapidement.

Si tu veux faire une bonne librairie de matrices, il faut que tu saches que tu n'auras pas le temps (ni la compétence?) de refaire une librairie complète à la uBlas.

Il faut que tu saches précisément ce dont tu as besoin et que tu l'implémentes.

En gros, la base, c'est d'avoir une classe vecteur rapide. A partir de ça, tu les utilises pour créer tes matrices comme tu le ferais avec les std::vector.
Ensuite il te faut une classe matrice, et puis une classe sparsematrix.
Il faut alors bien penser à mettre tous les operateurs qui peuvent t'être utiles, (produit matrice vecteur, etc...).

Une fois que tu auras ça, tu pourras réflechir à faire des algorithmes. Là il faut que tu te poses sérieusement la question du design.
Est ce qu'ils ont besoin d'être tous regroupés qq part? Est ce qu'ils doivent faire partie de tes classes de base? Est ce qu'ils doivent supporter tout types de données (si oui template ou pas?)

Mais le plus important, c'est vraiment ta classe vecteur. Je te donnerai le conseil de faire une classe du type :

Code:

---

1 2 3 4 5 6 7

class { int _size ; float* _data ; public: // some methods } ;

---

pour avoir un alignement mémoire de qualité, et aucune vérification de type STL_BOUND_CHECK qui pourrait plomber tes perfs.

Là il faut que tu rédiges des tests pour vérifier que ta mémoire est bien alouée/libérée...

Pour tes matrices, il faut (à mon avis) que ce soit un seul vecteur de données et les dimensions (largeur + hauteur)

Pour les sparsematrix, là il faut te documenter sur le stockage et la compression de matrices creuses. A un moment tu devras choisir entre temps d'accès aux données ou place mémoire. A voir ce dont tu as besoin...

Une fois que t'en es là : il faut que tu fasses des benchmarks sur toutes tes fonctions et que tu les optimises. (very sleepy peut être utile ou gprof sous linux). Tu seras surpris de gagner parfois x2 ou x3...

Pour les expressions templates dont parlait qqn ici, je pense que tu as 5 min avant d'y penses ;)

Et une fois que ta lib marche : achete le compilo intel et distribue ton .lib :ccool:

Citation:

---

Envoyé par regis.portalez

En gros, la base, c'est d'avoir une classe vecteur rapide. A partir de ça, tu les utilises pour créer tes matrices comme tu le ferais avec les std::vector.
Ensuite il te faut une classe matrice, et puis une classe sparsematrix.
Il faut alors bien penser à mettre tous les operateurs qui peuvent t'être utiles, (produit matrice vecteur, etc...).

---

Non jamais, un vector de vector va peter ta localite de cache.

Citation:

---

Envoyé par regis.portalez

Une fois que tu auras ça, tu pourras réflechir à faire des algorithmes. Là il faut que tu te poses sérieusement la question du design.
Est ce qu'ils ont besoin d'être tous regroupés qq part? Est ce qu'ils doivent faire partie de tes classes de base? Est ce qu'ils doivent supporter tout types de données (si oui template ou pas?)

---

La question ne se pose pas: separation algo/data OBLIGATOIRE

Citation:

---

Envoyé par regis.portalez

Mais le plus important, c'est vraiment ta classe vecteur. Je te donnerai le conseil de faire une classe du type :

Code:

---

1 2 3 4 5 6 7

class { int _size ; float* _data ; public: // some methods } ;

---

pour avoir un alignement mémoire de qualité, et aucune vérification de type STL_BOUND_CHECK qui pourrait plomber tes perfs.

---

Pure LU et FUD. Apprendre a compiler en -NDEBUG -O3 devrais t'apporter plein de chose.

Citation:

---

Envoyé par regis.portalez

Pour tes matrices, il faut (à mon avis) que ce soit un seul vecteur de données et les dimensions (largeur + hauteur)

---

non cf mes posts, la bonne maniere est la maniere NRC.

Citation:

---

Envoyé par regis.portalez

Pour les expressions templates dont parlait qqn ici, je pense que tu as 5 min avant d'y penses ;)

---

Non, c'est des le debut.

Petite remarque, parce que ça fait plusieurs fois qu'il est cité, le code de Joel pour les tableaux NRC est "faux" :
s/delete/delete []

lequel ?, je croyais avoir corriger dans un autre codepad.

Celui cité un peu plus haut sur hfr. (et je sais qu'il avait été recopié ici aussi)

Citation:

---

Envoyé par Goten

Celui cité un peu plus haut sur hfr. (et je sais qu'il avait été recopié ici aussi)

---

Ouias.

J'ai trouve mieux, ls dit code NRC + du bench: http://boost.codepad.org/PPN6uRYy

A Joel :

je n'ai jamais parlé de vector de vector... c'est effectivement catastrophique...

Et quand je dis que les expressions templates sont à avoir à la fin, c'est parce qu'avant de les faire, il vaut mieux savoir faire le reste... Surtout si c'est pour rédiger une bibliothèque utilisée en local.

Les expressions templates sont là pour corriger les défauts d'utilisation de la lib par ses utilisateurs.

Par exemple, il faut être naïf pour écrire :

A = (B*C) + D - 3*D.norm()*E...

Quand on écrit une lib pour soi, on sait qu'on ne va pas écrire ça.
On va faire soi même la boucle en dur, comme un grand

Quand on la livre à qqn, il faut prévoir le cas où le mec va pas s'embêter à réfléchir...

Citation:

---

Envoyé par regis.portalez

Par exemple, il faut être naïf pour écrire :

A = (B*C) + D - 3*D.norm()*E...

Quand on écrit une lib pour soi, on sait qu'on ne va pas écrire ça.
On va faire soi même la boucle en dur, comme un grand

---

non jamais ...

Citation:

---

Envoyé par regis.portalez

Et quand je dis que les expressions templates sont à avoir à la fin, c'est parce qu'avant de les faire, il vaut mieux savoir faire le reste... Surtout si c'est pour rédiger une bibliothèque utilisée en local.

---

Le problème, c'est que les expression template ne sont pas factorisable dans un code déjà écrit s'ils n'a pas été pensé pour. Si on veut s'en servir, on est obligé de s'en servir dès le début.

Citation:

---

Envoyé par regis.portalez

Les expressions templates sont là pour corriger les défauts d'utilisation de la lib par ses utilisateurs.

---

Non. Ils sont là pour corriger des problèmes d'instanciation de classes et de combinaison des opérateurs associés.

Citation:

---

Envoyé par regis.portalez

Par exemple, il faut être naïf pour écrire :

A = (B*C) + D - 3*D.norm()*E...

---

Naïf ? Je dirais qu'il faut être normal pour écrire un tel code. Le principe d'une lib, c'est justement de permettre d'écrire le code de manière claire et lisible.

Citation:

---

Envoyé par regis.portalez

Quand on écrit une lib pour soi, on sait qu'on ne va pas écrire ça.

---

Ben si, justement. C'est le but même :) Sinon, autant ne pas écrire de lib :)

Citation:

---

Envoyé par regis.portalez

On va faire soi même la boucle en dur, comme un grand

---

Pourquoi pas ? Ca n'est qu'une duplication de code comme une autre :)

Citation:

---

Envoyé par regis.portalez

Quand on la livre à qqn, il faut prévoir le cas où le mec va pas s'embêter à réfléchir...

---

Le client a droit à une lib simple d'utilisation et efficace, quand moi je devrais m'ennuyer avec un truc odieux ? C'est pas très logique tout ça :)

Youhou les gens, le type du début a pas dit qu'il voulait distribuer sa lib au monde entier, il veut juste avoir un truc qui marche et qui lui fasse LU en à peu près autant de temps que MATLAB.

C'est pas la peine de le gaver avec des templates de fous.

Il voudrait faire un vrai lib distribuable (ou même commercialisable) je dis pas, mais là c'est pas vraiment le cas...

à voir le deuxième message du mec : "j'ai découvert le principe des Template que récement, je suis pas sur d'en connaitre assez pour comprendre ce que tu veux me dire."

Si vous voulez vous faire un débat sur la meilleur lib de matrice de niveau pro allez-y, mais ça ne répond pas au sujet.

Citation:

---

Le problème, c'est que les expression template ne sont pas factorisable dans un code déjà écrit s'ils n'a pas été pensé pour. Si on veut s'en servir, on est obligé de s'en servir dès le début.

---

Oui, mais il faut d'abord savoir comment on rédige des classes et des templates avant de le faire.
Sinon, c'est comme si tu disais à un type qui sait pas comme marche un moteur de fabriquer une voiture en commençant par la boîte automatique...

Bonsoir,

je pense au contraire que chacun cherche à répondre à la problématique initiale qui spécifie bien un problème de temps de calcul pour une matrice de Fourier comparativement à ce que fait Matlab. Pour approcher les performances de ce dernier, il faut s'inspirer des bibliothèques externes qu'il utilise (blas & lapack optimisés pour la réutilisation des mémoires caches + umfpack si ma mémoire est bonne). La moindre des choses c'est bien d'inviter les "candidats à faire aussi bien en C++" à utiliser les expression templates pour faire du plus haut niveau que fortran tout en maintenant des performances de calcul honnêtes. En ce qui concerne la factorisation LU proprement dite, imbriquer trois boucles ne suffira pas.

Citation:

---

Envoyé par regis.portalez

Youhou les gens, le type du début a pas dit qu'il voulait distribuer sa lib au monde entier, il veut juste avoir un truc qui marche et qui lui fasse LU en à peu près autant de temps que MATLAB.

<snip>

à voir le deuxième message du mec : "j'ai découvert le principe des Template que récement, je suis pas sur d'en connaitre assez pour comprendre ce que tu veux me dire."

<snip>

Oui, mais il faut d'abord savoir comment on rédige des classes et des templates avant de le faire.

---

Ouai sauf que c'est totalement contradictoire, si tu maitrises pas ce genre de chose alors l'objectif fixé est pas réalisable.

Petite parenthèse : les NRC ont une licence fort pénible. Je ne sais pas comment s'interprète le modèle d'allocation vis à vis de la licence, mais c'est un détail auquel il faut faire attention.
J'ai eu un exemple de lib opensource (OTB) dernièrement qui a dégagé ses dépendances aux NRC.

pour avoir eu le probleme icic, les algos numeriques sont sous licenses, l'allocateur memoire ets "Public Domain"

Merci pour l'info. C'est bon à savoir.