使用 OpenMP 执行三重指针 (C) 卸载到 NVIDIA GPU

Question

我一直在研究传热代码。基本上，这段代码确定了立方体及其所有面的初始条件。六个面从不同的温度开始，然后代码将计算由于它们之间的热传递，所有面的温度如何变化。现在，我一直在尝试使用 OpenMP 指令卸载到 NVIDIA GPU。此代码使用三重指针初始化面部条件，这是一种数组数组。阅读了一些有关此事的资料后，我了解到 3D 架构不容易卸载到 GPU。所以我的问题是，是否可以将这个三重指针数组卸载到 GPU，或者我是否必须使用更扁平的数组形式。

这里我留下代码，它仍在 CPU 上工作。代码的并行版本。

#include <omp.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define N 25 //This defines the number of points per dimension (Cube = N*N*N)
#define NUM_STEPS 6000 //This is the number of simulations time steps

/*writeFile: this function writes simulation results into a file. 
 * A file is created for each iteration that's passed to the function 
 * as a parameter. It also takes the triple pointer to the simulation 
 * data*/
void writeFile(int iteration, double*** data){
    char filename[50];
    char itr[12];
    sprintf(itr, "%d", iteration);

    strcpy(filename, "heat_");
    strcat(filename, itr);
    strcat(filename, ".txt");
    //printf("Filename is %s\n", filename);

    FILE *fp;
    fp = fopen(filename, "w");

    fprintf(fp, "x,y,z,T\n");
    for(int i=0; i<N; i++){
        for(int j=0;j<N; j++){
            for(int k=0; k<N; k++){
                fprintf(fp,"%d,%d,%d,%f\n", i,j,k,data[i][j][k]);
            }
        }
    }
    fclose(fp);
}
void compute_heat_transfer(double ***arrayOld, double ***arrayNew){
    int i,j,k;
    /*Compute steady-state solution*/
    for(int nsteps=0; nsteps < NUM_STEPS; nsteps++){
        /*if(nsteps % 100 == 0){
            writeFile(nsteps, arrayOld);
        }*/
        #pragma omp parallel shared(arrayNew, arrayOld) private(i,j,k)
        {
            #pragma omp for
            for(i=1; i<N-1; i++){
                for(j=1; j<N-1; j++){
                    for(k=1;k<N-1;k++){
                        //This is the 6-neighbor stencil computation
                        arrayNew[i][j][k] = (arrayOld[i-1][j][k] + arrayOld[i+1][j][k] + arrayOld[i][j-1][k] + arrayOld[i][j+1][k] +
                                 arrayOld[i][j][k-1] + arrayOld[i][j][k+1])/6.0;
                    }
                }
            }
            #pragma omp for
            for(i=1; i<N-1; i++){
                for(j=1; j<N-1; j++){
                    for(k=1; k<N-1; k++){
                        arrayOld[i][j][k] = arrayNew[i][j][k];
                    }   
                }
            }
        }
    }
}

