我想写一个并行代码,在3D矩阵上工作,每个进程都有自己的子矩阵,但为了完成任务,他们需要一些关于相邻进程的子矩阵(只是边界平面)的信息。我通过点对点通信发送这些信息,但我知道对于大矩阵来说,这不是一个好主意,所以我决定使用派生数据类型进行通信。我对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();
}
对于3d矩阵,通常您必须使用一个矢量的矢量(因为涉及两个步长)-这是可能的,但更简单的是使用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中一直都会出现。
感谢Jonathan Dursi。在这里,我想发布完整的代码,创建一个三维矩阵,并使用派生的数据类型进行通信(只有带有常量y的平面才会从一个进程发送到另一个进程)。我使用了Jonathan Dursi在上面发布的函数。
#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();
}
很抱歉告诉您,您审阅的代码仍然无法正常工作。输出似乎是正确的原因是IE和KE相等。如果你把它们不同,你会看到,这些值是用交替的Y索引写的。
如果你看看Jonathan Durst的代码样本的内存分配,它是这样的:
[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);