lotto Soluzione CUDA di Axel a banda sparsa = b per varie b

Question

Ho una matrice a banda sparsa A e mi piacerebbe (diretto) risolvere Ax = b. Ho circa 500 vettori b, quindi mi piacerebbe risolvere per i corrispondenti 500 x. Sono nuovo di zecca per CUDA, quindi sono un po 'confuso su quali opzioni ho a disposizione.

cuSOLVER ha un risolutore diretto batch cuSolverSP per sparse A_i x_i = b_i utilizzando QR here. (Anche a me starei bene con LU perché A è decentemente condizionato.) Comunque, per quanto posso dire, non posso sfruttare il fatto che tutti i miei A_i sono uguali.

Un'opzione alternativa è quella di determinare prima una fattorizzazione LU sparsa (QR) sulla CPU o GPU, quindi eseguire in parallelo la backsubstitution (rispettivamente, backsub e matrix mult) sulla GPU? Se cusolverSp< t >csrlsvlu() è per un b_i, esiste un modo standard per eseguire questa operazione in batch per più b_i?

Infine, poiché non ho intuito per questo, dovrei aspettarmi una accelerazione su una GPU per una di queste opzioni, dato il sovraccarico necessario? x ha lunghezza ~ 10000-100000. Grazie.

Answer 1

Attualmente sto lavorando a qualcosa di simile. Ho deciso di avvolgere fondamentalmente il gradiente coniugato e il livello-0 incompleto di campioni di utilità per il solutore di gradienti precondizionati di Cholesky precondizionati forniti con l'SDK CUDA in una piccola classe.

Li potete trovare nella vostra directory CUDA_HOME nel percorso: samples/7_CUDALibraries/conjugateGradient e /Developer/NVIDIA/CUDA-samples/7_CUDALibraries/conjugateGradientPrecond

In sostanza, si dovrebbe caricare la matrice nella memoria del dispositivo una volta (e per ICCG, calcolare la/analisi corrispondente balsamo matrice), allora chiama il kernel di risoluzione con diversi vettori b.

Non so cosa sia l'aspetto della struttura della matrice a matrice, ma se è simmetrica e diagonalmente dominante (le bande diagonali lungo ogni riga e colonna sono di segno opposto rispetto alla diagonale e la loro somma è inferiore alla ingresso diagonale) o positivo definito (nessun autovettore con autovalore di zero), allora CG e ICCG dovrebbero essere utili. In alternativa, i vari algoritmi multigrid sono un'altra opzione se si è disposti a passare attraverso la codifica.

Se la matrice è solo semi-definita positiva (ad es. Ha almeno un autovettore con autovalore pari a zero), è comunque possibile evitare di utilizzare CG o ICCG a condizione che: 1) La mano destra lato (i vettori b) sono ortogonali allo spazio nullo (spazio nullo che significa autovettori con autovalore pari a zero). 2) La soluzione ottenuta è ortogonale allo spazio nullo.

È interessante notare che se si dispone di uno spazio Null non banale, i diversi risolutori numerici possono fornire risposte diverse per lo stesso sistema esatto. Le soluzioni finiranno per differire da una combinazione lineare dello spazio nullo ... Questo problema mi ha causato molte ore di debug e frustrazioni prima che finalmente venissi scoperto, quindi è bene esserne consapevole.

Infine, se la matrice ha un Circulant Band structure, si potrebbe prendere in considerazione l'utilizzo di un risolutore veloce basato su trasformata di Fourier (FFT). I risolutori numerici basati su FFT possono spesso fornire prestazioni superiori nei casi in cui sono applicabili.

Answer 2

Se non ti dispiace andare con una libreria open-source, si potrebbe anche verificare CUSP: CUSP Quick Start Page

Esso dispone di una suite abbastanza decente di solutori, tra cui alcuni metodi precondizionati: CUSP Preconditioner Examples

Il precondizionatore di aggregazione levigato (una variante del multigrid algebrico) sembra funzionare molto bene purché la GPU disponga di memoria interna sufficiente per esso.