Решение для определения собственных значений параллельно с использованием CUDA

Question

Я искал и искал в Интернете, и я не могу найти ответ, который я ищу. У меня есть особая проблема.

Я редактирую это, чтобы просто проблема и надеюсь, что это более читаемо и понятно.

Предположим, у меня 5000 симметричных плотных матриц 20x20. Я хотел бы создать ядро ​​в CUDA, которое будет иметь каждый поток, ответственный за вычисление собственных значений для каждой из симметричных матриц.

Пример кода ядра CUDA был бы замечательным, если это было возможно.

Любая помощь и предложения будут оценены!

Спасибо,

Johnathan

Answer 1

Я хотел бы создать ядро ​​в CUDA, который будет иметь каждый поток, отвечающий за вычисление собственных значений для каждого из симметричных матриц.

Для меня сомнительно, будет ли это самый быстрый подход, но это может быть для очень маленьких матриц. Даже в этой ситуации могут быть некоторые оптимизации хранения данных, которые могут быть сделаны (перемежение глобальных данных по потокам), но это может усложнить ситуацию.

Как указано, этот запрос может быть сопоставлен с алгоритмом «смущающей параллели», где каждый поток работает с полностью независимой задачей. Нам нужно найти подходящий однопоточный «код-донор». После быстрого поиска в Google я наткнулся на this. Достаточно просто изменить этот код для запуска в этом поточно-независимом режиме. Нам нужно только взять 3 подпрограммы (jacobi_eigenvalue, r8mat_diag_get_vector и r8mat_identity) и соответствующим образом украсить эти подпрограммы __host__ __device__ для использования на графическом процессоре, делая никаких других изменений.

Этот код является GNU LGPL, лицензированным J Burkardt в Университете штата Флорида. Поэтому, с учетом этого, и после conventional wisdom, я не включил значительную часть этого кода в этот ответ. Но вы должны быть в состоянии восстановить мои результаты экспериментально, используя инструкции, которые я даю.

ПРИМЕЧАНИЕ. Я не уверен, какие юридические последствия используют этот код, который утверждает, что он лицензирован GNU LGPL. Вы должны обязательно придерживаться any necessary requirements, если вы решите использовать этот код или его части. Моя основная цель использования здесь заключается в том, чтобы продемонстрировать концепцию относительно тривиального «неудобного параллельного» расширения однопоточного решателя проблем.

Должно быть тривиально восстановить мой полный код, перейдя в here и скопировав 3 указанных функции в места, указанные в остальном скелете кода. Но это не меняет ни одно из вышеупомянутых уведомлений/отказ от ответственности. Используйте его на свой страх и риск.

Опять же, никакие другие изменения не могут быть лучшей идеей с точки зрения производительности, но это приводит к тривиальной сумме усилий и может служить полезной отправной точкой.Некоторые возможные оптимизации могут быть:

1. искать стратегию перемежения данных так, чтобы смежные темы, более вероятно, будет читать соседние данные
2. устранить new и delete функции из кода потока, и заменить его неподвижный распределение (это легко сделать)
3. удалить ненужный код - например, тот, который вычисляет и сортирует собственные векторы, если эти данные не нужны

в любом случае, с вышеуказанным кодом украшен донорского, мы п eed только обертывают вокруг него тривиальное ядро ​​(je), чтобы запускать каждый поток, работающий на отдельных наборах данных (т. матрицы), и каждый поток создает собственный набор собственных значений (и собственных векторов - для этой конкретной кодовой базы).

Я разработал его для работы только с 3 нитями и 3 матрицами 4x4 для тестирования, но это должно быть тривиально, чтобы расширить его до скольких матриц/потоков, которые вы хотите.

Для краткости презентации я отказался от the usual error checking, но я рекомендую вам использовать его или хотя бы запустить код с cuda-memcheck, если вы внесете какие-либо изменения.

Я также построил код для настройки размера кучи устройства вверх для размещения в ядре new операций, в зависимости от количества матриц (то есть потоков) и размеров матрицы. Если вы работали над второй оптимизацией, упомянутой выше, вы, вероятно, могли бы удалить ее.

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; 
} 

компиляции и пример работы:

$ 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 
$ 

Выход кажется правдоподобным для меня, в основном соответствующий выход here.