0

我一直在使用传热代码。基本上,这段代码确定了立方体及其所有面的初始条件。这六个面从不同的温度开始,然后代码将计算由于它们之间的热传递而导致所有面的温度如何变化。现在,我一直在尝试使用 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;
}
4

1 回答 1

1

使用具有动态存储的可变长度数组:

  1. 分配:
double (*arr)[N][N] = calloc(N, sizeof *arr);
  1. 索引。使用良好的旧arr[i][j][k]语法

  2. 解除分配。

free(arr)
  1. 展平。
double *flat = (double*)arr;

请注意, C 标准保证此转换可以正常工作。尽管它很可能适用于所有能够使用 GPU 的平台。

  1. 传递给函数。VLA 可以是函数的参数。
void fun(int n, double arr[n][n][n]) {
  ...
}

示例性用法是:

foo(N, arr);

编辑

VLA 友好的变体compute_heat_transfer()

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

关键字restrictinarrNew[restrict n][n][n]用于让编译器假设arrNewarrOld不是别名。它应该让编译器使用更积极的优化。

请注意,arrNewandarrOld指向数组的指针。因此,与其复制arrNewarrOld您,不如简单地交换这些指针,形成一种简单的双缓冲。它应该使代码更快。

于 2021-09-28T18:46:18.570 回答