MPI矩陣乘法

去年學習了平行計算，接觸了MPI、Pthreads和OpenMP等常用的並行方法實現了並行的矩陣乘法，本章在此總結一下MPI的矩陣乘法使用。

• 使用簡單的MPI_Send和MPI_Recv實現
• 使用較高階的MPI_Scatter和MPI_Gather實現

MPI_Send和MPI_Recv實現

#include<stdio.h>
#include<stdlib.h>
#include<mpi.h>
#include<time.h>

int main(int argc,char *argv[])
{
    double
 
 start, stop;
    int i, j, k, l;
    int *a, *b, *c, *buffer, *ans;
    int size = 1000;
    int rank, numprocs, line;

    MPI_Init(NULL,NULL);//MPI Initialize
    MPI_Comm_rank(MPI_COMM_WORLD,&rank);//獲得當前程序號
    MPI_Comm_size(MPI_COMM_WORLD,&numprocs);//獲得程序個數

    line = size/numprocs;//將資料分為(程序數)個塊,主程序也要處理資料
 

    a = (int*)malloc(sizeof(int)*size*size);
    b = (int*)malloc(sizeof(int)*size*size);
    c = (int*)malloc(sizeof(int)*size*size);
    //快取大小大於等於要處理的資料大小，大於時只需關注實際資料那部分
    buffer = (int*)malloc(sizeof(int)*size*line);//資料分組大小
    ans = (int*)malloc(sizeof(int)*size*line);//儲存資料塊計算的結果

    //主程序對矩陣賦初值，並將矩陣N廣播到各程序,將矩陣M分組廣播到各程序
 

    if (rank==0)
    {
        //從檔案中讀入矩陣
        FILE *fp;

        fp=fopen("a.txt","r");//開啟檔案
        start = MPI_Wtime();
        for(i=0;i<1000;i++) //讀資料
            for(j=0;j<1000;j++)
                fscanf(fp,"%d",&a[i*size+j]);  
        fclose(fp);//關閉檔案

        fp=fopen("b.txt","r");

        for(i=0;i<1000;i++) 
            for(j=0;j<1000;j++)
                fscanf(fp,"%d",&b[i*size+j]);  
        fclose(fp);
        //將矩陣N傳送給其他從程序
        for (i=1;i<numprocs;i++)
        {
                MPI_Send(b,size*size,MPI_INT,i,0,MPI_COMM_WORLD);
        }
        //依次將a的各行傳送給各從程序
        for (l=1; l<numprocs; l++)
        {
            MPI_Send(a+(l-1)*line*size,size*line,MPI_INT,l,1,MPI_COMM_WORLD);
        }
        //接收從程序計算的結果
        for (k=1;k<numprocs;k++)
        {
            MPI_Recv(ans,line*size,MPI_INT,k,3,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
            //將結果傳遞給陣列c
            for (i=0;i<line;i++)
            {
                for (j=0;j<size;j++)
                {
                    c[((k-1)*line+i)*size+j] = ans[i*size+j];
                }

            }
        }
        //計算a剩下的資料
        for (i=(numprocs-1)*line;i<size;i++)
        {
            for (j=0;j<size;j++)
            {
                int temp=0;
                for (k=0;k<size;k++)
                    temp += a[i*size+k]*b[k*size+j];
                c[i*size+j] = temp;
            }
        }

        fp=fopen("c.txt","w");
        for(i=0; i<size; i++){
            for(j=0; j<size; j++)
                fprintf(fp,"%d ",c[i*size+j]);
            fputc('\n',fp);
        }
        fclose(fp);
        //結果測試
        //統計時間
        stop = MPI_Wtime();

        printf("rank:%d time:%lfs\n",rank,stop-start); 

        free(a);
        free(b);
        free(c);
        free(buffer);
        free(ans);
    }

    //其他程序接收資料，計算結果後，傳送給主程序
    else
    {
        //接收廣播的資料(矩陣b)
        MPI_Recv(b,size*size,MPI_INT,0,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE);

        MPI_Recv(buffer,size*line,MPI_INT,0,1,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
        //計算乘積結果，並將結果傳送給主程序
        for (i=0;i<line;i++)
        {
            for (j=0;j<size;j++)
            {
                int temp=0;
                for(k=0;k<size;k++)
                    temp += buffer[i*size+k]*b[k*size+j];
                ans[i*size+j]=temp;
            }
        }
        //將計算結果傳送給主程序
        MPI_Send(ans,line*size,MPI_INT,0,3,MPI_COMM_WORLD);
    }

    MPI_Finalize();//結束

    return 0;
}

MPI_Scatter和MPI_Gather實現

#include<stdio.h>
#include<mpi.h>
#include <malloc.h>
#define M 1000
#define N 1000
int main()
{
int my_rank;/*My process rank*/
int comm_sz;/*Number of processes*/
int local_M;
int i,j,k;
double start,finish;/*timer*/
int tem;
//初始化MPI
MPI_Init(NULL,NULL);
MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
MPI_Comm_size(MPI_COMM_WORLD,&comm_sz);
//每個矩陣分配到的行數
local_M=M/comm_sz;
//分配到每個程序的矩陣
int *local_Matrix_one=(int*)malloc(local_M*N*sizeof(int));
//定義兩個矩陣
int *Matrix_one=NULL;
int *Matrix_two=(int*)malloc(M*N*sizeof(int));
//每個程序裡的結果矩陣
int *local_result=(int*)malloc(local_M*N*sizeof(int));
//結果矩陣
int *result_Matrix=NULL;
if(my_rank==0)
{
//printf("process %d of %d\n",my_rank,comm_sz);
FILE * fp;
//讀取第一個矩陣
Matrix_one=(int*)malloc(M*N*sizeof(int));
//Matrix_one[M][N]={0};
fp=fopen("a.txt","r");
for(i=0;i<M;i++)
{
for(j=0;j<N;j++)
fscanf(fp,"%d ",&Matrix_one[i*N+j]);
fscanf(fp,"\n");
}
fclose(fp);
/*for(i=0;i<M;i++)
{
for(j=0;j<N;j++)
   printf("%d ",Matrix_one[i*N+j]);
   printf("\n");
}*/
//讀取第二個矩陣
start=MPI_Wtime();
fp=fopen("b.txt","r");
for(j=0;j<N;j++)
{
for(i=0;i<M;i++)
fscanf(fp,"%d ",&Matrix_two[i*N+j]);
fscanf(fp,"\n");
}
fclose(fp);
/*for(j=0;j<N;j++)
{
for(i=0;i<M;i++)
   printf("%d ",Matrix_two[i*M+j]);
   printf("\n");
}*/
//資料分發
MPI_Scatter(Matrix_one,local_M*N,MPI_INT,local_Matrix_one,local_M*N,MPI_INT,0,MPI_COMM_WORLD);

//資料廣播
MPI_Bcast(Matrix_two,M*N,MPI_INT,0,MPI_COMM_WORLD);
//計算local結果
for(i=0;i<local_M;i++)
for(j=0;j<M;j++){
tem=0;
for(k=0;k<N;k++)
tem +=local_Matrix_one[i*M+k]*Matrix_two[j*M+k];
local_result[i*M+j]=tem;
}
free(local_Matrix_one);
result_Matrix=(int*)malloc(M*N*sizeof(int));
//結果聚集
MPI_Gather(local_result,local_M*N,MPI_INT,result_Matrix,local_M*N,MPI_INT,0,MPI_COMM_WORLD);
//剩餘行處理（處理不能整除的情況）
int rest=M%comm_sz;
if(rest!=0)
for(i=M-rest-1;i<M;i++)
for(j=0;j<M;j++){
tem=0;
for(k=0;k<N;k++)
tem +=Matrix_one[i*M+k]*Matrix_two[j*M+k];
result_Matrix[i*M+j]=tem;
}
finish=MPI_Wtime();
free(Matrix_one);
free(Matrix_two);
free(local_result);

printf("Proc %d > Elapsed time = %e seconds\n",my_rank,finish-start);
//將結果寫入檔案
fp=fopen("c.txt","w");
for(i=0;i<M;i++)
{
for(j=0;j<N;j++)
fprintf(fp,"%d ",result_Matrix[i*N+j]);
fscanf(fp,"\n");
}
fclose(fp);
/*for(i=0;i<local_M;i++)
{
for(j=0;j<N;j++)
   printf("%d ",local_result[i*N+j]);
   printf("\n");
}*/
}
else{
//printf("process %d of %d\n",my_rank,comm_sz);
//資料分發
MPI_Scatter(Matrix_one,local_M*N,MPI_INT,local_Matrix_one,local_M*N,MPI_INT,0,MPI_COMM_WORLD);
//資料廣播
MPI_Bcast(Matrix_two,M*N,MPI_INT,0,MPI_COMM_WORLD);
//計算local結果
for(i=0;i<local_M;i++)
for(j=0;j<M;j++){
tem=0;
for(k=0;k<N;k++)
tem +=local_Matrix_one[i*M+k]*Matrix_two[j*M+k];
local_result[i*M+j]=tem;
}
free(local_Matrix_one);
free(Matrix_two);
//結果聚集
MPI_Gather(local_result,local_M*N,MPI_INT,result_Matrix,local_M*N,MPI_INT,0,MPI_COMM_WORLD);
free(local_result);
//printf("%d %d\n",local_M,my_rank);
}
MPI_Finalize();
return 0;
}

結果分析

① 執行時間分析：
並行時，隨著程序數目的增多，平行計算的時間越來越短；當達到一定的程序數時，執行時間小到最小值；然後再隨著程序數的增多，執行時間反而越來越長。
② 加速比分析：
隨著程序數的增大，加速比也是逐漸增大到最大值；再隨著程序數的增大，加速比逐漸減小。
③ 執行效率分析：
隨著程序數的增大，程式執行效率不斷降低
④ 原因分析：
MPI並行程式的測試平臺為Intel Core i5 CPU，為雙核CPU，即在一個處理器上整合兩個運算核心，提高了運算效率，因此會比序列的執行時間要短。由於一個程序只能在一個核上執行，因此只能有兩個程序並行執行，又因為多程序執行在兩個CPU上，會有程序切換等操作，所以才會出現程序數增加而執行時間增加的情況。