Ich möchte einen Kernel in CUDA erstellen, die jeden Thread verantwortlich für die Berechnung der Eigenwerte für jede der symmetrischen Matrizen haben.
Es ist fraglich für mich, ob dies der schnellste Ansatz wäre, aber es könnte für sehr kleine Matrizen sein. Sogar in dieser Situation könnten einige Datenspeicheroptimierungen vorgenommen werden (Verschachtelung globaler Daten über Threads hinweg), dies würde jedoch die Dinge komplizieren.
Wie bereits erwähnt, könnte diese Anfrage in einen "peinlich parallelen" Algorithmus abgebildet werden, bei dem jeder Thread nach einem völlig unabhängigen Problem arbeitet. Wir brauchen nur einen geeigneten "Donor Code" mit einem einzigen Gewinde zu finden. Nach einer schnellen Google-Suche stieß ich auf this. Es ist ziemlich einfach, diesen Code so zu ändern, dass er auf diese Thread-unabhängige Weise ausgeführt wird. Wir benötigen nur 3 Routinen (jacobi_eigenvalue, r8mat_diag_get_vector und r8mat_identity), und verzieren diese Routinen mit __host__ __device__ für die Verwendung auf der GPU, während keine anderen Änderungen.
Der fragliche Code scheint GNU LGPL zu sein, lizenziert von J Burkardt an der Florida State University. Daher habe ich in dieser Antwort und nach conventional wisdom keine signifikante Menge dieses Codes in diese Antwort aufgenommen. Aber Sie sollten in der Lage sein, meine Ergebnisse experimentell anhand der Anweisungen zu rekonstruieren, die ich gebe.
HINWEIS: Ich bin mir nicht sicher, welche rechtlichen Konsequenzen es gibt, diesen Code zu verwenden, der laut GNU LGPL lizenziert sein soll. Sie sollten sich an any necessary requirements halten, wenn Sie diesen Code oder Teile davon verwenden. Mein Hauptzweck ist hier, das Konzept einer relativ trivialen "peinlich parallelen" Erweiterung eines single-threaded Problemlösers zu demonstrieren.
Es sollte trivial sein, meinen vollständigen Code zu rekonstruieren, indem Sie here gehen und die 3 angegebenen Funktionen an die Stellen kopieren, die im restlichen Code-Skelett angegeben sind. Dies ändert jedoch nichts an den zuvor genannten Hinweisen/Disclaimern. Verwenden Sie es auf eigene Gefahr.
Auch hier sind keine anderen Änderungen vom Standpunkt der Performance aus gesehen nicht die beste Idee, aber sie führen zu einem trivialen Aufwand und können als möglicher Ausgangspunkt dienen.Einige mögliche Optimierungen könnten sein:
- eine Daten Verschachtelung Strategie suchen, so dass benachbarte Fäden sind eher benachbarten Daten
- beseitigen die
new und delete Funktionen aus dem Thread-Code zu lesen, und ersetzen Sie es mit einem festen Zuordnung (dies ist einfach zu tun)
- Entfernen Sie unnötige Code - zum Beispiel das, was berechnet und die Eigenvektoren sortiert, wenn diese Daten nicht benötigte ist
auf jeden Fall mit dem oben Spender Code eingerichtet, wir n Umschließen Sie einfach einen trivialen Kernel (je), um jeden Thread zu starten, der auf separaten Datensätzen arbeitet (d. h. Matrizen), und jeder Thread erzeugt seine eigene Menge von Eigenwerten (und Eigenvektoren - für diese spezielle Codebasis).
Ich habe es erstellt, um mit nur 3 Threads und 3 4x4 Matrizen für Testzwecke zu arbeiten, aber es sollte trivial sein, es zu beliebig vielen Matrizen/Threads zu erweitern.
Der Kürze der Darstellung halber habe ich auf the usual error checking verzichtet, aber ich empfehle, dass Sie es verwenden oder zumindest Ihren Code mit cuda-memcheck ausführen, wenn Sie irgendwelche Änderungen vornehmen.
Ich habe auch den Code zum Anpassen der Heap-Größe des Geräts nach oben, um die in-Kernel new Operationen, abhängig von der Anzahl der Matrizen (dh Threads) und Matrix-Dimensionen anzupassen. Wenn Sie an der zweiten oben erwähnten Optimierung gearbeitet haben, könnten Sie dies wahrscheinlich entfernen.
t1177.cu:
#include <stdio.h>
#include <iostream>
const int num_mat = 3; // total number of matrices = total number of threads
const int N = 4; // square symmetric matrix dimension
const int nTPB = 256; // threads per block
// test symmetric matrices
double a1[N*N] = {
4.0, -30.0, 60.0, -35.0,
-30.0, 300.0, -675.0, 420.0,
60.0, -675.0, 1620.0, -1050.0,
-35.0, 420.0, -1050.0, 700.0 };
double a2[N*N] = {
4.0, 0.0, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0,
0.0, 0.0, 3.0, 0.0,
0.0, 0.0, 0.0, 2.0 };
double a3[N*N] = {
-2.0, 1.0, 0.0, 0.0,
1.0, -2.0, 1.0, 0.0,
0.0, 1.0, -2.0, 1.0,
0.0, 0.0, 1.0, -2.0 };
/* ---------------------------------------------------------------- */
//
// the following functions come from here:
//
// https://people.sc.fsu.edu/~jburkardt/cpp_src/jacobi_eigenvalue/jacobi_eigenvalue.cpp
//
// attributed to j. burkardt, FSU
// they are unmodified except to add __host__ __device__ decorations
//
//****************************************************************************80
__host__ __device__
void r8mat_diag_get_vector (int n, double a[], double v[])
/* PASTE IN THE CODE HERE, FROM THE ABOVE LINK, FOR THIS FUNCTION */
//****************************************************************************80
__host__ __device__
void r8mat_identity (int n, double a[])
/* PASTE IN THE CODE HERE, FROM THE ABOVE LINK, FOR THIS FUNCTION */
//****************************************************************************80
__host__ __device__
void jacobi_eigenvalue (int n, double a[], int it_max, double v[],
double d[], int &it_num, int &rot_num)
/* PASTE IN THE CODE HERE, FROM THE ABOVE LINK, FOR THIS FUNCTION */
// end of FSU code
/* ---------------------------------------------------------------- */
__global__ void je(int num_matr, int n, double *a, int it_max, double *v, double *d){
int idx = threadIdx.x+blockDim.x*blockIdx.x;
int it_num;
int rot_num;
if (idx < num_matr){
jacobi_eigenvalue(n, a+(idx*n*n), it_max, v+(idx*n*n), d+(idx*n), it_num, rot_num);
}
}
void initialize_matrix(int mat_id, int n, double *mat, double *v){
for (int i = 0; i < n*n; i++) *(v+(mat_id*n*n)+i) = mat[i];
}
void print_vec(int vec_id, int n, double *d){
std::cout << "matrix " << vec_id << " eigenvalues: " << std::endl;
for (int i = 0; i < n; i++) std::cout << i << ": " << *(d+(n*vec_id)+i) << std::endl;
std::cout << std::endl;
}
int main(){
// make sure device heap has enough space for in-kernel new allocations
const int heapsize = num_mat*N*sizeof(double)*2;
const int chunks = heapsize/(8192*1024) + 1;
cudaError_t cudaStatus = cudaDeviceSetLimit(cudaLimitMallocHeapSize, (8192*1024) * chunks);
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "set device heap limit failed!");
}
const int max_iter = 1000;
double *h_a, *d_a, *h_v, *d_v, *h_d, *d_d;
h_a = (double *)malloc(num_mat*N*N*sizeof(double));
h_v = (double *)malloc(num_mat*N*N*sizeof(double));
h_d = (double *)malloc(num_mat* N*sizeof(double));
cudaMalloc(&d_a, num_mat*N*N*sizeof(double));
cudaMalloc(&d_v, num_mat*N*N*sizeof(double));
cudaMalloc(&d_d, num_mat* N*sizeof(double));
memset(h_a, 0, num_mat*N*N*sizeof(double));
memset(h_v, 0, num_mat*N*N*sizeof(double));
memset(h_d, 0, num_mat* N*sizeof(double));
initialize_matrix(0, N, a1, h_a);
initialize_matrix(1, N, a2, h_a);
initialize_matrix(2, N, a3, h_a);
cudaMemcpy(d_a, h_a, num_mat*N*N*sizeof(double), cudaMemcpyHostToDevice);
cudaMemcpy(d_v, h_v, num_mat*N*N*sizeof(double), cudaMemcpyHostToDevice);
cudaMemcpy(d_d, h_d, num_mat* N*sizeof(double), cudaMemcpyHostToDevice);
je<<<(num_mat+nTPB-1)/nTPB, nTPB>>>(num_mat, N, d_a, max_iter, d_v, d_d);
cudaMemcpy(h_d, d_d, num_mat*N*sizeof(double), cudaMemcpyDeviceToHost);
print_vec(0, N, h_d);
print_vec(1, N, h_d);
print_vec(2, N, h_d);
return 0;
}
Kompilierung und Probelauf:
$ nvcc -o t1177 t1177.cu
$ cuda-memcheck ./t1177
========= CUDA-MEMCHECK
matrix 0 eigenvalues:
0: 0.166643
1: 1.47805
2: 37.1015
3: 2585.25
matrix 1 eigenvalues:
0: 1
1: 2
2: 3
3: 4
matrix 2 eigenvalues:
0: -3.61803
1: -2.61803
2: -1.38197
3: -0.381966
========= ERROR SUMMARY: 0 errors
$
Der Ausgang mir plausibel erscheint, vor allem die Ausgabe passend here.
Jacobi ist auch iterativ. Vielleicht ist das Problem Ihre erste Vermutung. Wenn Ihre Schätzung gut ist, benötigen Sie weniger Iterationen, um zu einer Lösung zu konvergieren. Sind Sie nur an Eigenwerten interessiert oder benötigen Sie Eigenvektoren? – duffymo
Ich brauche nur die Eigenwerte. Es gibt Artikel, die über die parallele Durchführung von Jacobi geschrieben wurden. Darauf habe ich oben Bezug genommen. Entschuldigung für die Verwirrung. Um meine Situation besser zu erklären, führe ich vorher andere Berechnungen durch, die sich mit bestimmten Merkmalen jedes Punktes beschäftigen. Ich berechne dann Eigenwerte basierend auf diesen Merkmalen. Der iterative Teil besteht nur darin, diese Eigenwerte zu berechnen. Ich hoffe das hilft. Lassen Sie mich wissen, ob das hilft! – Johnathan
Ich bin mir nicht sicher, ob ich verstehe was "Punkte" bedeutet. Ist das ein Synonym für die Anzahl der Zeilen/Spalten in der quadratischen Matrix? Größere Matrizen bedeuten mehr Berechnungen - nichts wird das ändern. – duffymo