int main (int argc, char *argv[]) {
    int i,j,k,nsteps; 
    double mean;
    double ***arrayOld; //Variable that will hold the data of the past iteration
    double ***arrayNew; //Variable where newly computed data will be stored
    arrayOld = (double***)malloc(N*sizeof(double**));
    arrayNew = (double***)malloc(N*sizeof(double**));
    if(arrayOld== NULL){
        fprintf(stderr, "Out of memory");
        exit(0);
    }
    for(i=0; i<N;i++){
        arrayOld[i] = (double**)malloc(N*sizeof(double*));
        arrayNew[i] = (double**)malloc(N*sizeof(double*));
        if(arrayOld[i]==NULL){
            fprintf(stderr, "Out of memory");
            exit(0);
        }
        for(int j=0;j<N;j++){
            arrayOld[i][j] = (double*)malloc(N*sizeof(double));
            arrayNew[i][j] = (double*)malloc(N*sizeof(double));
            if(arrayOld[i][j]==NULL){
                fprintf(stderr,"Out of memory");
                exit(0);
            }
        }
    }

    /*Set boundary values and compute mean boundary values*/
    mean = 0.0;

    for(i=0; i<N; i++){
        for(j=0;j<N;j++){
            arrayOld[i][j][0] = 100.0;
            mean += arrayOld[i][j][0];
        }
    }

    for(i=0; i<N; i++){
        for(j=0;j<N;j++){
            arrayOld[i][j][N-1] = 100.0;
            mean += arrayOld[i][j][N-1];
        }
    }

    for(j=0; j<N; j++){
        for(k=0;k<N;k++){
            arrayOld[0][j][k] = 100.0;
            mean += arrayOld[0][j][k];
        }
    }
    
    for(j=0; j<N; j++){
        for(k=0;k<N;k++){
            arrayOld[N-1][j][k] = 100.0;
            mean += arrayOld[N-1][j][k];
        }
    }

    for(i=0; i<N; i++){
        for(k=0;k<N;k++){
            arrayOld[i][0][k] = 100.0;
            mean += arrayOld[i][0][k];
        }
    }
    
    for(i=0; i<N; i++){
        for(k=0;k<N;k++){
            arrayOld[i][N-1][k] = 0.0;
            mean += arrayOld[i][N-1][k];
        }
    }

    mean /= (6.0 * (N*N));

    /*Initialize interior values*/
    for(i=1; i<N-1; i++){
        for(j=1; j<N-1; j++){
            for(k=1; k<N-1;k++){
                arrayOld[i][j][k] = mean;
            }
        }
    }

    double tdata = omp_get_wtime();

    compute_heat_transfer(arrayOld, arrayNew);

    tdata = omp_get_wtime()-tdata;

    printf("Execution time was %f secs\n", tdata);
    
    for(i=0; i<N;i++){
        for(int j=0;j<N;j++){
            free(arrayOld[i][j]);
            free(arrayNew[i][j]);
        }
        free(arrayOld[i]);
        free(arrayNew[i]);
    }
    free(arrayOld);
    free(arrayNew);
    
    return 0;
}

Answer 1

使用动态存储的可变长度数组：

分配：

double (*arr)[N][N] = calloc(N, sizeof *arr);

正在编制索引。使用很好的旧 arr[i][j][k] 语法
取消分配。

free(arr)

扁平化。

double *flat = (double*)arr;

请注意，C 标准不保证此转换有效。尽管它很可能适用于所有能够使用 GPU 的平台。

传递给函数。 VLA 可以是函数的参数。

void fun(int n, double arr[n][n][n]) {
  ...
}

示例用法为：

foo(N, arr);

编辑

compute_heat_transfer() 的 VLA 友好变体：

void compute_heat_transfer(int n, double arrayOld[restrict n][n][n], double arrayNew[restrict n][n][n]) {
    int i,j,k;
    /*Compute steady-state solution*/
    for(int nsteps=0; nsteps < NUM_STEPS; nsteps++){
        /*if(nsteps % 100 == 0){
            writeFile(nsteps, arrayOld);
        }*/
        #pragma omp parallel for collapse(3) 
        for(i=1; i<n-1; i++){
        for(j=1; j<n-1; j++){
        for(k=1; k<n-1; k++){
           //This is the 6-neighbor stencil computation
          arrayNew[i][j][k] = (arrayOld[i-1][j][k] + arrayOld[i+1][j][k] + arrayOld[i][j-1][k] + arrayOld[i][j+1][k] +
                                 arrayOld[i][j][k-1] + arrayOld[i][j][k+1])/6.0;
        }}}
        #pragma omp parallel for collapse(3)
        for(i=1; i<n-1; i++){
        for(j=1; j<n-1; j++){
        for(k=1; k<n-1; k++){
          arrayOld[i][j][k] = arrayNew[i][j][k];
        }}}
    }
}

arrNew[restrict n][n][n] 中的关键字 restrict 用于让编译器假设 arrNew 和 arrOld 没有别名。它应该让编译器使用更积极的优化。

请注意，arrNew 和 arrOld 是指向数组的指针。因此，与其将 arrNew 复制到 arrOld，您还可以简单地交换这些指针，形成一种简单的双缓冲。它应该使代码更快。

使用 OpenMP 执行三重指针 (C) 卸载到 NVIDIA GPU

Perform a triple pointer (C) offloading to NVIDIA GPU with OpenMP

c

gcc

gpu

openmp

offloading