mpi学习(2)--MPI实现一般矩阵乘法

 

思路：

有两个矩阵A,B。从0号进程读入A,B,并且对B做转置（转置算法不太好，这里做转置主要是为了一会计算方便可以提高CACHE的命中率，但不清楚是否以转置浪费的时间为代价是否合理），之后对A矩阵按照按块分割，发送到各个进程。这里调用了MPI_Scatter，因此矩阵大小要可以被进程数整除。然后对每个进程的私有变量做计算。最后由0号进程回收各个部分的值。

代码：github_mpi_multiply

/*
 * =====================================================================================
 *
 *       Filename:  multiply.c
 *
 *    Description:  矩阵相乘
 *
 *        Version:  1.0
 *        Created:  2015年08月08日 14时54分19秒
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  shiyan (), shiyan233@hotmail.com
 *   Organization:  
 *
 * =====================================================================================
 */

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

/*-----------------------------------------------------------------------------
 *  n*m矩阵,n行m列
 *-----------------------------------------------------------------------------*/
#define N 4
#define M 4

void Transpose(double *b,double *tmp ,int n, int row, int col);
/* 
 * ===  FUNCTION  ======================================================================
 *         Name:  Read_vec
 *  Description:  
 * =====================================================================================
 */
    void
Read_vec (double *local_a,
        double *tmp,
         int local_n,
         int n,
         int myrank,
         MPI_Comm comm)
{
    double *a = NULL;
    double *b = NULL;
    if (myrank == 0)
    {
        srand(time(NULL));
        a = malloc( n*sizeof(double) );
        int i = 0;
        printf("a is\n");
        for (i=0; i<n;i++)
        {
            a[i] = rand()%10;
            if((i+1)%M == 0)
                printf("%lf\n", a[i]);
            else printf("%lf, ", a[i]);
        }
        printf("\nb is \n");    
        b = malloc( n*sizeof(double) );
        i = 0;
        for (i=0; i<n;i++)
        {
            b[i] = rand()%10 ;
            if((i+1) % N == 0)
                printf("%lf\n",b[i]);
            else printf("%lf, ", b[i]);
        }
        putchar('\n');
        Transpose(b, tmp, n, M, N);//将b转置之后放入tmp中
        MPI_Bcast(tmp, M*N,MPI_DOUBLE, 0,MPI_COMM_WORLD);
        MPI_Scatter(a, local_n, MPI_DOUBLE, local_a, local_n, MPI_DOUBLE, 0, comm);
        free(a);
        free(b);    
    }else{
        MPI_Bcast(tmp, N*M, MPI_DOUBLE, 0,MPI_COMM_WORLD);
        MPI_Scatter(a, local_n, MPI_DOUBLE, local_a, local_n, MPI_DOUBLE, 0, comm);    
    }    
}        /* -----  end of function Read_vec  ----- */

/* 
 * ===  FUNCTION  ======================================================================
 *         Name:  Gather_vec
 *  Description:  
 * =====================================================================================
 */
    void
Gather_vec (
        double *local_c,
        int local_n,
        int n,
        int myrank,
        MPI_Comm comm
        )
{
    double *t = NULL;
    int i = 0;
    if (myrank == 0){
        t = malloc(N*N*sizeof(double));
        MPI_Gather(local_c, N, MPI_DOUBLE, t, N, MPI_DOUBLE, 0, comm);
        printf("\nresult:\n");
        for(i=0;i<N*N;i++)if((i+1)%4 == 0){
            printf("%lf\n", t[i]);
        }else printf("%lf ", t[i]);
        putchar('\n');
        free(t);
    }else{
        MPI_Gather(local_c, N, MPI_DOUBLE, t, N, MPI_DOUBLE, 0, comm);
    }
}        /* -----  end of function Gather_vec  ----- */
/* 
 * ===  FUNCTION  ======================================================================
 *         Name:  transpose
 *  Description:  空间O(n)的转置
 * =====================================================================================
 */
void Transpose(double *b,double *tmp ,int n, int row, int col)
{
    int k = 0;
    for (k=0;k<n;k++)
    {
        tmp[k] = b[ (k%row)*col + k/row ];
    }
}        /* -----  end of function transpose  ----- */



/* 
 * ===  FUNCTION  ======================================================================
 *         Name:  Multiply
 *  Description:  
 * =====================================================================================
 */
    void
Multiply ()
{
    MPI_Init(NULL,NULL);

    int myrank;
    int size;
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

    int local_n = (M*N)/size;
    double *local_a = malloc(local_n*sizeof(double));
    double *local_b = malloc(M*N*sizeof(double));
    double *local_c = malloc(N*sizeof(double));
    Read_vec(local_a, local_b, local_n, (M*N), myrank, MPI_COMM_WORLD);
    int i = 0;
    int j = 0;
    double temp = 0;
    for (i=0;i<local_n;i++)
        local_c[i] = 0.0;
    for(i =0 ;i<N; i++){
        temp = 0;
        for(j=0;j<local_n;j++)
            temp += local_a[j]*local_b[i*M+j];
        local_c[i] = temp;
    }
    Gather_vec(local_c, local_n, M*N, myrank, MPI_COMM_WORLD);    
    free(local_a);
    free(local_b);
    free(local_c);
    MPI_Finalize();
}        /* -----  end of function Multiply  ----- */
/* 
 * ===  FUNCTION  ======================================================================
 *         Name:  main
 *  Description:  
 * =====================================================================================
 */
    int
main ( int argc, char *argv[] )
{
    Multiply();
    return EXIT_SUCCESS;
}                /* ----------  end of function main  ---------- */

其中使用的几个函数的参数及作用：

1.MPI_Bcast：MPI的广播函数，作用是从源进程向在通信子内的各个进程发送数据。

MPI_Bcast(void *sendbuf, int size, MPI_Type, int src, MPI_Comm);

第一个参数是要发送的数据，第二个为大小，第三个是MPI类型，第四个是源进程号， 第五个是通信子

2.MPI_Scatter(void *send_buf, int send_size, MPI_Type, void *receive_buf,int receive_size,MPI_Type, int src, MPI_Comm )

sendbuf,发送数据，send_size发送的大小，其中这个发送大小是指发送给每个进程的大小，并不是总的大小。

另外之所以在每个进程都执行了MPI_Scatter是因为这样可以降低复杂度，因为每个进程调用可以将复杂度从O(N)降低到O(logn)，对于MPI_Bcast也是同样道理。

3.MPI_Gather(void *send_buf, int send_size, MPI_Type, void *receive_buf, int receive_size, MPI_Type, int src, MPI_Comm)

将通信子下每个进程的都接收到src进程中，注意接受的大小一定要大于等于发送的大小。

tip:注意内存泄漏问题。