c - 编写混合 MPI/OpenACC 程序-6ren

c - 编写混合 MPI/OpenACC 程序

转载作者：太空宇宙更新时间：2023-11-04 06:21:58

我正在尝试编写混合 MPI/OpenACC 代码，其中代码需要执行 8 项不同的工作(在本例中为 8 项不同的扫描)。这 8 个作业使用 MPI 划分为 [1-8] 个进程/节点，并使用 OpenACC 并行化这 8 个作业中需要完成的计算。

在每个进程完成其计算后，我减少解决方案并将最小值传递给进程 0，这是最终解决方案。

下面是生成 .txt 输出文件的完整代码 (test.c) 的 MCVE

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "mpi.h"

#define min(a,b) (a > b) ? b : a
#define max(a,b) (a < b) ? b : a

#define NPES 8 // max number of PEs allowed
#define DEFAULT_BORDER_LOCATION   -1
#define DEFAULT_BORDER_DISTANCE   INFINITY
#define DEFAULT_INTERIOR_DISTANCE 90000

typedef struct {
  int order;
  int firstLevel, lastLevel, level;
  int xDim, yDim, zDim;
  int xSweepOff, ySweepOff, zSweepOff;
  double dx, dy, dz;
} SweepInfo;

typedef struct {
  double   dx, dy, dz;
  int * location;
  double * distance;
} Phi;

typedef struct {
  int x, y, z;
} Grid3D;



void calc_dist_field( Phi * p, int totalNodes );
void write_to_file(double * dist);
static SweepInfo make_sweepInfo( Phi * p, int my_rank );
static void fast_sweep( Phi * p, SweepInfo * s );
static double solveEikonal(Phi * p, int index, int max_x, int max_y);
static void update_distance(Phi * p, int totalNodes);
static void set_distance_negative_inside(Phi * p, int totalNodes);
static void adjust_boundary( Phi * p );


// public method declarations
Grid3D make_grid3D(int x, int y, int z);
void vti_get_dimensions(FILE *vti, double *d);
void vti_get_data(FILE *vti, int *l, int b_l, double *d, double b_d, Grid3D g);

// private method declarations
static void move_file_pointer(FILE *file_ptr, int lineNumber, int r);
static void get_location(FILE *vti, int *l, int b_l, Grid3D g);
static void get_distance(FILE *vti, double *d, double b_d, Grid3D g);

static int npes;          // Number of PEs
static int my_rank;       // Rank of the PE
static char * fileName;
static char * outfileName;

static int NX, NY, NZ, totalNodes;

int main(int argc, char *argv[]) {

  // MPI startup routine
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &npes);

  fileName    = argv[1];
  outfileName = argv[2];
  FILE *f = fopen(fileName, "r");

  double dims[6];
  vti_get_dimensions(f, dims);
  NX = dims[0] + 3;
  NY = dims[1] + 3;
  NZ = dims[2] + 3;

  totalNodes = NX * NY * NZ;
  Phi *p = (Phi *) malloc(sizeof(Phi));
  p->location = (int *) malloc(sizeof(int) * totalNodes);
  p->distance = (double *) malloc(sizeof(double) * totalNodes);

  p->dx = dims[3]; p->dy = dims[4]; p->dz = dims[5];

  vti_get_data( f, p->location, DEFAULT_BORDER_LOCATION,
                   p->distance, DEFAULT_BORDER_DISTANCE,
                make_grid3D(NX, NY, NZ));


  update_distance(p, totalNodes);

  calc_dist_field(p, totalNodes);

  MPI_Finalize();
  return 0;
}

void calc_dist_field( Phi * p, int totalNodes ) {

  int sweepNumber = my_rank + 1;
  double * tmp_dist;

  MPI_Barrier(MPI_COMM_WORLD);
  if(my_rank == 0){
    tmp_dist = (double *) malloc( totalNodes * sizeof(double) );
  }

  // sn represents the sweep number
  for( int sn = sweepNumber; sn <= NPES; sn += npes) {
    SweepInfo s = make_sweepInfo(p, sn);

    printf("PE: [%d] - performing sweep number ..... [%d/%d]\n", my_rank, sn, NPES);

    fast_sweep(p, &s);

    printf("PE: [%d] - completed sweep number ...... [%d/%d]\n", my_rank, sn, NPES);

  }

  MPI_Barrier(MPI_COMM_WORLD);
  #pragma acc update host(p->distance[0:totalNodes])
  MPI_Reduce(p->distance, tmp_dist, totalNodes, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);

  if( my_rank == 0 ) {
    free( p->distance );
    p->distance = tmp_dist;
    set_distance_negative_inside(p, totalNodes);
    adjust_boundary(p);
    write_to_file(p->distance);
    printf("%s file created\n", outfileName);
  }

}

static void update_distance(Phi * p, int totalNodes) {

  int    *l = &p->location[0];
  double *d = &p->distance[0];

  for(int i = 0; i < totalNodes; i++) {
    if(*l != DEFAULT_BORDER_LOCATION && *d != DEFAULT_BORDER_DISTANCE ) {
      *d = (*l == 1 && *d == INFINITY) ? -1 : (*d > 0.0 || *d < 0.0) ? *d : DEFAULT_INTERIOR_DISTANCE;
    }
    l++; d++;
  }

}



void write_to_file(double * dist) {

  int x = NX;
  int y = NY;
  int z = NZ;
  char fname[255];
  sprintf(fname, "%s.txt", outfileName);
  FILE *fp = fopen(fname, "w");

  int i,j,k;
  double *t = &dist[0];
  for(i = 0; i < z; i++){
    for(j = 0; j < y; j++){
      for(k = 0; k < x; k++) {
    fprintf(fp, "%f  ", *(t++));
      }
      fprintf(fp, "\n");
    }
    fprintf(fp, "\n");
  }

}


