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

要点

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应用示例：模拟二维网格热传输。加速图结构半监督学习。

参阅其他 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 行代码），编程部分应该非常简单。 处理这些类型的程序时的主要问题是确定并行化点。

参阅一：计算思维

参阅二：亚图跨际