ALGLIB liberia open source per il calcolo numerico

[dimaios]

Con il presente post vi informo che è stata rilasciata una nuova release della ALGLIB.

27.03.2012 ALGLIB 3.5.0 is released
http://www.alglib.net/

E' una libreria molto interessante che uso spesso con profitto e soddisfazione.
In passato ho valutato molte alternative ma alla fine ho deciso di usarla e non mi sono ancora pentito.
E' cross-platform e supporta diversi linguaggi di programmazione.

[DarwinNE]

E' molto interessante. Che tipo di calcoli ti trovi a fare?
Io utilizzo LAPACK per trattare alcuni problemi agli autovalori. La cosa un po' delicata è che il mio programma è in C++ e bisogna quindi interfacciarsi al Fortran, cosa che si fa bene, ma che alle volte è un tantino scomoda. Sai mica come si comporta ALGLIB per gli autovalori? E' facilmente scalabile su sistemi multicore? Come potenza pura e dura, è comparabile a LAPACK accoppiato a ATLAS?

[dimaios]

Ciao DarwinNE
rispondo brevemente alle domande.

E' molto interessante. Che tipo di calcoli ti trovi a fare?

In genere algoritmi real-time per la misura di oggetti in campo industriale.
Le immagini possono provenire da telecamere o più frequentemente da laser scanner 2D.
Inoltre mi capita spesso di dover risolvere problemi di ottimizzazione matematica , controlli automatici e risoluzione di equazioni differenziali non lineari.

Io utilizzo LAPACK per trattare alcuni problemi agli autovalori. La cosa un po' delicata è che il mio programma è in C++ e bisogna quindi interfacciarsi al Fortran, cosa che si fa bene, ma che alle volte è un tantino scomoda.

Se provieni da LAPACK non c'è problema perché la relazione tra Alglib e Lapack è strettissima.
http://www.kitware.com/InfovisWiki/index.php/Linear_Algebra_Packages
Ho impiegato la Alglib sia in C++ che in C# e non ho avuto alcun problema.
Anche io ho utilizzato in passato Lapack.

Sai mica come si comporta ALGLIB per gli autovalori?

E' difficile dare una risposta univoca.
Per quanto riguarda le prove che ho condotto si è comportata molto bene ma le matrici non erano enormi e non erano sparse.
Hanno migliorato molto gli algoritmi nelle ultime versioni.
Se ho dubbi di condizionamento utilizzo l' SVD per evitare sorprese.

E' facilmente scalabile su sistemi multicore? Come potenza pura e dura, è comparabile a LAPACK accoppiato a ATLAS?

Questo è un argomento molto ostico.
Nel forum dell' Alglib trovi alcune risposte più o meno esplicite tipo questa http://forum.alglib.net/viewtopic.php?f=2&t=377 ma non ho esperienze tali da poter esprimere un giudizio autorevole.

Nei casi di mio interesse eseguo lo split sui frame provenienti dai sensori e non sul singolo algoritmo dell'AlgLib per cui il parallelismo risulta di semplice implementazione.
Non ho esperienza con AlgLib su come parallelizzare per esempio la ricerca degli autovalori di una matrice sparsa di grandi dimensioni ( se è quello che vorresti fare ).

Per quanto riguarda la velocità di calcolo ho avuto grosse soddisfazioni ( compilato in modalità release ovviamente ).
Non saprei compararlo con Lapack accoppiato con Atlas, secondo me ti conviene effettuare di persona alcuni test. Se provieni dal mondo Lapack impieghi poco tempo per scrivere qualche riga di codice in C++ per caricare una matrice "cattiva" e fare un confronto in termini di performance.

[DarwinNE]

Per quanto riguarda le prove che ho condotto si è comportata molto bene ma le matrici non erano enormi e non erano sparse.
Hanno migliorato molto gli algoritmi nelle ultime versioni.

Molto interessante. Nel mio caso, però le matrici sono enormi e sono piene (complesse e non simmetriche), tant'é che il problema principale che limita l'algoritmo è l'occupazione di memoria, seguita dal tempo di calcolo.
Attualmente, utilizzo ZGEEVX. Quando ho da diagonalizzare più matrici alla volta, potrei parallelizzare il calcolo lanciando più diagonalizzazioni, ma ci sono altri problemi legati all'occupazione della memoria (e alla struttura del programma, che sarebbe da rivedere parzialmente). Attivando il multithreading a livello di ATLAS si riesce a dare una svegliata al tutto sui sistemi multiprocessore, e si risparmia un 20% sul tempo di calcolo, il che è già un ottimo risultato.
Utilizzo anche FFTW3, anche se i problemi con la FFT sono molto meno stringenti e potrei adottare algoritmi anche meno ottimizzati. Soluzioni come ALGLIB mi interessano, perché ho notato che una delle principali difficoltà che i colleghi incontrano usando il mio codice è che il tutto è piuttosto difficile da compilare.

Se provieni dal mondo Lapack impieghi poco tempo per scrivere qualche riga di codice in C++ per caricare una matrice "cattiva" e fare un confronto in termini di performance.

Mi sembra un buon consiglio. La strategia con cui gli elementi delle matrici sono memorizzati è la stessa del Fortran e di Lapack? Altrimenti, devo rivedere un bel po' di roba...

[dimaios]

Per avere un'idea della sintassi per dichiarare le matrici ( esempio con una principal component analysis ) è del tipo :

