🥑CPU和GPU(C代码)并行计算及二维网格热传输

关键词

CPU | GPU | C | 并行计算 | 二维 | 网格 | CUDA | 热传输 | 线程 | 雅可比迭代 | 图结构 | 半监督学习 | Python | C++

CPU和GPU(C代码)并行计算及二维网格热传输 | Notion亞圖跨際 on Notion

🏈指点迷津 | Brief

要点

1. CPU并行计算示例（底层C代码）

   1. C代码读取命令行参数，垂直或水平翻转输入图像，使用定时函数加速线程执行，多线程创建，在多线程中分割任务和数据

   2. 分析多线程水平翻转图像内存消耗和改善

   3. 创建一个既占用内存又占用内核的旋转图片的C 代码，再创建一个充满零图像的函数，主函数中，旋转图片代码将用户提供的旋转度数转换为弧度，再调用旋转函数旋转图片，主函数中再创建一个空白图像并启动多个线程。创建新代码，量化分析线程效率。

   4. 线程管理和同步：创建一个占内核、占内存、并且由多个独立操作（函数）组成的程序：高斯过滤算法函数，边缘检测内核放大边缘函数，将灰度图像转换为二值（黑白）图像的操作函数。使用预计算减少带宽：预先高斯过滤函数，预先边缘检测函数，预先灰度转换函数。

2. GPU和CUDA并行计算

   1. cuda代码：将图像读入 CPU 端数组，使用Nvidia API 初始化GPU，启用GPU内核，使用cuda Malloc()创建GPU端内存，CPU至GPU数据传输，主函数中GPU完成处理，将数据返回CPU。示例：线性索引读写BMP文件，GPU内核函数垂直和水平翻转图像，内核函数拷贝图像

   2. CUDA 主机/设备编程模型：修改上述cuda代码，显示执行时支持的最大块数。在多块GPU上执行代码观察Kernel性能，观察GPU Core性及其内存，CUDA代码执行观察CPU和GPU之间数据流。

3. CUDA应用示例：模拟二维网格热传输。加速图结构半监督学习。

GPU模拟二维网格热传输示例

我使用了两种类似的 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 行代码），编程部分应该非常简单。 处理这些类型的程序时的主要问题是确定并行化点。