static SweepInfo make_sweepInfo( Phi * p, int my_rank ) {
  SweepInfo s;

  s.order      = my_rank;
  s.firstLevel = 3;
  s.lastLevel  = (NX + NY + NZ) - 6;

  s.xDim = NX-2; s.dx = p->dx;
  s.yDim = NY-2; s.dy = p->dy;
  s.zDim = NZ-2; s.dz = p->dz;

  s.xSweepOff = (s.order == 4 || s.order == 8 ) ? s.xDim + 1 : 0;
  s.ySweepOff = (s.order == 2 || s.order == 6 ) ? s.yDim + 1 : 0;
  s.zSweepOff = (s.order == 3 || s.order == 7 ) ? s.zDim + 1 : 0;

  return s;
}

static void fast_sweep( Phi * p, SweepInfo * s ) {

  int start, end, incr;

  start = ( s->order == 2 || s->order == 5 || s->order == 7 || s->order == 8 ) ? s->lastLevel : s->firstLevel;

  if ( start == s->firstLevel ) {
    end  = s->lastLevel + 1;
    incr = 1;
  }
  else {
    end  = s->firstLevel - 1;
    incr = 0;
  }

  int max_x  = s->xDim + 2;
  int max_y  = s->yDim + 2;
  int max_xy = max_x * max_y;

  #pragma acc data create(p[0:1]) copy(p->distance[0:totalNodes])
  for(int level = start; level != end; level = (incr) ? level+1 : level-1) {
    // s - start, e - end
    int xs, xe, ys, ye;

    xs = max(1, level-(s->yDim + s->zDim))    , ys = max(1,level-(s->xDim + s->zDim));
    xe = min(s->xDim, level-(s->firstLevel-1)), ye = min(s->yDim, level-(s->firstLevel-1));

    int x, y, z, i, j, k, index;
    #pragma acc parallel
    {
      #pragma acc loop independent
      for(x = xs; x <= xe; x++) {
        #pragma acc loop independent
        for(y = ys; y <= ye; y++) {
          z = level - (x+y);
          if(z > 0 && z <= NZ-2) {
            i = abs(z-s->zSweepOff);
            j = abs(y-s->ySweepOff);
            k = abs(x-s->xSweepOff);
            index = i * max_xy + j * max_x + k;
            p->distance[index] = solveEikonal(p, index, NX, NY);
          }
        }
      } // end of acc parallel
    }
  }
}

#pragma acc routine seq
static double solveEikonal(Phi * p, int index, int max_x, int max_y) {

  int max_xy = max_x * max_y;

  double dist_new = 0;
  double dist_old = p->distance[index];

  double dx = p->dx, dy = p->dy, dz = p->dz;
  double minX = min(p->distance[index-1], p->distance[index+1]);
  double minY = min(p->distance[abs(index-max_x)], p->distance[abs(index+max_x)]);
  double minZ = min(p->distance[abs(index-max_xy)],p->distance[abs(index+max_xy)]);

  double m[] = { minX, minY, minZ} ;
  double d[] = { dx, dy, dz};

  // sort the mins 
  for(int i = 1; i < 3; i++){
    for(int j = 0; j < 3-i; j++) {
      if(m[j] > m[j+1]) {
        double tmp_m = m[j];
        double tmp_d = d[j];
        m[j] = m[j+1]; d[j] = d[j+1];
        m[j+1] = tmp_m; d[j+1] = tmp_d;
      }
    }
  }

  // simplifying the variables
  double m_0 = m[0], m_1 = m[1], m_2 = m[2];
  double d_0 = d[0], d_1 = d[1], d_2 = d[2]; 
  double m2_0 = m_0 * m_0, m2_1 = m_1 * m_1, m2_2 = m_2 * m_2;
  double d2_0 = d_0 * d_0, d2_1 = d_1 * d_1, d2_2 = d_2 * d_2;

  dist_new = m_0 + d_0;
  if(dist_new > m_1) {

    double s = sqrt(- m2_0 + 2 * m_0 * m_1 - m2_1 + d2_0 + d2_1); 
    dist_new = ( m_1 * d2_0 + m_0 * d2_1 + d_0 * d_1 * s) / (d2_0 + d2_1);

    if(dist_new > m_2) {

      double a = sqrt(- m2_0 * d2_1 - m2_0 * d2_2 + 2 * m_0 * m_1 * d2_2
                      - m2_1 * d2_0 - m2_1 * d2_2 + 2 * m_0 * m_2 * d2_1
                      - m2_2 * d2_0 - m2_2 * d2_1 + 2 * m_1 * m_2 * d2_0
                      + d2_0 * d2_1 + d2_0 * d2_2 + d2_1 * d2_2);

      dist_new = (m_2 * d2_0 * d2_1 + m_1 * d2_0 * d2_2 + m_0 * d2_1 * d2_2 + d_0 * d_1 * d_2 * a) /
                  (d2_0 * d2_1 + d2_0 * d2_2 + d2_1 * d2_2);
    }
  }

  return min(dist_old, dist_new);
}



static void set_distance_negative_inside(Phi * p, int totalNodes) {

  int    *l = &p->location[0];
  double *d = &p->distance[0];

  for(int i = 0; i < totalNodes; i++) {
    if(*l != DEFAULT_BORDER_LOCATION && *d != DEFAULT_BORDER_DISTANCE ) {
      if( *l == 1) *d = -1;
    }
    l++; d++;
  }

}

static void adjust_boundary( Phi * p ) {

  int x, y, z, xy, i, j, k;
  x  = NX;
  y  = NY;
  z  = NZ;
  xy = x * y;

  for(i = 0; i < z; i++){
    for(j = 0; j < y; j++){
      for(k = 0; k < x; k++){
        int I = i, J = j, K = k;
        I = (i == z-1) ? I-1 : (!i) ? I+1 : I;
        J = (j == y-1) ? J-1 : (!j) ? J+1 : J;
        K = (k == x-1) ? K-1 : (!k) ? K+1 : K;
        if( i != I || j != J || k != K) {
            int l_index = i * xy + j * x + k;
            int r_index = I * xy + J * x + K;
            p->distance[l_index] = p->distance[r_index];
        }
      }
    }
  }
}



/**************** vti_parser ********************************/

