Batch CUDA Lösung von spärlich gebändert Ax = b für verschiedene b

Question

Ich habe eine spärlich gebänderte Matrix A und ich möchte (direkt) lösen Ax = b. Ich habe ungefähr 500 Vektoren b, also würde ich gerne für die entsprechenden 500 x's auflösen. Ich bin neu bei CUDA, daher bin ich ein wenig verwirrt darüber, welche Möglichkeiten ich habe.

cuSOLVER hat einen Batch-Direktlöser cuSolverSP für Sparse A_i x_i = b_i mit QR here. (Ich würde auch mit LU in Ordnung sein, da A anständig konditioniert ist.) Allerdings kann ich, soweit ich das beurteilen kann, nicht die Tatsache ausnutzen, dass alle meine A_i gleich sind.

Wäre es eine Alternative, zuerst eine spärliche LU (QR) -Faktorisierung auf der CPU oder der GPU zu ermitteln und dann parallel die Backsubstitution (Backsub und Matrix mult) auf der GPU durchzuführen? Wenn cusolverSp< t >csrlsvlu() für einen b_i ist, gibt es eine Standardmethode, um diese Operation für mehrere b_i durchzuführen?

Schließlich, da ich keine Intuition dafür habe, sollte ich eine Beschleunigung auf einer GPU für eine dieser Optionen erwarten, angesichts der notwendigen Overhead? x hat eine Länge von ~ 10000-100000. Vielen Dank.

Answer 1

Ich arbeite derzeit an etwas ähnliches selbst. Ich entschied mich dafür, die mit dem CUDA SDK gelieferten konjugierten Gradienten- und Level-0-unvollständigen, cholesky-vorkonditionierten konjugierten Gradienten-Solver-Utility-Beispiele in eine kleine Klasse zu packen.

Sie können sie in Ihrem CUDA_HOME Verzeichnis unter dem Pfad: samples/7_CUDALibraries/conjugateGradient und /Developer/NVIDIA/CUDA-samples/7_CUDALibraries/conjugateGradientPrecond

Grundsätzlich würden Sie die Matrix in den Gerätespeicher geladen werden einmal (und für ICCG, berechnen die entsprechende Anlage/Matrix-Analyse), dann rufe den Solve-Kernel mit verschiedenen b-Vektoren auf.

Ich weiß nicht, wie Sie Ihre Matrix Bandstruktur aussehen, aber wenn es symmetrisch und entweder diagonal dominant ist (außerhalb diagonalen Bändern entlang jeder Zeile und Spalte sind entgegengesetzte Vorzeichen der Diagonale und ihre Summe ist kleiner als die Diagonaleintrag) oder positiv definit (keine Eigenvektoren mit einem Eigenwert von 0), dann sollten CG und ICCG nützlich sein. Alternativ dazu sind die verschiedenen Mehrgitter-Algorithmen eine weitere Option, wenn Sie bereit sind, sie zu codieren.

Wenn Ihre Matrix nur positiv semi-definit ist (zB hat mindestens einen Eigenvektor mit einem Eigenwert von Null), können Sie trotzdem mit CG oder ICCG durchkommen, solange Sie sicherstellen, dass: 1) Die rechte Hand Seiten (b Vektoren) sind orthogonal zum Nullraum (Nullraum bedeutet Eigenvektoren mit einem Eigenwert von Null). 2) Die Lösung, die Sie erhalten, ist orthogonal zum Nullraum.

Es ist interessant zu bemerken, dass wenn Sie einen nicht-trivialen Nullraum haben, verschiedene numerische Solver Ihnen unterschiedliche Antworten für das gleiche genaue System geben können. Die Lösungen werden sich am Ende durch eine lineare Kombination des Nullraums unterscheiden ... Dieses Problem hat mich viele viele Stunden des Debuggens und der Frustration verursacht, bevor ich mich schließlich durchgesetzt habe, also ist es gut, sich dessen bewusst zu sein.

Schließlich, wenn Ihre Matrix eine Circulant Band structure hat, können Sie einen schnellen Fourier-Transformation (FFT) basierten Löser verwenden. FFT-basierte numerische Solver können in Fällen, in denen sie anwendbar sind, oft eine überlegene Leistung erbringen.

Answer 2

Wenn Sie nicht mit einem Open-Source-Bibliothek geht nichts ausmacht, können Sie auch SPITZE Besuche: CUSP Quick Start Page

Es hat eine recht ordentlich Suite von Solver, darunter ein paar vorkonditioniert Methoden: CUSP Preconditioner Examples

Der geglättete Aggregation Preconditioner (eine Variante des algebraischen Multigrids) scheint sehr gut zu funktionieren, solange Ihre GPU genug internen Speicher dafür hat.