[C++] Librerie di algebra lineare & HPC
Ciao a tutti,
spero sia la sezione giusta.
Per la mia tesi magistrale (analisi numerica) mi sto ritrovando ad utilizzare in modo intensivo C++ (con openMP per parallelizzare su CPUs), in particolare ho a che fare con matrici sparse. Come libreria di algebra lineare sto usando Eigen, ma so che in rete ce ne sono altre di popolari come PETSc e Armadillo che non ho mai usato personalmente.
Volevo chiedere, per chi ha avuto la possibilità oppure è del settore, quale preferite/consigliate/usate e chiaramente perché.
spero sia la sezione giusta.
Per la mia tesi magistrale (analisi numerica) mi sto ritrovando ad utilizzare in modo intensivo C++ (con openMP per parallelizzare su CPUs), in particolare ho a che fare con matrici sparse. Come libreria di algebra lineare sto usando Eigen, ma so che in rete ce ne sono altre di popolari come PETSc e Armadillo che non ho mai usato personalmente.
Volevo chiedere, per chi ha avuto la possibilità oppure è del settore, quale preferite/consigliate/usate e chiaramente perché.
Risposte
Non ho molta esperienza con le matrici sparse, e di quelle librarie conosco solo Eigen (ma sono anni che non ho a che fare con BLAS). Da quel che ho capito, internamente funzionano in modo simile e personalmente immagino che abbiano performance similari.
Se già conosci la sintassi di Eigen è probabilmente più conveniente continuare ad usarlo. Ma immagino che dipenda da quale sia quella più usata nel tuo specifico settore. Le altre due librerie mi sembra che siano qualcosa di più di librerie di algebra lineare.
Se già conosci la sintassi di Eigen è probabilmente più conveniente continuare ad usarlo. Ma immagino che dipenda da quale sia quella più usata nel tuo specifico settore. Le altre due librerie mi sembra che siano qualcosa di più di librerie di algebra lineare.
Io ho un po' di esperienza con queste librerie.
Eigen è una libreria con un bel design, ma non molto facile da usare, almeno all'inizio. Costrutti apparentemente innocui possono generare errori molto criptici. Eigen non supporta il calcolo parallelo a memoria distribuita [MPI], quindi sarà sempre limitato ad un solo nodo.
PETSc è una libreria ad oggetti scritta in C ma si può usare anche da C++ e da Python. È molto conosciuta, ha generalmente buone prestazioni ed un'interfaccia decisamente più user-friendly di Eigen. È una libreria per simulazioni giganti e supporta solo MPI, non OpenMP. Chiaramente la puoi usare anche su un solo nodo con più processi.
Trilinos è il rivale più noto di PETSc. Trilinos è scritto in C++ ed è una libreria molto più grande di PETSc per quanto riguarda il range di utilizzo, infatti contiene sia moduli di algebra lineare sia moduli di ottimizzazione, partizionamento di grafi e tante altre cose. A giudicare dalla documentazione, direi che Trilinos supporta sia OpenMP sia MPI, ma non l'ho mai usata personalmente.
PETSc e Trilinos sono i due giganti del settore, quindi magari vale la pena esplorare una di queste due.
Eigen è una libreria con un bel design, ma non molto facile da usare, almeno all'inizio. Costrutti apparentemente innocui possono generare errori molto criptici. Eigen non supporta il calcolo parallelo a memoria distribuita [MPI], quindi sarà sempre limitato ad un solo nodo.
PETSc è una libreria ad oggetti scritta in C ma si può usare anche da C++ e da Python. È molto conosciuta, ha generalmente buone prestazioni ed un'interfaccia decisamente più user-friendly di Eigen. È una libreria per simulazioni giganti e supporta solo MPI, non OpenMP. Chiaramente la puoi usare anche su un solo nodo con più processi.
Trilinos è il rivale più noto di PETSc. Trilinos è scritto in C++ ed è una libreria molto più grande di PETSc per quanto riguarda il range di utilizzo, infatti contiene sia moduli di algebra lineare sia moduli di ottimizzazione, partizionamento di grafi e tante altre cose. A giudicare dalla documentazione, direi che Trilinos supporta sia OpenMP sia MPI, ma non l'ho mai usata personalmente.
PETSc e Trilinos sono i due giganti del settore, quindi magari vale la pena esplorare una di queste due.
Grazie ad entrambi per le risposte 
@vict Sì ho molta familiarità con Eigen, anche perché come prima librerai di algebra lineare mi sembrava quella più user-friendly.
@Raptorista Cavolo, non sapevo che Eigen non supportasse MPI/OpenMP. Questo è un problema. Credo che opterò per PETSc, nonostante supporti solo MPI. Anche perché mi pare di vedere (leggendo su scicomp) che per quanto riguarda sistemi lineari sparsi ecc. sia molto usata da gente del settore. Se posso chiedere, hai avuto problemi nell'utilizzare MPI con PETSc?