static void move_file_pointer(FILE *file_ptr, int lineNumber, int r) {
    char tmpStr[512];
    if(r) rewind(file_ptr);
    while (lineNumber > 0){
        fgets (tmpStr, 511, file_ptr);
        lineNumber--;
    }
}

void vti_get_dimensions(FILE *vti, double *d) {
    char tmpStr[512];
    rewind(vti);
    while (1) {
        fgets (tmpStr, 511, vti);
        if ( strstr(tmpStr, "ImageData WholeExtent") ) {
            sscanf(tmpStr, "    <ImageData WholeExtent=\"0 %lf 0 %lf 0 %lf\" Spacing=\"%lf %lf %lf\">",
                    &d[0], &d[1], &d[2], &d[3], &d[4], &d[5]);
            break;
        }
    }
}


void vti_get_data(FILE *vti, int *l, int b_l, double *d, double b_d, Grid3D g) {

    // move the file pointer to
    // line 6 from beginning
    move_file_pointer(vti, 6, 1);

    get_location(vti, l, b_l, g);

    // move the file pointer 2 lines
    // forward from its last position
    move_file_pointer(vti, 2, 0);

    get_distance(vti, d, b_d, g);

}

static void get_location(FILE *vti, int *l, int b_l, Grid3D g) {
    int i, j, k, *t = &l[0];
    for (i = 0; i < g.z; i++){
        for (j = 0; j < g.y; j++) {
            for (k = 0; k < g.x; k++) {
                // Border
                if (k == 0 || k == g.x-1 || j == 0 || j == g.y-1 || i == 0 || i == g.z-1 ) {
                    *(t++) = b_l;
                }
                else{ // Interior
                    fscanf(vti, "%d ", t++);
                }
            }
        }
    }
}



static void get_distance(FILE *vti, double *d, double b_d, Grid3D g) {
    int i, j, k;
    double *t = &d[0];
    for (i = 0; i < g.z; i++){
        for (j = 0; j < g.y; j++) {
            for (k = 0; k < g.x; k++) {
                // Border distance
                if (k == 0 || k == g.x-1 || j == 0 || j == g.y-1 || i == 0 || i == g.z-1 ) {
                    *(t++) = b_d;
                }
                else{ // Interior distance
                    fscanf(vti, "%lf ", t++);
                }
            }
        }
    }
}

Grid3D make_grid3D(int x, int y, int z){
    Grid3D g;
    g.x = x; g.y = y; g.z = z;

    return g;
}

当我放弃 openacc 指令并使用 [1-8] 进程运行它时，代码可以正常工作，但是当使用 open acc 编译器时，我得到一个 cudaError。

调用 cuStreamSynchronize 返回错误 700:内核执行期间地址非法

MPI 编译:

mpicc -Wall -g -std=c99 -I/cm/shared/apps/openmpi/gcc/64/1.8.5_wocuda/include -L/cm/shared/apps/openmpi/gcc/64/1.8.5_wocuda/lib -lmpi test.c -o mpi_exec.out

OpenACC 编译:

pgcc -acc -ta=tesla:managed -Minfo=accel -g -lm -I/cm/shared/apps/openmpi/gcc/64/1.8.5_wocuda/include -L/cm/shared/apps/openmpi/gcc/64/1.8.5_wocuda/lib -lmpi test.c -o oacc_exec.out

要运行可执行文件，您需要传入一个输入 vti 文件和输出文件名。

mpirun -np <1-8> <executable> input.vti outputName

输入文件的链接 input.vti

我希望这段代码非常灵活，我想让它可以在具有 1 个 GPU 的单个节点上运行，同时运行 [1-8] 个进程，也可以在 [1-8] 个节点上运行，每个节点具有[1-2] 图形处理器。而且我没有使用 CUDA MPS。

我的规范

GNU/Linux x86_64
NVIDIA GeForce GTX Titan CC: 3.5

pgcc 15.7-0 64-bit target on x86-64 Linux -tp sandybridge 
gcc (GCC) 4.8.1

对此有任何帮助或建议，我们将不胜感激。

编辑:

** Compiling with OpenACC

`$ pgcc -fast -ta=tesla:managed -Minfo=accel -   I/cm/shared/apps/openmpi/gcc/64/1.8.5_wocuda/include -L/cm/shared/apps/openmpi/gcc/64/1.8.5_wocuda/lib -lmpi rcrovella.c -o withacc
PGC-W-0129-Floating point overflow. Check constants and constant expressions (rcrovella.c: 88)
PGC-W-0129-Floating point overflow. Check constants and constant expressions (rcrovella.c: 142)
PGC-W-0129-Floating point overflow. Check constants and constant expressions (rcrovella.c: 143)
PGC-W-0129-Floating point overflow. Check constants and constant expressions (rcrovella.c: 308)
fast_sweep:
    225, Generating copy(p[:1])
    228, Loop is parallelizable
    230, Loop is parallelizable
         Accelerator kernel generated
         Generating Tesla code
        228, #pragma acc loop gang /* blockIdx.y */
        230, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */
solveEikonal:
    246, Generating acc routine seq
    262, Loop is parallelizable
    263, Loop carried dependence of m prevents parallelization
         Loop carried backward dependence of m prevents vectorization
         Loop carried dependence of d prevents parallelization
         Loop carried backward dependence of d prevents vectorization
PGC/x86-64 Linux 15.7-0: compilation completed with warnings`

** Compiling without OpenACC
pgcc -I/cm/shared/apps/openmpi/gcc/64/1.8.5_wocuda/include -L/cm/shared/apps/openmpi/gcc/64/1.8.5_wocuda/lib -lmpi rcrovella.c -o   noaccPGC-W-0129-Floating point overflow. Check constants and constant    expressions (rcrovella.c: 88)
PGC-W-0129-Floating point overflow. Check constants and constant expressions (rcrovella.c: 142)
PGC-W-0129-Floating point overflow. Check constants and constant expressions (rcrovella.c: 143)
PGC-W-0129-Floating point overflow. Check constants and constant expressions (rcrovella.c: 308)
PGC/x86-64 Linux 15.7-0: compilation completed with warnings