Codice: Seleziona tutto

// 6 samples of 3D data (x, y, z)
double values[18] = {1.0, 2.0, 3.0,
                     1.0, 2.0, 3.0,
                     1.0, 2.0, 3.0,
                     1.0, 2.0, 3.0,
                     1.0, 2.0, 3.0,
                     1.0, 2.0, 3.0};
 
// add this data to ALGLIB's format
alglib::real_2d_array ptInput;
ptInput.setcontent(6, 3, values);
 
// this guy gets passed in and is filled in with an integer status code
alglib::ae_int_t info;
 
// scalar values that describe variances along each eigenvector
alglib::real_1d_array eigValues;
 
// unit eigenvectors which are the orthogonal basis that we want
alglib::real_2d_array eigVectors;
 
// perform the analysis
pcabuildbasis(ptInput, 6, 3, info, eigValues, eigVectors);
 
// now the vectors can be accessed as follows:
double basis0_x = eigVectors[0][0];
double basis0_y = eigVectors[1][0];
double basis0_z = eigVectors[2][0];
 
double basis1_x = eigVectors[0][1];
double basis1_y = eigVectors[1][1];
double basis1_z = eigVectors[2][1];
 
double basis2_x = eigVectors[0][2];
double basis2_y = eigVectors[1][2];
double basis2_z = eigVectors[2][2];

Quindi passi facilmente dalla notazione vettoriale classica del C++ alla forma impiegata dalla Alglib.

Prima di tutto però chiederei al guru alcune infomazioni perché da quello che ho capito hai a che fare con matrici di numeri complessi : http://forum.alglib.net/viewtopic.php?f=2&t=285.

Sergey Bochkanov è sempre molto preciso e puntuale nelle risposte per cui domanderei a lui se il tuo problema è risolvibile in modo efficiente tramite Alglib.

[DarwinNE]

Quindi passi facilmente dalla notazione vettoriale classica del C++ alla forma impiegata dalla Alglib.

Molto interessante. Cerco di essere più esplicito sul mio problema, sperando di non essere OT e di non sfruttare troppo il tuo tempo e la tua disponibilità. Immagina di avere fra le mani una matrice nonsimmetrica complessa che occupa da sola il 70% della memoria della macchina che sta facendo girare il problema. Questa matrice è memorizzata sfruttando le convenzioni del Fortran. Se per convertirla nella forma richiesta da ALGLIB devo farne una copia, mi trovo ad avere un problema non facilmente risolvibile...

Nel caso, chiedo al forum di Alglib...

P.S. se vuoi saperne di più su quello che sto facendo e se hai accesso ad articoli universitari, è uscito da poco un articolo che ti segnalo:

http://www.opticsinfobase.org/josaa/abs ... a-29-3-367

[dimaios]

..... immagina di avere fra le mani una matrice nonsimmetrica complessa che occupa da sola il 70% della memoria della macchina che sta facendo girare il problema. Questa matrice è memorizzata sfruttando le convenzioni del Fortran. Se per convertirla nella forma richiesta da ALGLIB devo farne una copia, mi trovo ad avere un problema non facilmente risolvibile...

La copia fisica è assolutamente da evitare.
In realtà se dovessi affrontare il problema procederei nel seguente modo :
[1] Chiedo a Sergey Bochkanov se il problema è risolvibile con la AlgLib dando indicazioni su dimensioni e tipologia della matrice
[2] Scrivo qualche riga di codice offline per eseguire un test funzionale sulla ricerca di autovalori della matrice e confronto il risultato con quello generato da Lapack.
( salvo su disco la matrice generata dal codice Lapack ed i risultati e la ricarico con il programma C++ Alglib ).
[3] Lancio diverse simulazioni facendo eseguire i calcoli dal codice Lapack e Alglib per avere la certezza che le due soluzioni siano equivalenti.

Il codice sarebbe da riscrivere ed è per questo che una analisi preliminare costi-benefici deve essere assolutamente fatta a priori.

P.S. se vuoi saperne di più su quello che sto facendo e se hai accesso ad articoli universitari, è uscito da poco un articolo che ti segnalo:

http://www.opticsinfobase.org/josaa/abs ... a-29-3-367

Purtroppo non ho accesso al testo completo ma dall' abstract deduco quanto l'argomento possa essere interessante dal punto di vista teorico ed anche applicativo.
Complimenti ... è davvero un bel settore per la ricerca.

Fammi sapere se hai trovato un riscontro positivo con l' Alglib.

[DirtyDeeds]

Nel mio caso, però le matrici sono enormi e sono piene

Per curiosità, enormi quanto?

[DarwinNE]

Per curiosità, enormi quanto?

Beh, un 3000 elementi di lato. Anche di più, in certi casi. Dovrei controllare, ma mi pare di aver già utilizzato matrici che fanno 10k elementi di lato. In realtà, temo che il mio codice sia subottimale e che in qualche caso vengano utilizzate delle matrici temporanee che non riesco ad evitare

Intanto ringrazio dimaios per i consigli e le utili informazioni.

[DirtyDeeds]

Con 3000 per lato non dovresti avere grossi problemi di memoria, almeno credo. Il massimo con cui mi è capitato di lavorare è stato 4096*4096 sotto matlab, però erano matrici ragionevolmente sparse (non l'inversa, però).