@vict Sì ho molta familiarità con Eigen, anche perché come prima librerai di algebra lineare mi sembrava quella più user-friendly.
@Raptorista Cavolo, non sapevo che Eigen non supportasse MPI/OpenMP. Questo è un problema. Credo che opterò per PETSc, nonostante supporti solo MPI. Anche perché mi pare di vedere (leggendo su scicomp) che per quanto riguarda sistemi lineari sparsi ecc. sia molto usata da gente del settore. Se posso chiedere, hai avuto problemi nell'utilizzare MPI con PETSc?
"feddy":
@Raptorista Cavolo, non sapevo che Eigen non supportasse MPI/OpenMP. Questo è un problema.
Supporta solo OpenMP.
"feddy":
Credo che opterò per PETSc, nonostante supporti solo MPI.
È una buona scelta. Su un solo nodo non noterai la differenza comunque.
"feddy":
Anche perché mi pare di vedere (leggendo su scicomp) che per quanto riguarda sistemi lineari sparsi ecc. sia molto usata da gente del settore.
Né più né meno di Trilinos.
"feddy":
Se posso chiedere, hai avuto problemi nell'utilizzare MPI con PETSc?
Non propriamente problemi, però in certi casi devi fare attenzione a rendere i vettori compatibili con le matrici che vuoi usare, avere il giusto sparsity pattern, cose così. È una libreria in C, quindi molte cose vanno fatte a mano.
"Raptorista":
Supporta solo OpenMP.
Ok, a questo punto provo ad iniziare con Eigen giusto per scrupolo (ho molta più familiarità con openMP), ma ho deciso che utilizzerò PETSc alla fine.
"Raptorista":
È una buona scelta. Su un solo nodo non noterai la differenza comunque.
Certamente, il fatto è che posso utilizzare il cluster di ateneo, mi ero dimenticato di dirlo in effetti

Un'ultima cosa: su MPI ho veramente poca esperienza, conosci per caso qualche buon testo "introduttivo" sul suo utilizzo (con c++) ?
"feddy":
Un'ultima cosa: su MPI ho veramente poca esperienza, conosci per caso qualche buon testo "introduttivo" sul suo utilizzo (con c++) ?
A suo tempo ho usato "Pacheco - Parallel programming with MPI". È un po' datato ma dovrebbe essere ancora utile. Le implementazioni più comuni di MPI sono in C, non in C++, ma si possono usare da codice C++.
Avrai effettivamente a che fare con più processori? Insomma su che sistema intendi usarlo?
Grazie Raptorista 
@vict Sì ho a che fare con un Beowulf cluster

@vict Sì ho a che fare con un Beowulf cluster
Ok, quindi devi usare MPI. D'altra parte non è strettamente necessario che la tua libraria supporti internamente MPI. Insomma puoi dividere il problema tra i nodi manualmente e poi chiamare le funzioni di algebra lineare internamente ai nodi. Se invece la libreria lo usa internamente potrebbe non essere strettamente necessario che tu sappia usare bene MPI. Dovresti cominciare con i tutorial di PETSc secondo me.
Ho guardato i tutorial/esempi/esericizi fino ad ora. Molto molto interessanti e belli, anche se devo spenderci un bel po' ora per impratichirmi.
Riguardo all'inizio del tuo ultimo post, quello che devo fare idealmente ora è risolvere sistemi lineari di grandi dimensioni (taglia 10^7) che sono dentro un ciclo for. Questi sistemi sono tutti indipendenti gli uni dagli altri e alla fine mi serve la somma delle varie soluzioni, perciò il problema è banalmente parallelizzabile. Sono inoltre relativamente pochi (sono $8$).
Credo che la cosa da fare qui sia semplicemente risolvere su ciascun nodo un sistema lineare e poi alla fine sommare il tutto. Da neofita di MPI mi viene da dire una cosa del tipo "8 processi in cui vengono calcolate le soluzioni e alla fine mando tutto al processo $0$ per sommare".
Oppure c'è un modo più efficiente (almeno in teoria) ?
Riguardo all'inizio del tuo ultimo post, quello che devo fare idealmente ora è risolvere sistemi lineari di grandi dimensioni (taglia 10^7) che sono dentro un ciclo for. Questi sistemi sono tutti indipendenti gli uni dagli altri e alla fine mi serve la somma delle varie soluzioni, perciò il problema è banalmente parallelizzabile. Sono inoltre relativamente pochi (sono $8$).
Credo che la cosa da fare qui sia semplicemente risolvere su ciascun nodo un sistema lineare e poi alla fine sommare il tutto. Da neofita di MPI mi viene da dire una cosa del tipo "8 processi in cui vengono calcolate le soluzioni e alla fine mando tutto al processo $0$ per sommare".
Oppure c'è un modo più efficiente (almeno in teoria) ?
Va bene calcolare le soluzioni degli otto sistemi lineari indipendentemente, ma usando ben più di un processo per sistema.
La cosa veramente più banale è probabilmente risolvere gli 8 sistemi in serie, dividendo ciascun sistema su tutti i processi [mille processi mi sembra un numero ragionevole] e su ciascun processo accumulare la soluzione, poi salvarla in parallelo. Puoi fare questa cosa senza mai sporcarti le mani con MPI, usando le funzionalità di I/O parallelo di PETSc.
Se hai risorse in abbondanza, allora puoi usare 8000 processi, divisi in gruppi da mille, e fare tutto in parallelo, e questo sarà probabilmente più efficiente dell'altro approccio dividendo ciascun sistema tra 8000 processi, però è un po' più complicato da implementare.
La cosa veramente più banale è probabilmente risolvere gli 8 sistemi in serie, dividendo ciascun sistema su tutti i processi [mille processi mi sembra un numero ragionevole] e su ciascun processo accumulare la soluzione, poi salvarla in parallelo. Puoi fare questa cosa senza mai sporcarti le mani con MPI, usando le funzionalità di I/O parallelo di PETSc.
Se hai risorse in abbondanza, allora puoi usare 8000 processi, divisi in gruppi da mille, e fare tutto in parallelo, e questo sarà probabilmente più efficiente dell'altro approccio dividendo ciascun sistema tra 8000 processi, però è un po' più complicato da implementare.