Как использовать производный тип данных MPI для 3D-массива?

Question

Я хочу написать параллельный код, который работает на трехмерной матрице, где каждый процесс имеет свою собственную подматрицу, но для выполнения своей работы им нужна некоторая информация о подматрице соседних процессов (только граничные плоскости).Я посылаю эту информацию с двухточечной связью, но я знаю, что для большой матрицы это не очень хорошая идея, поэтому я решил использовать производный тип данных для связи.У меня проблема с mpi_type_vector: например, у меня есть матрица NX*NY*NZ, и я хочу отправить плоскость с константой NY другому процессу. Я пишу эти строки для этого:

MPI_Datatype sub;

MPI_Type_vector(NX, NZ, NY*NZ, MPI_DOUBLE, &sub);

MPI_Type_commit(&sub);

, но этоне работает (не могу отправить желаемый самолет).Что случилось?мой тестовый код здесь:

#include <mpi.h>
#include <iostream>

using namespace std;

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

    int const IE=100,JE=25,KE=100;
    int size,rank;
    MPI_Status status;

    MPI_Init(&argc,&argv);
    MPI_Comm_size(MPI_COMM_WORLD,&size);
    MPI_Comm_rank(MPI_COMM_WORLD,&rank);
    MPI_Datatype sub;
    MPI_Type_vector(KE,IE,IE+(JE-1)*IE,MPI_DOUBLE,&sub);
    MPI_Type_commit(&sub);

    if (rank==0){

        double*** a=new double**[IE];

        for(int i=0;i<IE;i++){
            a[i]=new double *[JE];
            for(int j=0;j<JE;j++){
                a[i][j]=new double [KE];
            }
        }

        for(int i=0;i<IE;i++){
            for(int j=0;j<JE;j++){
                for(int k=0;k<KE;k++){
                    a[i][j][k]=2;
                }}}

        for(int i=0;i<IE;i++){
            for(int j=0;j<JE;j++){
                a[i][j][0]=2;
            }}

        MPI_Send(&a[0][0][0],1,sub,1,52,MPI_COMM_WORLD);

    }

    if (rank==1){

        double*** b=new double**[IE];

        for(int i=0;i<IE;i++){
            b[i]=new double *[JE];
            for(int j=0;j<JE;j++){
                b[i][j]=new double [KE];
            }
        }

        for(int i=0;i<IE;i++){
            for(int j=0;j<JE;j++){
                for(int k=0;k<KE;k++){
                    b[i][j][k]=0;
                }}}

        MPI_Recv(&b[0][0][0][0],1,sub,0,52,MPI_COMM_WORLD,&status);

        for(int i=0;i<IE;i++){
            for(int j=0;j<JE;j++){
                for(int k=0;k<KE;k++){
                    if(b[i][j][k]>0){
                        cout<<"b["<<i<<"]["<<j<<"]["<<k<<"]="<<b[i][j][k]<<endl;
                    }}}}

    }

    MPI_Finalize();

}

Answer 1

В случае с трехмерной матрицей, как правило, вам придется использовать вектор векторов (поскольку в нем задействованы два шага) - что возможно, но гораздо проще использовать MPI_Type_create_subarray () , который простопозволяет вам вырезать кусок многомерного массива, который вам нужен.

Обновление : Одна из проблем в приведенном выше коде состоит в том, что выделенный вами 3d-массив не является смежным;это коллекция IE * JE, выделенных 1d массивов, которые могут быть или не быть рядом друг с другом.Поэтому нет надежного способа извлечь из него плоскость данных.

Вам нужно сделать что-то вроде этого:

double ***alloc3d(int l, int m, int n) {
    double *data = new double [l*m*n];
    double ***array = new double **[l];
    for (int i=0; i<l; i++) {
        array[i] = new double *[m];
        for (int j=0; j<m; j++) {
            array[i][j] = &(data[(i*m+j)*n]);
        }
    }
    return array;
}

Тогда данные будут в одном большом кубе, как выЯ ожидаю, с массивом указателей, указывающих на это.Это - тот факт, что C не имеет реальных многомерных массивов - все время возникает с C + MPI.

Answer 2

Сожалею, что ваш проверенный код все еще работает неправильно. Причина, по которой вывод кажется правильным, заключается в том, что IE и KE равны. Если вы их отличите, вы увидите, что значения записаны в чередующихся Y-индексах.

Если вы посмотрите на распределение памяти в примере кода Джонатана Дерста, которое выглядит следующим образом:

[x0y0z0] [x0y0z1] [x0y1z0] [x0y1z1] [x1y0z0] [x1y0z1] [x1y1z0] [x1y1z1]   //or
{x0:(y0:[z0,z1]) ; (y1:[z0,z1])} ; {x1:(y0:[z0,z1]) ; (y1:[z0,z1])} //nx=ny=nz=2 
        <bl.len>  
 X             count                X  
     |<-          stride               ->| 

вы увидите, что у вас есть количество блоков nx с длиной блока nz значений и шаг между ними ny * nz.

Ваш код работает правильно, если вы измените тип данных на:

MPI_Type_vector(IE,KE,KE*JE,MPI_DOUBLE,&sub);

Answer 3

Спасибо Джонатану Дурси.Здесь я хочу опубликовать полный код, который создает трехмерную матрицу и использует производные типы данных для связи (только плоскость с константой y будет передаваться из одного процесса в другой). Я использовал функцию Джонатана Дурси, опубликованную выше.

#include <mpi.h>
#include <iostream>
#include <math.h>
#include <fstream>
#include <vector>
using namespace std;
#define IE 100
#define JE 50
#define KE 100

#define JE_loc 52
double ***alloc3d(int l, int m, int n) {
double *data = new double [l*m*n];
double ***array = new double **[l];
for (int i=0; i<l; i++) {
    array[i] = new double *[m];
    for (int j=0; j<m; j++) {
        array[i][j] = &(data[(i*m+j)*n]);
    }
}
return array;
}




int main(int argc ,char ** argv)
{
//////////////////////declartion/////////////////////////////
int const NFREQS=100,ia=7,ja=7,ka=7;
double const pi=3.14159;
int i,j,size,rank,k;
//MPI_Status status[10];
MPI_Status status;
MPI_Request request[10];
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Datatype sub;
MPI_Type_vector(KE,IE,IE+(JE-1)*IE,MPI_DOUBLE,&sub);
MPI_Type_commit(&sub);


double ***a=alloc3d(IE,JE,KE);
for (i=0; i<IE; i++) {
    for (j=0; j<JE; j++) {
        for (k=0; k<KE; k++) {
            a[i][j][k]=0.0;
        }
    }
}
if (rank==0) {
    for (i=0; i<IE; i++) {
        for (j=0; j<JE; j++) {
            for (k=0; k<KE; k++) {
                a[i][j][k]=2;
            }
        }
    }
    MPI_Send(&a[0][0][0],1,sub,1,52,MPI_COMM_WORLD);

}
if (rank==1) {
    MPI_Recv(&a[0][49][0],1,sub,0,52,MPI_COMM_WORLD,&status);
    for (i=0; i<IE; i++) {
        for (j=0; j<JE; j++) {
            for (k=0; k<KE; k++) {
                if (a[i][j][k]>0) {
                    cout<<"a["<<i<<"]["<<j<<"]["<<k<<"]="<<a[i][j][k]<<endl;
                }
            }
        }
    }

}
MPI_Finalize();
}