并行编程——MPI

MPI是Message Passing Interface的简称，通过这个协议可以在各个进程——尤其是分布式内存进程——间能够进行通信，交流消息共同完成一个任务。

进行mpi编程的基本流程如下

首先要载入头文件

Fortran 77:
    include 'mpif.h' 
Fortran 90:
    use mpi 
C/C++:
    #include "mpi.h" 

第二步是初始化MPI环境

Fortran 77/90:
    Call MPI_INIT(ierror)
其中integer ierror
C:
    int MPI_Init(int *argc, char ***argv); 
C++:
    int MPI::Init(int *argc, char ***argv); 

第三步是获知参与并行的核的总数

Fortran 77/90:
    Call MPI_COMM_SIZE(comm, size, ierror)
其中integer comm, size, ierror
C:
    int MPI_Comm_size (MPI_Comm comm, int *size); 
C++:
    int MPI::COMM_WORLD.Get_size( ); 

第四步是得知自己所在的进程的序列号

Fortran 77/90:
    Call MPI_COMM_RANK(comm, rank, ierror)
其中 integer comm, rank, ierror
C:
    int MPI_Comm_rank(MPI_Comm comm, int *rank); 
C++:
    int MPI::COMM_WORLD.Get_rank( ); 

进行相关计算后

第五步是结束MPI环境

Fortran 77/90:
     Call MPI_FINALIZE(ierror)
其中integer ierror
C:
     int MPI_Finalize(); 
C++:
     MPI::Finalize(); 

还有两个不常用的方法

计算逝去的时间

Fortran 77/90:
    double precision MPI_WTIME() 
C:
    double MPI_Wtime(); 
C++:
    double MPI::Wtime(); 

和中止MPI环境

Fortran:
    Call MPI_ABORT(comm, errorcode, ierror)
其中integer comm, errorcode, ierror
C:
    int MPI_Abort(MPI_Comm comm int errorcode ); 
C++:
    MPI::COMM_WORLD.Abort( int errorcode ); 

MPI的数据类型则包括以下几种

C data type	MPI data type
char	MPI_CHAR
short int	MPI_SHORT
int	MPI_INT
long int	MPI_LONG
float	MPI_FLOAT
double	MPI_DOUBLE
long double	MPI_LONG_DOUBLE

Fortran data type	MPI data type
INTEGER	MPI_INTEGER
REAL	MPI_REAL
REAL*8	MPI_REAL8
DOUBLE PRECISION	MPI_DOUBLE_PRECISION
COMPLEX	MPI_COMPLEX
LOGICAL	MPI_LOGICAL
CHARACTER	MPI_CHARACTER

通常所用的消息交换方法有以下几种

MPI_SEND&MPI_RECV

Fortran 77/90:
    Call MPI_SEND(BUF, COUNT, DATATYPE, DEST, TAG, COMM, IERROR)
其中<type> BUF(*)；INTEGER COUNT, DATATYPE, DEST, TAG, COMM, IERROR
C:
    int MPI_Send(void* buf, int count, MPI_Datatype datatype, int dest,int tag, MPI_Comm comm)
C++:
    MPI::Comm::Send(const void* buf, int count, const MPI::Datatype& datatype, int dest, int tag)

Fortran 77/90:
    Call MPI_RECV(BUF, COUNT, DATATYPE, SOURCE, TAG, COMM, STATUS, IERROR)
其中<type> BUF(*);INTEGER COUNT, DATATYPE, SOURCE, TAG, COMM, STATUS(MPI_STATUS_SIZE),IERROR

C:
    int MPI_Recv(void* buf, int count, MPI_Datatype datatype, int source,int tag, MPI_Comm comm, MPI_Status *status)
C++:
    MPI::Comm::Recv(void* buf, int count, const MPI::Datatype& datatype,int source, int tag, MPI::Status& status)

这一组方法是阻塞式的，也就是会导致程序等待直到收到对应消息，所以有可能会发生程序死锁例如

if (myid==0)
{
 MPI_RECV(a,100,MPI_INTEGER,0,10,comm, status)
 MPI_SEND(b,100,MPI_INTEGER,1,11,comm,status)
}
else
{
 MPI_RECV(a,100,MPI_INTEGER,0,11,comm,status)
 MPI_SEND(b,100,MPI_INTEGER,0,10,comm,status)
}

还有Scatter&Gather

Fortran 77/90:
    MPI_SCATTER(SENDBUF, SENDCOUNT, SENDTYPE, RECVBUF, RECVCOUNT, RECVTYPE, ROOT, COMM, IERROR)
其中<type> SENDBUF(*), RECVBUF(*); INTEGER SENDCOUNT, SENDTYPE, RECVCOUNT, RECVTYPE, ROOT, COMM, IERROR
C:
    int MPI_Scatter(void* sendbuf, int sendcount, MPI_Datatype sendtype, void* recvbuf, int recvcount, MPI_Datatype recvtype, int root,MPI_Comm comm)
C++:
    MPI::Comm::Scatter(const void* sendbuf, int sendcount, const MPI::Datatype& sendtype, void* recvbuf, int recvcount,
const MPI::Datatype& recvtype, int root)