** Running with OpenACC
$ mpirun -n 1 withacc ../my_test/input.vti withacc1
PE: [0] - performing sweep number ..... [1/8]
PE: [0] - completed sweep number ...... [1/8]
PE: [0] - performing sweep number ..... [2/8]
PE: [0] - completed sweep number ...... [2/8]
PE: [0] - performing sweep number ..... [3/8]
PE: [0] - completed sweep number ...... [3/8]
PE: [0] - performing sweep number ..... [4/8]
PE: [0] - completed sweep number ...... [4/8]
PE: [0] - performing sweep number ..... [5/8]
PE: [0] - completed sweep number ...... [5/8]
PE: [0] - performing sweep number ..... [6/8]
PE: [0] - completed sweep number ...... [6/8]
PE: [0] - performing sweep number ..... [7/8]
PE: [0] - completed sweep number ...... [7/8]
PE: [0] - performing sweep number ..... [8/8]
PE: [0] - completed sweep number ...... [8/8]
withacc1 file created

** Running without OpenACC
$ mpirun -n 1 noacc ../my_test/input.vti noacc1
PE: [0] - performing sweep number ..... [1/8]
PE: [0] - completed sweep number ...... [1/8]
PE: [0] - performing sweep number ..... [2/8]
PE: [0] - completed sweep number ...... [2/8]
PE: [0] - performing sweep number ..... [3/8]
PE: [0] - completed sweep number ...... [3/8]
PE: [0] - performing sweep number ..... [4/8]
PE: [0] - completed sweep number ...... [4/8]
PE: [0] - performing sweep number ..... [5/8]
PE: [0] - completed sweep number ...... [5/8]
PE: [0] - performing sweep number ..... [6/8]
PE: [0] - completed sweep number ...... [6/8]
PE: [0] - performing sweep number ..... [7/8]
PE: [0] - completed sweep number ...... [7/8]
PE: [0] - performing sweep number ..... [8/8]
PE: [0] - completed sweep number ...... [8/8]
noacc1 file created

** Compare
$ diff -q noacc1.txt withacc1.txt
Files noacc1.txt and withacc1.txt differ

最佳答案

Also in this version I couldn't get openacc to work at all however solution to this should help me a lot.

这是我发现的:

当您使用托管内存设施时:
```
-ta=tesla:managed
```
通常我们不在代码中包含data 指令或子句。这个想法是让 cuda 管理的内存运行时为我们管理数据移动。所以我注释掉了我认为是两个“无关的”数据指令。
我认为您的parallel 加速器指令格式不正确。我的编译器 (PGI 15.7) 告诉我并行区域中的 independent 指令不正确:
```
PGCC-S-0155-Illegal context(parallel) for independent  (t2.c: 228)
```
将 #pragma acc parallel 更改为 #pragma acc kernels 是一种可能的解决方法。
您的代码在使用 INFINITY 时发出一些编译器警告。由于这些只是警告，我没有费心去解决它们。

出于某种原因，我发现编译器在进入加速器区域时没有正确处理 SweepInfo 结构 (s)。为了解决这个问题，我修改了这个:

