我使用了两种类似的 CUDA(统一计算设备架构)方法。本文讨论了这两种解决方案、它们的优缺点,以及一些可应用于类似并行化问题的理论基础。首先,我们简单讨论一下问题:传热拉普拉斯方程,更具体地说,是它使用雅可比迭代法的解。
其中,它使用以下方程计算热量如何随时间局部传递:
雅可比迭代法是一种通过不断细化解来迭代计算微分方程解的方法,直到答案收敛于稳定解或已完成一些固定数量的步骤,并且答案被认为足够好或不收敛 。示例代码表示材料的 2D 平面,该平面已被划分为大小相等的单元网格。 当热量施加到该平面的中心时,拉普拉斯方程指示热量如何随时间从网格点传递到网格点。
为了计算下一次迭代的给定网格点的温度,需要计算当前迭代中相邻网格点的温度平均值。该方程(C 代码)的主要雅可比迭代循环最终如下:
while ( hasError && iter < iter_max_num ) {
hasError = 0;
for( int i = 1; i < size-1; i++) {
for(int j = 1; j < size-1; j++) {
new_grid[i][j] = 0.25 * (grid[i][j+1] + grid[i][j-1] +
grid[i-1][j] + grid[i+1][j]);
double errorCurrent = fabs(new_grid[i][j] - grid[i][j]);
if(errorCurrent > CONV_THRESHOLD ){
hasError = 1;
}
}
}
aux = grid;
grid = new_grid;
new_grid = aux;
iter++;
}
完整的解可以在下面的要点中找到:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
double **grid;
double **new_grid;
double **aux;
int size;
int iter_max_num;
void allocate_memory(){
grid = (double **) malloc(size * sizeof(double *));
new_grid = (double **) malloc(size * sizeof(double *));
for(int i = 0; i < size; i++){
grid[i] = (double *) malloc(size * sizeof(double));
new_grid[i] = (double *) malloc(size * sizeof(double));
}
}
void initialize_grid(){
int center = size / 2;
int linf = center - (size / 10);
int lsup = center + (size / 10);
for (int i = 0; i < size; i++){
for (int j = 0; j < size; j++){
if ( i >= linf && i <= lsup && j >= linf && j <= lsup)
grid[i][j] = 100;
else
grid[i][j] = 0;
new_grid[i][j] = 0.0;
}
}
}
void save_grid(){
for(int i = 0; i < size; i++){
for(int j = 0; j < size; j++){
fprintf(stdout, "%lf ", grid[i][j]);
}
fprintf(stdout, "\\n");
}
}
int main(int argc, char *argv[]){
if (fscanf(stdin, "%d", &size) != 1) {
printf("Error: could not read input Size.\\n");
exit(1);
}
if (fscanf(stdin, "%d", &iter_max_num) != 1) {
printf("Error: could not read Iter_Max_Num Size.\\n");
exit(1);
}
int hasError = 1;
int iter = 0;
struct timeval time_start;
struct timeval time_end;
allocate_memory();
initialize_grid();
printf("Jacobi relaxation calculation: %d x %d grid\\n", size, size);
gettimeofday(&time_start, NULL);
while ( hasError && iter < iter_max_num ) {
hasError = 0;
for( int i = 1; i < size-1; i++) {
for(int j = 1; j < size-1; j++) {
new_grid[i][j] = 0.25 * (grid[i][j+1] + grid[i][j-1] +
grid[i-1][j] + grid[i+1][j]);
double errorCurrent = fabs(new_grid[i][j] - grid[i][j]);
if(errorCurrent > CONV_THRESHOLD ){
hasError = 1;
}
}
}
aux = grid;
grid = new_grid;
new_grid = aux;
iter++;
}
gettimeofday(&time_end, NULL);
double exec_time = (double) (time_end.tv_sec - time_start.tv_sec) +
(double) (time_end.tv_usec - time_start.tv_usec) / 1000000.0;
save_grid();
printf("\\nKernel executed in %lf seconds with %d iterations.\\n", exec_time, iter);
return 0;
}
对于 2048 x 2048 矩阵和 15000 次迭代的收敛阈值,该算法大约需要 2.43 秒(是的,这很多)。这意味着这是不可扩展的,并指向可能的并行优化。有很多选项可以使程序并行。 可以轻松应用于 C 程序的一个简单解决方案是使用 OpenMP(开放式多处理)API 。
但如果我们想要更进一步,我们可以利用 GPU 来实现更大的加速。 鉴于程序整体较小(总共只有 140 行代码),编程部分应该非常简单。 处理这些类型的程序时的主要问题是确定并行化点。
https://embed.notionlytics.com/wt/ZXlKM2IzSnJjM0JoWTJWVWNtRmphMlZ5U1dRaU9pSlhiRWhvWlV4VVQxbHNjMlZYV2tKbU9URndaU0lzSW5CaFoyVkpaQ0k2SWpBMU9EVXlOemRtTTJGa016UmlaREU0WTJObE9XVXlOek5sWmpNNVpUVmpJbjA9