Fortran 77/90:
    Call MPI_GATHER(SENDBUF, SENDCOUNT, SENDTYPE, RECVBUF, RECVCOUNT, RECVTYPE, ROOT, COMM, IERROR)
其中<type> SENDBUF(*), RECVBUF(*);INTEGER SENDCOUNT, SENDTYPE, RECVCOUNT, RECVTYPE, ROOT, COMM, IERROR
C:
    int MPI_Gather(void* sendbuf, int sendcount, MPI_Datatype sendtype, void* recvbuf, int recvcount, MPI_Datatype recvtype, int root, MPI_Comm comm)
C++:
    MPI::Comm::Gather(const void* sendbuf, int sendcount, const MPI::Datatype& sendtype, void* recvbuf, int recvcount,
const MPI::Datatype& recvtype, int root)

以下是例子

Fortran的

Program mpitest
  use mpi
  Implicit None
  Integer,Parameter :: row=5000
  Integer,Parameter :: col=5000
  Real(Kind=8) :: mat_a(row,col),mat_b(row,col),mat_c(row,col)
  Integer :: myid,numprocs,rc,ierr,iprovided
  Integer(kind=sp) :: startcol,endcol,colsn
  Integer(kind=sp) :: i,j,k
  Call MPI_INIT(ierr)
  Call MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr)
  Call MPI_COMM_SIZE(MPI_COMM_WORLD,numprocs,ierr)
  mat_a=0.0d0
  mat_b=0.0d0
  mat_c=0.0d0
  colsn=col/numprocs
  startcol=colsn*myid+1
  endcol=colsn*(myid+1)
  If (myid==0) Then
     mat_a=100.0d0
     mat_b=80.0d0
  End If
  Call MPI_SCATTER(mat_a,colsn*row,MPI_DOUBLE_PRECISION,mat_a,colsn*row,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
  Call MPI_SCATTER(mat_b,colsn*row,MPI_DOUBLE_PRECISION,mat_b,colsn*row,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
   mat_c(:,1:colsn)=mat_a(:,1:colsn)-mat_b(:,1:colsn)
  Call MPI_GATHER(mat_c(:,1:colsn),colsn*row,MPI_DOUBLE_PRECISION,mat_c,colsn*row,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)

  Call MPI_FINALIZE(rc)
  Stop
End Program mpitest

C++的

# include <cstdlib>
# include <iostream>
# include <ctime>
# include "mpi.h"

using namespace std;

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

//****************************************************************************80

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

{
    const int row=5000;
    const int col=5000;
    int myid;
    int nprocs;
    int colsn;
    int i,j,k;
     MPI::Init(argc,argv);
    myid=MPI::COMM_WORLD.Get_rank();
    nprocs=MPI::COMM_WORLD.Get_size();
    colsn=row/nprocs;
    double (*mat_a)[col];
    double (*mat_b)[col];
    double (*mat_c)[col];
    double (*tmp_a)[col];
    double (*tmp_b)[col];
    double (*tmp_c)[col];
    tmp_a=new double[colsn][col];
    tmp_b=new double[colsn][col];
    tmp_c=new double[colsn][col];
    if (myid==0)
    {
        mat_a=new double[row][col];
        mat_b=new double[row][col];
        mat_c=new double[row][col];
        for ( i = 0; i < row; i++ )
        {
            for ( j = 0; j < col; j++ )
            {
                mat_a[i][j]=100.0;
            }
        }
        for ( i = 0; i < row; i++ )
        {
            for ( j = 0; j < col; j++ )
            {
                mat_b[i][j]=80.0;
            }
        }
        for ( i = 0; i < row; i++ )
        {
            for ( j = 0; j < col; j++ )
            {
                mat_c[i][j]=0.0;
            }
        }
    }
    else
    {
        mat_a=new double[1][col];
        mat_b=new double[1][col];
        mat_c=new double[1][col];
        for ( i = 0; i < colsn; i++ )
        {
            for ( j = 0; j < col; j++ )
            {
                tmp_a[i][j]=0.0;
            }
        }
        for ( i = 0; i < colsn; i++ )
        {
            for ( j = 0; j < col; j++ )
            {
                tmp_b[i][j]=0.0;
            }
        }
        for ( i = 0; i < colsn; i++ )
        {
            for ( j = 0; j < col; j++ )
            {
                tmp_c[i][j]=0.0;
            }
        }


    }

    MPI::COMM_WORLD.Scatter(mat_a, colsn*col, MPI::DOUBLE, tmp_a, colsn*col,MPI::DOUBLE,0);
    MPI::COMM_WORLD.Scatter(mat_b, colsn*col, MPI::DOUBLE, tmp_b, colsn*col,MPI::DOUBLE,0);

    for (i=0;i<colsn;i++)
    {
         for (j=0;j<col;j++)
         {
             tmp_c[i][j]=tmp_a[i][j]+tmp_b[i][j];
         }
    }
  
    MPI::COMM_WORLD.Gather(tmp_c, colsn*col, MPI::DOUBLE, mat_c, colsn*col,MPI::DOUBLE,0);
    MPI::Finalize();
    return 0;
}

教程

http://www.lam-mpi.org/tutorials/bindings/

http://www.amazon.com/Introduction-Parallel-Computing-Ananth-Grama/dp/0201648652

http://www.megashare.com/3146342