int x, y, z, i, j, k, index;
#pragma acc parallel
{
  #pragma acc loop independent
  for(x = xs; x <= xe; x++) {
    #pragma acc loop independent
    for(y = ys; y <= ye; y++) {
      z = level - (x+y);
      if(z > 0 && z <= NZ-2) {
        i = abs(z-s->zSweepOff);
        j = abs(y-s->ySweepOff);
        k = abs(x-s->xSweepOff);

为此:

int x, y, z, i, j, k, index;
int xSO = s->xSweepOff;
int ySO = s->ySweepOff;
int zSO = s->zSweepOff;
#pragma acc kernels
{
  #pragma acc loop independent
  for(x = xs; x <= xe; x++) {
    #pragma acc loop independent
    for(y = ys; y <= ye; y++) {
      z = level - (x+y);
      if(z > 0 && z <= NZ-2) {
        i = abs(z-zSO);
        j = abs(y-ySO);
        k = abs(x-xSO);

我可能会再仔细研究一下。我认为这里要么存在我不理解的限制，要么存在编译器错误。

通过上述更改，我能够让您的代码运行完成，没有任何明显的问题。这是我修改后的代码:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "mpi.h"

#define min(a,b) (a > b) ? b : a
#define max(a,b) (a < b) ? b : a

#define NPES 8 // max number of PEs allowed
#define DEFAULT_BORDER_LOCATION   -1
#define DEFAULT_BORDER_DISTANCE   INFINITY
#define DEFAULT_INTERIOR_DISTANCE 90000

typedef struct {
  int order;
  int firstLevel, lastLevel, level;
  int xDim, yDim, zDim;
  int xSweepOff, ySweepOff, zSweepOff;
  double dx, dy, dz;
} SweepInfo;

typedef struct {
  double   dx, dy, dz;
  int * location;
  double * distance;
} Phi;

typedef struct {
  int x, y, z;
} Grid3D;



void calc_dist_field( Phi * p, int totalNodes );
void write_to_file(double * dist);
static SweepInfo make_sweepInfo( Phi * p, int my_rank );
static void fast_sweep( Phi * p, SweepInfo * s );
static double solveEikonal(Phi * p, int index, int max_x, int max_y);
static void update_distance(Phi * p, int totalNodes);
static void set_distance_negative_inside(Phi * p, int totalNodes);
static void adjust_boundary( Phi * p );


// public method declarations
Grid3D make_grid3D(int x, int y, int z);
void vti_get_dimensions(FILE *vti, double *d);
void vti_get_data(FILE *vti, int *l, int b_l, double *d, double b_d, Grid3D g);

// private method declarations
static void move_file_pointer(FILE *file_ptr, int lineNumber, int r);
static void get_location(FILE *vti, int *l, int b_l, Grid3D g);
static void get_distance(FILE *vti, double *d, double b_d, Grid3D g);

static int npes;          // Number of PEs
static int my_rank;       // Rank of the PE
static char * fileName;
static char * outfileName;

static int NX, NY, NZ, totalNodes;

int main(int argc, char *argv[]) {

  // MPI startup routine
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &npes);

  fileName    = argv[1];
  outfileName = argv[2];
  FILE *f = fopen(fileName, "r");

  double dims[6];
  vti_get_dimensions(f, dims);
  NX = dims[0] + 3;
  NY = dims[1] + 3;
  NZ = dims[2] + 3;

  totalNodes = NX * NY * NZ;
  Phi *p = (Phi *) malloc(sizeof(Phi));
  p->location = (int *) malloc(sizeof(int) * totalNodes);
  p->distance = (double *) malloc(sizeof(double) * totalNodes);

  p->dx = dims[3]; p->dy = dims[4]; p->dz = dims[5];

  vti_get_data( f, p->location, DEFAULT_BORDER_LOCATION,
                   p->distance, DEFAULT_BORDER_DISTANCE,
                make_grid3D(NX, NY, NZ));


  update_distance(p, totalNodes);

  calc_dist_field(p, totalNodes);

  MPI_Finalize();
  return 0;
}

void calc_dist_field( Phi * p, int totalNodes ) {

  int sweepNumber = my_rank + 1;
  double * tmp_dist;

  MPI_Barrier(MPI_COMM_WORLD);
  if(my_rank == 0){
    tmp_dist = (double *) malloc( totalNodes * sizeof(double) );
  }

  // sn represents the sweep number
  for( int sn = sweepNumber; sn <= NPES; sn += npes) {
    SweepInfo s = make_sweepInfo(p, sn);

    printf("PE: [%d] - performing sweep number ..... [%d/%d]\n", my_rank, sn, NPES);

    fast_sweep(p, &s);

    printf("PE: [%d] - completed sweep number ...... [%d/%d]\n", my_rank, sn, NPES);

  }

  MPI_Barrier(MPI_COMM_WORLD);
//  #pragma acc update host(p->distance[0:totalNodes])
  MPI_Reduce(p->distance, tmp_dist, totalNodes, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);

  if( my_rank == 0 ) {
    free( p->distance );
    p->distance = tmp_dist;
    set_distance_negative_inside(p, totalNodes);
    adjust_boundary(p);
    write_to_file(p->distance);
    printf("%s file created\n", outfileName);
  }

}

static void update_distance(Phi * p, int totalNodes) {

  int    *l = &p->location[0];
  double *d = &p->distance[0];

  for(int i = 0; i < totalNodes; i++) {
    if(*l != DEFAULT_BORDER_LOCATION && *d != DEFAULT_BORDER_DISTANCE ) {
      *d = (*l == 1 && *d == INFINITY) ? -1 : (*d > 0.0 || *d < 0.0) ? *d : DEFAULT_INTERIOR_DISTANCE;
    }
    l++; d++;
  }

}



void write_to_file(double * dist) {

  int x = NX;
  int y = NY;
  int z = NZ;
  char fname[255];
  sprintf(fname, "%s.txt", outfileName);
  FILE *fp = fopen(fname, "w");

  int i,j,k;
  double *t = &dist[0];
  for(i = 0; i < z; i++){
    for(j = 0; j < y; j++){
      for(k = 0; k < x; k++) {
    fprintf(fp, "%f  ", *(t++));
      }
      fprintf(fp, "\n");
    }
    fprintf(fp, "\n");
  }

}


static SweepInfo make_sweepInfo( Phi * p, int my_rank ) {
  SweepInfo s;

  s.order      = my_rank;
  s.firstLevel = 3;
  s.lastLevel  = (NX + NY + NZ) - 6;

  s.xDim = NX-2; s.dx = p->dx;
  s.yDim = NY-2; s.dy = p->dy;
  s.zDim = NZ-2; s.dz = p->dz;

  s.xSweepOff = (s.order == 4 || s.order == 8 ) ? s.xDim + 1 : 0;
  s.ySweepOff = (s.order == 2 || s.order == 6 ) ? s.yDim + 1 : 0;
  s.zSweepOff = (s.order == 3 || s.order == 7 ) ? s.zDim + 1 : 0;

  return s;
}

static void fast_sweep( Phi * p, SweepInfo * s ) {

  int start, end, incr;

  start = ( s->order == 2 || s->order == 5 || s->order == 7 || s->order == 8 ) ? s->lastLevel : s->firstLevel;

  if ( start == s->firstLevel ) {
    end  = s->lastLevel + 1;
    incr = 1;
  }
  else {
    end  = s->firstLevel - 1;
    incr = 0;
  }

  int max_x  = s->xDim + 2;
  int max_y  = s->yDim + 2;
  int max_xy = max_x * max_y;

  //#pragma acc data create(p[0:1]) copy(p->distance[0:totalNodes])
  for(int level = start; level != end; level = (incr) ? level+1 : level-1) {
    // s - start, e - end
    int xs, xe, ys, ye;

    xs = max(1, level-(s->yDim + s->zDim))    , ys = max(1,level-(s->xDim + s->zDim));
    xe = min(s->xDim, level-(s->firstLevel-1)), ye = min(s->yDim, level-(s->firstLevel-1));

    int x, y, z, i, j, k, index;
    int xSO = s->xSweepOff;
    int ySO = s->ySweepOff;
    int zSO = s->zSweepOff;
    #pragma acc kernels
    {
      #pragma acc loop independent
      for(x = xs; x <= xe; x++) {
        #pragma acc loop independent
        for(y = ys; y <= ye; y++) {
          z = level - (x+y);
          if(z > 0 && z <= NZ-2) {
            i = abs(z-zSO);
            j = abs(y-ySO);
            k = abs(x-xSO);
            index = i * max_xy + j * max_x + k;
            p->distance[index] = solveEikonal(p, index, NX, NY);
          }
        }
      } // end of acc parallel
    }
  }
}

#pragma acc routine seq
static double solveEikonal(Phi * p, int index, int max_x, int max_y) {

  int max_xy = max_x * max_y;

  double dist_new = 0;
  double dist_old = p->distance[index];

  double dx = p->dx, dy = p->dy, dz = p->dz;
  double minX = min(p->distance[index-1], p->distance[index+1]);
  double minY = min(p->distance[abs(index-max_x)], p->distance[abs(index+max_x)]);
  double minZ = min(p->distance[abs(index-max_xy)],p->distance[abs(index+max_xy)]);

  double m[] = { minX, minY, minZ} ;
  double d[] = { dx, dy, dz};

  // sort the mins 
  for(int i = 1; i < 3; i++){
    for(int j = 0; j < 3-i; j++) {
      if(m[j] > m[j+1]) {
        double tmp_m = m[j];
        double tmp_d = d[j];
        m[j] = m[j+1]; d[j] = d[j+1];
        m[j+1] = tmp_m; d[j+1] = tmp_d;
      }
    }
  }

  // simplifying the variables
  double m_0 = m[0], m_1 = m[1], m_2 = m[2];
  double d_0 = d[0], d_1 = d[1], d_2 = d[2]; 
  double m2_0 = m_0 * m_0, m2_1 = m_1 * m_1, m2_2 = m_2 * m_2;
  double d2_0 = d_0 * d_0, d2_1 = d_1 * d_1, d2_2 = d_2 * d_2;

  dist_new = m_0 + d_0;
  if(dist_new > m_1) {

    double s = sqrt(- m2_0 + 2 * m_0 * m_1 - m2_1 + d2_0 + d2_1); 
    dist_new = ( m_1 * d2_0 + m_0 * d2_1 + d_0 * d_1 * s) / (d2_0 + d2_1);

    if(dist_new > m_2) {

      double a = sqrt(- m2_0 * d2_1 - m2_0 * d2_2 + 2 * m_0 * m_1 * d2_2
                      - m2_1 * d2_0 - m2_1 * d2_2 + 2 * m_0 * m_2 * d2_1
                      - m2_2 * d2_0 - m2_2 * d2_1 + 2 * m_1 * m_2 * d2_0
                      + d2_0 * d2_1 + d2_0 * d2_2 + d2_1 * d2_2);

      dist_new = (m_2 * d2_0 * d2_1 + m_1 * d2_0 * d2_2 + m_0 * d2_1 * d2_2 + d_0 * d_1 * d_2 * a) /
                  (d2_0 * d2_1 + d2_0 * d2_2 + d2_1 * d2_2);
    }
  }

  return min(dist_old, dist_new);
}



static void set_distance_negative_inside(Phi * p, int totalNodes) {

  int    *l = &p->location[0];
  double *d = &p->distance[0];

  for(int i = 0; i < totalNodes; i++) {
    if(*l != DEFAULT_BORDER_LOCATION && *d != DEFAULT_BORDER_DISTANCE ) {
      if( *l == 1) *d = -1;
    }
    l++; d++;
  }

}

static void adjust_boundary( Phi * p ) {

  int x, y, z, xy, i, j, k;
  x  = NX;
  y  = NY;
  z  = NZ;
  xy = x * y;

  for(i = 0; i < z; i++){
    for(j = 0; j < y; j++){
      for(k = 0; k < x; k++){
        int I = i, J = j, K = k;
        I = (i == z-1) ? I-1 : (!i) ? I+1 : I;
        J = (j == y-1) ? J-1 : (!j) ? J+1 : J;
        K = (k == x-1) ? K-1 : (!k) ? K+1 : K;
        if( i != I || j != J || k != K) {
            int l_index = i * xy + j * x + k;
            int r_index = I * xy + J * x + K;
            p->distance[l_index] = p->distance[r_index];
        }
      }
    }
  }
}



/**************** vti_parser ********************************/

static void move_file_pointer(FILE *file_ptr, int lineNumber, int r) {
    char tmpStr[512];
    if(r) rewind(file_ptr);
    while (lineNumber > 0){
        fgets (tmpStr, 511, file_ptr);
        lineNumber--;
    }
}

void vti_get_dimensions(FILE *vti, double *d) {
    char tmpStr[512];
    rewind(vti);
    while (1) {
        fgets (tmpStr, 511, vti);
        if ( strstr(tmpStr, "ImageData WholeExtent") ) {
            sscanf(tmpStr, "    <ImageData WholeExtent=\"0 %lf 0 %lf 0 %lf\" Spacing=\"%lf %lf %lf\">",
                    &d[0], &d[1], &d[2], &d[3], &d[4], &d[5]);
            break;
        }
    }
}


void vti_get_data(FILE *vti, int *l, int b_l, double *d, double b_d, Grid3D g) {

    // move the file pointer to
    // line 6 from beginning
    move_file_pointer(vti, 6, 1);

    get_location(vti, l, b_l, g);

    // move the file pointer 2 lines
    // forward from its last position
    move_file_pointer(vti, 2, 0);

    get_distance(vti, d, b_d, g);

}

static void get_location(FILE *vti, int *l, int b_l, Grid3D g) {
    int i, j, k, *t = &l[0];
    for (i = 0; i < g.z; i++){
        for (j = 0; j < g.y; j++) {
            for (k = 0; k < g.x; k++) {
                // Border
                if (k == 0 || k == g.x-1 || j == 0 || j == g.y-1 || i == 0 || i == g.z-1 ) {
                    *(t++) = b_l;
                }
                else{ // Interior
                    fscanf(vti, "%d ", t++);
                }
            }
        }
    }
}



static void get_distance(FILE *vti, double *d, double b_d, Grid3D g) {
    int i, j, k;
    double *t = &d[0];
    for (i = 0; i < g.z; i++){
        for (j = 0; j < g.y; j++) {
            for (k = 0; k < g.x; k++) {
                // Border distance
                if (k == 0 || k == g.x-1 || j == 0 || j == g.y-1 || i == 0 || i == g.z-1 ) {
                    *(t++) = b_d;
                }
                else{ // Interior distance
                    fscanf(vti, "%lf ", t++);
                }
            }
        }
    }
}

Grid3D make_grid3D(int x, int y, int z){
    Grid3D g;
    g.x = x; g.y = y; g.z = z;

    return g;
}

这是我的编译命令:

pgc++ -fast -acc -ta=tesla:managed -Minfo=accel -I/opt/pgi/linux86-64/15.7/mpi/mpich/include -L/opt/pgi/linux86-64/15.7/mpi/mpich/lib -lmpi t2.c -o t2

这是输出:

$ LD_LIBRARY_PATH=/opt/pgi/linux86-64/15.7/mpi/mpich/lib ./t2 input.vti output
PE: [0] - performing sweep number ..... [1/8]
PE: [0] - completed sweep number ...... [1/8]
PE: [0] - performing sweep number ..... [2/8]
PE: [0] - completed sweep number ...... [2/8]
PE: [0] - performing sweep number ..... [3/8]
PE: [0] - completed sweep number ...... [3/8]
PE: [0] - performing sweep number ..... [4/8]
PE: [0] - completed sweep number ...... [4/8]
PE: [0] - performing sweep number ..... [5/8]
PE: [0] - completed sweep number ...... [5/8]
PE: [0] - performing sweep number ..... [6/8]
PE: [0] - completed sweep number ...... [6/8]
PE: [0] - performing sweep number ..... [7/8]
PE: [0] - completed sweep number ...... [7/8]
PE: [0] - performing sweep number ..... [8/8]
PE: [0] - completed sweep number ...... [8/8]
output file created
$

在我的例子中，它似乎比非 OpenACC 版本运行得更快(编译时没有 -acc -ta=tesla:managed -Minfo=accel)并且我 diff-ed 创建的输出文件，它们在非 OpenACC 和 OpenACC 版本之间是相同的。

我还尝试使用 2 个 MPI 等级运行此代码。它似乎运行没有崩溃:

$ LD_LIBRARY_PATH=/opt/pgi/linux86-64/15.7/mpi/mpich/lib /opt/pgi/linux86-64/15.7/mpi/mpich/bin/mpirun -n 2 ./t2 input.vti output
PE: [1] - performing sweep number ..... [2/8]
PE: [0] - performing sweep number ..... [1/8]
PE: [0] - completed sweep number ...... [1/8]
PE: [0] - performing sweep number ..... [3/8]
PE: [1] - completed sweep number ...... [2/8]
PE: [1] - performing sweep number ..... [4/8]
PE: [0] - completed sweep number ...... [3/8]
PE: [0] - performing sweep number ..... [5/8]
PE: [1] - completed sweep number ...... [4/8]
PE: [1] - performing sweep number ..... [6/8]
PE: [0] - completed sweep number ...... [5/8]
PE: [0] - performing sweep number ..... [7/8]
PE: [1] - completed sweep number ...... [6/8]
PE: [1] - performing sweep number ..... [8/8]
PE: [0] - completed sweep number ...... [7/8]
PE: [1] - completed sweep number ...... [8/8]
output file created
$

输出数据文件与单秩版本不同，但与双秩非 OpenACC 运行生成的版本相匹配。因此，如果还有任何遗留问题，我认为它们与 MPI 相关，而不与 OpenACC 相关。

编辑:

为了扩展最后一点，让我们将 OpenACC 排除在外，而只使用 PGI 15.7 工具链附带的 MPI (MPICH):

$ /opt/pgi/linux86-64/15.7/mpi/mpich/bin/mpicc t2.c -o t2 -lmpi
PGC-W-0129-Floating point overflow. Check constants and constant expressions (t2.c: 88)
PGC-W-0129-Floating point overflow. Check constants and constant expressions (t2.c: 142)
PGC-W-0129-Floating point overflow. Check constants and constant expressions (t2.c: 143)
PGC-W-0129-Floating point overflow. Check constants and constant expressions (t2.c: 308)
PGC/x86-64 Linux 15.7-0: compilation completed with warnings
$ LD_LIBRARY_PATH=/opt/pgi/linux86-64/15.7/mpi/mpich/lib /opt/pgi/linux86-64/15.7/mpi/mpich/bin/mpirun -n 1 ./t2 input.vti output1rank
PE: [0] - performing sweep number ..... [1/8]
PE: [0] - completed sweep number ...... [1/8]
PE: [0] - performing sweep number ..... [2/8]
PE: [0] - completed sweep number ...... [2/8]
PE: [0] - performing sweep number ..... [3/8]
PE: [0] - completed sweep number ...... [3/8]
PE: [0] - performing sweep number ..... [4/8]
PE: [0] - completed sweep number ...... [4/8]
PE: [0] - performing sweep number ..... [5/8]
PE: [0] - completed sweep number ...... [5/8]
PE: [0] - performing sweep number ..... [6/8]
PE: [0] - completed sweep number ...... [6/8]
PE: [0] - performing sweep number ..... [7/8]
PE: [0] - completed sweep number ...... [7/8]
PE: [0] - performing sweep number ..... [8/8]
PE: [0] - completed sweep number ...... [8/8]
output1rank file created
$ LD_LIBRARY_PATH=/opt/pgi/linux86-64/15.7/mpi/mpich/lib /opt/pgi/linux86-64/15.7/mpi/mpich/bin/mpirun -n 2 ./t2 input.vti output2rank
PE: [0] - performing sweep number ..... [1/8]
PE: [1] - performing sweep number ..... [2/8]
PE: [1] - completed sweep number ...... [2/8]
PE: [1] - performing sweep number ..... [4/8]
PE: [0] - completed sweep number ...... [1/8]
PE: [0] - performing sweep number ..... [3/8]
PE: [1] - completed sweep number ...... [4/8]
PE: [1] - performing sweep number ..... [6/8]
PE: [0] - completed sweep number ...... [3/8]
PE: [0] - performing sweep number ..... [5/8]
PE: [1] - completed sweep number ...... [6/8]
PE: [1] - performing sweep number ..... [8/8]
PE: [1] - completed sweep number ...... [8/8]
PE: [0] - completed sweep number ...... [5/8]
PE: [0] - performing sweep number ..... [7/8]
PE: [0] - completed sweep number ...... [7/8]
output2rank file created
$ diff -q output1rank.txt output2rank.txt
Files output1rank.txt and output2rank.txt differ
$

关于c - 编写混合 MPI/OpenACC 程序，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/33130676/

文章推荐： html - 导航栏在跳跃

文章推荐： c - setlocale() 没有按预期工作

文章推荐： java - 如何在 netbeans JFrame 编辑器中显示自定义对话框？

文章推荐： c - 使用 CURL 测试时，我的 C 程序在接受连接时被卡住

c - 如何从客户端(用 C 编写)接收 int 数组到服务器(用 python 编写)
我只想从客户端向服务器发送数组 adc_array=[w, x, y, z]。下面是客户端代码，而我的服务器是在只接受 json 的 python 中。编译代码时我没有收到任何错误，但收到 2 条警告
node.js - 如何连接我的移动应用程序(用 lua 编写)和我的服务器(用 node.js 编写)？
我是 lua 和 Node js 的新手，我正在尝试将我正在开发的移动应用程序连接到服务器。问题是它连接到服务器，但我尝试传递的数据丢失或无法到达服务器。对我正在做的事情有什么问题有什么想法吗？ th
Haskell 编写 myLength
我在这个页面上工作 http://www.haskell.org/haskellwiki/99_questions/Solutions/4 我理解每个函数的含义，看到一个函数可以像这样以多种方式定义，
Java CSV 编写
我目前正在尝试将数据写入 excel 以生成报告。我可以将数据写入 csv 文件，但它不会按照我想要的顺序出现在 excel 中。我需要数据在每列的最佳和最差适应性下打印，而不是全部打印在平均值下。这
Java - 编写、读取和修改带参数的字符串
所以，我正在做一个项目，现在我有一个问题，所以我想得到你的帮助:) 首先，我已经知道如何编写和读取 .txt 文件，但我想要的不仅仅是 x.hasNext()。我想知道如何像 .ini 那样编写、读
javascript - 编写 For 循环来计算阶乘
我正在尝试编写一个函数，该函数将返回作为输入给出的任何数字的阶乘。现在，我的代码绝对是一团糟。请帮忙。 function factorialize(num) { for (var i=num, i
Javascript，编写 if 条件的更好方法
这个问题已经有答案了: Check variable equality against a list of values (16 个回答) 已关闭 4 年前。有没有一种简洁或更好的方法来编写这个条件
aframe - 编写 A 型框架的测试规范
我对 VR 完全陌生，正在 AFrame 中为一个类(class)项目开发 VR 太空射击游戏，并且想知道 AFrame 中是否有 TDD 的任何文档/标准。有人能指出我正确的方向吗？最佳答案几乎
javascript - 编写 for 循环以使用数组创建多个方法
我正在尝试创建一个 for 循环，它将重现以下功能代码块，但以一种更具吸引力的方式。这是与 Soundcould 小部件 API 实现一起使用的 here on stackoverflow $(doc
Java 编写/编辑属性文件
我有一个非常令人困惑的问题。我正在尝试更改属性文件中的属性，但它只是没有更改... 这是代码: package config; import java.io.FileNotFoundException
aframe - 编写 A 型框架的测试规范
我对 VR 完全陌生，正在 AFrame 中为一个类(class)项目开发 VR 太空射击游戏，并且想知道 AFrame 中是否有 TDD 的任何文档/标准。有人能指出我正确的方向吗？最佳答案几乎
.net - 编写.NET互操作调试器
我正在开发一个用户模式(Ring3)代码级调试器。它还应支持.NET可执行文件的本机(x86)调试。基本上，我需要执行以下操作: 1).NET在隐身模式下加载某些模块，而没有LOAD_DLL_DEBU
python - 编写 if 语句以避免某些列表项的更好方法是什么？
我有一个列表，我知道有些项目是不必要打印的，我正在尝试通过 if 语句来做到这一点...但是它变得非常复杂，所以有没有什么方法可以在 if 语句中包含多个索引而无需打印重写整个声明。看起来像这样的东
c# - 编写 if 语句是否会以不同方式影响程序的速度和效率？
我很好奇以不同方式编写 if 语句是否会影响程序的速度和效率。所以，例如写一个这样的: bool isActive = true; bool isResponding = false; if (isA
javascript - 编写 if 语句的新方法
我在搜索网站的源代码时找到了一种以另一种方式(我认为)编写 if 语句的方法。代替: if(a)b; 或: a?b:''; 我读了: !a||b; 第三种方式和前两种方式一样吗？如果是，为什么我们要
Java + 编写 XML
我的数据采用以下格式(HashMap的列表) {TeamName=India, Name=Sachin, Score=170} {TeamName=India, Name=Sehwag, Score=
mysql - 编写 HAVING 条件的最有效方法
我目前正在完成 More JOIN operations sqlzoo 的教程，遇到了下面的代码作为#12 的答案: SELECT yr,COUNT(title) FROM movie JOIN ca
ruby - 编写 && 检查列表的更好方法？
我正试图找到一种更好的方法来编写这段代码: def down_up(array, player) 7.downto(3).each do |row| 8.times do |col
由 C++ 编写
出于某种原因，我的缓冲区中充满了乱码，我不确定为什么。我什至用十六进制编辑器检查了我的文件，以验证我的字符是否以 2 字节的 unicode 格式保存。我不确定出了什么问题。 [打开文件] fseek
c# - 编写 FizzBuzz
阅读编码恐怖片时，我刚刚又遇到了 FizzBuzz。原帖在这里:Coding Horror: Why Can't Programmers.. Program? 对于那些不知道的人:FizzBu

太空宇宙

个人简介

我是一名优秀的程序员,十分优秀！

作者热门文章

滴滴打车优惠券免费领取

全站热门文章

首页

博学

6Ren·AI

商城

c - 编写混合 MPI/OpenACC 程序