基于都志辉老师《MPI并行程序设计模式》第14章内容。
前面接触到的MPI发送的数据类型都是连续型的数据。非连续类型的数据，MPI也可以发送，但是需要预先处理，大概有两类方法：
（1）用户自定义新的数据类型，又称派生类型（类似定义结构体类型，但是比结构体复杂，需要考虑<类型，偏移量>两方面的内容）
（2）数据的打包和解包（将不连续的数据给压缩打包到连续的区域，然后再发送；接受到打包数据后，先解包再使用）
这样做的好处，我猜一个是可以有效减少通信的次数，提高程序效率；另一方面可以减轻程序员设计程序的负担，降低维护成本。
下面记录学习中的一些注意的点：
1. 通用的数据类型描述方法——类型图。MPI自定义派生类型的变量，往往由一些不同类型（MPI_INT, MPI_FLOAT, MPI_DOUBLE等）的变量按一定组织形式（每个变量偏移量是多少）组合。因此，引入类型图来形象描述MPI的派生变量是由那些不同的基础类型，按照怎么样的偏移量组合在一起的。
类型图 = {<基类型0，偏移0>，<基类型0，偏移0>，...，<基类型0，偏移0>}
先看书理解好什么是类型图，为什么要有类型图，学习MPI派生类型就顺畅一些。
另外，与结构体类似，MPI派生的数据类型，需要考虑内存对齐的因素：不同基变量的偏移量设计是有内存对齐讲究的。具体有什么讲究，为什么要有讲究，可以复习之前的学习blog（http://www.cnblogs.com/xbf9xbf/p/5121748.html）
2. 几种常见的新类型数据的定义方法
可以直接看书上给的实际例子，看例子就懂了什么意思了。在看例子的时候，注意一下偏移量与内存对齐的关系就可以了。
另外，在定义好新类型之后，还需要在MPI中注册新定义好的数据类型。
直接看下面这个综合的例子。
代码实现的内容是测试MPI用不同的方式发送非连续数据的传输效率。其中用到了几种新类型数据的定义方法。
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h> #define NUMBER_OF_TEST 10 int main(int argc, char *argv[])
{
MPI_Datatype vec1, vec_n;
int blocklens[];
MPI_Aint indices[]; // A=array
MPI_Datatype old_types[];
double *buf, *lbuf;
register double *in_p, *out_p;
int rank;
int n, stride;
double t1, t2, tmin;
int i,j,k,nloop;
MPI_Status status; MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
n = ;
stride = ;
nloop = /n;
buf = (double*)malloc(n*stride*sizeof(double));
if (!buf) MPI_Abort(MPI_COMM_WORLD,);
lbuf = (double*)malloc(n*sizeof(double));
if (!lbuf) MPI_Abort(MPI_COMM_WORLD,);
if (==rank) printf("Kind\tn\tstride\ttime(sec)\tRate(MB/sec)\n");
// 1. Vector数据传输测试: Vector是大结构 里面包含n个vec1结构
// 构造一个double vector
MPI_Type_vector(n, , stride, MPI_DOUBLE, &vec1);
MPI_Type_commit(&vec1);
if (==rank) {
MPI_Aint ext[];
MPI_Type_extent(vec1, ext);
printf("extent of vec1 : %d\n",(int)(*ext));
MPI_Type_extent(MPI_DOUBLE, ext);
printf("exten of MPI_DOUBLE : %d\n",(int)(*ext));
}
tmin = ;
for (k=; k<NUMBER_OF_TEST; k++) {
if (==rank) {
// 保证通信双方都ready
MPI_Sendrecv(MPI_BOTTOM, , MPI_INT, , , MPI_BOTTOM, , MPI_INT, , , MPI_COMM_WORLD, &status);
t1 = MPI_Wtime();
for (j=; j<nloop; j++) {
MPI_Send(buf, , vec1, , k, MPI_COMM_WORLD);
MPI_Recv(buf, , vec1, , k, MPI_COMM_WORLD, &status);
}
t2 = (MPI_Wtime()-t1) / nloop;
tmin = tmin>t2 ? t2 : tmin;
}
else if (==rank) {
// 保证通信双方都ready
MPI_Sendrecv(MPI_BOTTOM, , MPI_INT, , , MPI_BOTTOM, , MPI_INT, , , MPI_COMM_WORLD, &status);
for (j=; j<nloop; j++) {
MPI_Recv(buf, , vec1, , k, MPI_COMM_WORLD, &status);
MPI_Send(buf, , vec1, , k, MPI_COMM_WORLD);
}
}
}
tmin = tmin / ;
if (==rank) printf("Vector\t%d\t%d\t%f\t%f\n",n,stride,tmin,n*sizeof(double)*1.0e-6/tmin);
MPI_Type_free(&vec1);
// 2. 可变向量类型传输测试: Struct是小结构 每个struct由vec_n构成
blocklens[] = ;
blocklens[] = ;
indices[] = ;
indices[] = stride*sizeof(double);
old_types[] = MPI_DOUBLE;
old_types[] = MPI_UB; // 上限区间占位符 不占大小 只占位置
MPI_Type_struct(, blocklens, indices, old_types, &vec_n);
MPI_Type_commit(&vec_n);
if (==rank) {
MPI_Aint ext[];
MPI_Type_extent(vec_n, ext);
printf("extent of vec_n : %d\n",(int)(*ext));
}
tmin = ;
for (k=; k<NUMBER_OF_TEST; k++) {
if (==rank) {
MPI_Sendrecv(MPI_BOTTOM, , MPI_INT, , , MPI_BOTTOM, , MPI_INT, , , MPI_COMM_WORLD, &status);
t1 = MPI_Wtime();
for (j=; j<nloop; j++) {
MPI_Send(buf, n, vec_n, , k, MPI_COMM_WORLD);
MPI_Recv(buf, n, vec_n, , k, MPI_COMM_WORLD, &status);
}
t2 = (MPI_Wtime()-t1) / nloop;
tmin = tmin>t2 ? t2 : tmin;
}
else if (==rank) {
MPI_Sendrecv(MPI_BOTTOM, , MPI_INT, , , MPI_BOTTOM, , MPI_INT, , , MPI_COMM_WORLD, &status);
for (j=; j<nloop; j++) {
MPI_Recv(buf, n, vec_n, , k, MPI_COMM_WORLD, &status);
MPI_Send(buf, n, vec_n, , k, MPI_COMM_WORLD);
}
}
}
tmin = tmin / ;
if (==rank) printf("Struct\t%d\t%d\t%f\t%f\n",n,stride,tmin,n*sizeof(double)*1.0e-6/tmin);
MPI_Type_free(&vec_n);
// 3.User
tmin = ;
for (k=; k<NUMBER_OF_TEST; k++) {
if (==rank) {
MPI_Sendrecv(MPI_BOTTOM, , MPI_INT, , , MPI_BOTTOM, , MPI_INT, , , MPI_COMM_WORLD, &status);
t1 = MPI_Wtime();
for(j=; j<nloop; j++)
{
for (i=; i<n; i++) {
lbuf[i] = buf[i*stride];
}
MPI_Send(lbuf, n, MPI_DOUBLE, , k, MPI_COMM_WORLD);
MPI_Recv(lbuf, n, MPI_DOUBLE, , k, MPI_COMM_WORLD, &status);
for (i=; i<n; i++) {
buf[i*stride] = lbuf[i];
}
}
t2 = (MPI_Wtime()-t1) / nloop;
tmin = tmin>t2 ? t2 : tmin;
}
else if (==rank) {
MPI_Sendrecv(MPI_BOTTOM, , MPI_INT, , , MPI_BOTTOM, , MPI_INT, , , MPI_COMM_WORLD, &status);
for (j=; j<nloop; j++) {
MPI_Recv(lbuf, n, MPI_DOUBLE, , k, MPI_COMM_WORLD, &status);
for (i=; i<n; i++) {
buf[i*stride] = lbuf[i];
}
for (i=; i<n; i++) {
lbuf[i] = buf[i*stride];
}
MPI_Send(lbuf, n, MPI_DOUBLE, , k, MPI_COMM_WORLD);
}
}
}
tmin = tmin / 2.0;
if (==rank) printf("User(1)\t%d\t%d\t%f\t%f\n",n,stride,tmin,n*sizeof(double)*1.0e-6/tmin);
// 4. user-packing
tmin = ;
for (k=; k<NUMBER_OF_TEST; k++) {
if (==rank) {
MPI_Sendrecv(MPI_BOTTOM, , MPI_INT, , , MPI_BOTTOM, , MPI_INT, , , MPI_COMM_WORLD, &status);
t1 = MPI_Wtime();
for(j=; j<nloop; j++){
in_p = buf;
out_p = lbuf;
for(i=; i<n; i++){
out_p[i] = *in_p;
in_p += stride;
}
MPI_Send(lbuf, n, MPI_DOUBLE, , k, MPI_COMM_WORLD);
MPI_Recv(lbuf, n, MPI_DOUBLE, , k, MPI_COMM_WORLD, &status);
out_p = buf;
in_p = lbuf;
for (i=; i<n; i++) {
*out_p = in_p[i];
out_p += stride;
}
}
t2 = (MPI_Wtime()-t1) / nloop;
tmin = tmin>t2 ? t2 : tmin;
}
else if (==rank) {
MPI_Sendrecv(MPI_BOTTOM, , MPI_INT, , , MPI_BOTTOM, , MPI_INT, , , MPI_COMM_WORLD, &status);
for(j=; j<nloop; j++)
{
MPI_Recv(lbuf, n, MPI_DOUBLE, , k, MPI_COMM_WORLD, &status);
in_p = lbuf;
out_p = buf;
for (i=; i<n; i++) {
*out_p = in_p[i];
out_p += stride;
}
out_p = lbuf;
in_p = buf;
for (i=; i<n; i++) {
out_p[i] = *in_p;
in_p += stride;
}
MPI_Send(lbuf, n, MPI_DOUBLE, , k, MPI_COMM_WORLD);
}
}
}
tmin = tmin / ;
if (==rank) printf("User(2)\t%d\t%d\t%f\t%f\n",n,stride,tmin,n*sizeof(double)*1.0e-6/tmin);
MPI_Finalize();
}
代码执行结果如下：
结果分析；
（1）传输效率。可以看到前面vector和struct两种MPI提供的派生数据的方式（利用MPI_Type_vector和MPI_Type_struct），是要由于后面用户自己派生数据的方式的（程序员自己自己将1000个double数据放到一个buf中再发送）。
（2）extent的算法。每个MPI的派生类型变量，都有一个extent的概念。这里的extent就是考虑了内存对齐后的派生类型的新变量跨度大小。这里以MPI_Type_vector(1000, 1, 24, MPI_DOUBLE, &vec1)为例，为什么vec1的extent是191816。
具体算法如下：vec1的基础类型是MPI_DOUBLE，每个MPI_DOUBLE占8个字节；每个变量段的便宜量是24个MPI_DOUBLE（即，第1个double+23个double占位+第2个double+23个double占位,....）；这样重复999次，直到1000次的时候，最后一个double设定之后，后面就没有下一个doulbe了，因此也没有23个doulbe空间占位了（注意在算extent的时候不要把最后的23个double空间占位给算进去）。extent = 1000*24*8 - 23*8 = 191816
3.构造新数据类型时，计算地址偏移量的函数。
前面提到过，构造MPI派生类型时候，需要程序员指定每个变量相对于派生类型的入口地址（MPI_BOTTOM）偏移量是多少。前面的例子是通过人工设定的需要多少偏移量，MPI提供了一个函数MPI_Address(void *location, MPI_Aint *address)来帮我们解决这个问题。另外，还有个MPI_Type_size函数用来统计MPI派生变量中有用内容的大小，MPI_Type_extent函数测试MPI派生变量的跨度大小，下面的例子中也对这两个函数加以区分。
看如下的例子：
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h> int main(int argc, char *argv[])
{
int rank;
struct{
int a;
double b;
} value;
MPI_Datatype mystruct;
int blocklens[];
MPI_Aint indices[];
MPI_Datatype old_types[]; MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
blocklens[] = ;
blocklens[] = ;
old_types[] = MPI_INT;
old_types[] = MPI_DOUBLE;
MPI_Address(&value.a, &indices[]);
MPI_Address(&value.b, &indices[]);
indices[] = indices[]-indices[];
indices[] = ;
MPI_Type_struct(, blocklens, indices, old_types, &mystruct);
MPI_Type_commit(&mystruct);
MPI_Aint extent[];
MPI_Type_extent(mystruct, extent);
int size;
MPI_Type_size(mystruct, &size);
if (==rank) {
printf("stride:%d\n",(int)indices[]);
printf("extent:%d\n",(int)*extent);
printf("size of struct:%d\n",(int)sizeof(value));
printf("size of mystruct:%d\n", size);
} while (value.a>=) {
if (==rank) {
scanf("%d %lf",&value.a,&value.b);
}
MPI_Bcast(&value, , mystruct, , MPI_COMM_WORLD);
printf("Process %d got %d and %lf\n",rank,value.a,value.b);
}
MPI_Type_free(&mystruct);
MPI_Finalize();
}
代码执行结果如下：
分析如下：
（1）line23~25，就是用MPI_Address来计算派生类型中各个基础类型变量起始位置的绝对偏移量（间接通过一个{int, double}结构体类型作为标靶，来找到MPI目标派生类型各个基变量的偏移量的）
（2）定义的结构体value中虽然只有一个int和double，但是考虑内存对齐（double大小8的整数倍），所以size是16。
（3）同理，MPI派生类型mystruct经过内存对齐要求后，extent也是16。
（4）MPI_Type_size只要求有效数据大小，因此是int+double加一起是12个bytes。
4. 打包与解包
打包与解包并不是构造新的MPI数据类型，而是将不同类型的数据利用打包函数压缩到连续的发送缓冲区中；再通过解包函数，按照解包的规则，从接受缓冲区中将数据解包，进而后续使用。
打包函数和解包函数：
MPI_Pack(void *inbuf, int incount, MPI_datatype, void *outbuf, int outcount, int *position, MPI_Comm comm)
MPI_Unpack(void *inbuf, int insize, int *position, void *outbuf, int outcount, MPI_Datatype datatype, MPI_Comm comm)
打包后发送时的数据类型：
经过打包后的数据，也是存放在一个缓冲区中；此时再用MPI_Send函数发送时，需要指定发送类型为MPI_PACKED。
打包和解包在MPI通信中的位置：
a. 给定发送缓冲区、给定接受缓冲区
b. 将数据打包到发送缓冲区中（一个数据一个数据打包）
c. 将打包后的数据发送到目标进程
d. 将数据从目标进程的接收缓冲区中解包（一个数据一个数据解包，解包的顺序与打包的顺序相同）
可以看到，打包和解包在通信过程外围一层，是服务于MPI通信的（这里的通信可以是Send Recv也可以是Bcast这种广播模式）
下面看一个代码例子，root进程将一个整数和双精度数打包，然后广播给所有的进程，各进程分别将数据解包后再打印。
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h> int main(int argc, char *argv[])
{
int rank;
int packsize, position;
int a;
double b;
char packbuf[]; MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
while (a>=) {
if (==rank) { // 在进程0中准备发送的数据
scanf("%d %lf",&a,&b);
packsize = ;
MPI_Pack(&a, , MPI_INT, packbuf, , &packsize, MPI_COMM_WORLD);
MPI_Pack(&b, , MPI_DOUBLE, packbuf, , &packsize, MPI_COMM_WORLD);
}
MPI_Bcast(&packsize, , MPI_INT, , MPI_COMM_WORLD);
MPI_Bcast(packbuf, packsize, MPI_PACKED, , MPI_COMM_WORLD);
if (!=rank) { // 在进程1中处理进程0发送来的打包数据
position = ;
MPI_Unpack(packbuf, packsize, &position, &a, , MPI_INT, MPI_COMM_WORLD);
MPI_Unpack(packbuf, packsize, &position, &b, , MPI_DOUBLE, MPI_COMM_WORLD);
}
printf("Process %d got %d and %lf\n", rank, a, b);
}
MPI_Finalize();
return ;
}
执行结果如下：
总结一下，如果MPI不提供（1）自定义派生数据类型（2）打包和解包 两种非连续的数据发送方法，则发送非连续的组合数据，需要程序员写非常多的Send Recv通信对来完成，既降低了通信的效率，又增加了程序维护的成本。以后遇到类似的问题时候，可以回顾MPI这块的知识。