MPI Matrix Multiplication

                                                           

                                                                 

有两个矩阵A和B，计算两矩阵的乘积存于C中；

方法一：每个进程分配到矩阵A的部分行，并且有矩阵B的所有元素，每个进程计算矩阵A的部分行与矩阵B的乘积，得到矩阵C的部分行，所有进程计算完毕后，将结果传递给根进程进行汇总，得到最后的结果；

 

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


int main(int argc, char**argv){
  int size, rank, i, j, k;
  int N, len, indexC = 0;
  int *A, *B, *C, *localA, *localC;//原矩阵和部分矩阵
  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &size);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  if(rank == 0){
    printf("the Matrix of N:");
    scanf("%d", &N);//输入矩阵的阶数
    C = (int*)malloc(sizeof(int)*N*N);
    A = (int*)malloc(sizeof(int)*N*N);
    for(i = 0; i < N*N; i++){
      srand((unsigned int)i);
      A[i] = rand() % 10;
    }//赋值
  }
  MPI_Bcast(&N, 1, MPI_INT, 0, MPI_COMM_WORLD);
  B = (int*)malloc(sizeof(int)*N*N);
  for(i = 0; i < N*N; i++){
    srand((unsigned int)i);
    B[i] = rand() % 10;
  }
  len = N / size;//每个进程所得到的矩阵A的行数
  localA = (int*)malloc(sizeof(int)*len*N);
  localC = (int*)malloc(sizeof(int)*len*N);
  MPI_Scatter(A, len*N, MPI_INT, localA, len*N, MPI_INT, 0, MPI_COMM_WORLD);
  for(i = 0; i < len; i++){
    for(j = 0; j < N; j++){
      localC[indexC] = 0;
      for(k = 0; k < N; k++){
        localC[indexC] += localA[N*i+k] * B[N*k+j];
      }
      indexC++;
    }
  }//计算矩阵C的部分行
  MPI_Gather(localC, len*N, MPI_INT, C, len*N, MPI_INT, 0, MPI_COMM_WORLD);
  if(rank == 0){
    printf("the Matrix of A:\n");
    for(i = 0; i < N*N; i++)
      printf("%3d%s", A[i], (i+1) % N == 0 ? "\n" : " ");
    printf("the Matrix of B;\n");
    for(i = 0; i < N*N; i++)
      printf("%3d%s", B[i], (i+1) % N == 0 ? "\n" : " ");
    printf("the Matrix of C:\n");
    for(i = 0; i < N*N; i++)
      printf("%3d%s", C[i], (i+1) % N == 0 ? "\n" : " ");
  }
  MPI_Finalize();
  return 0;
}

 

方法二：每个进程分配到矩阵A的部分行和矩阵B的部分列，每一次的计算都将得到矩阵C的部分元素；每次计算完之后，一次将每个进程存储的矩阵A的部分行传递给下一个进程，再进行同样的计算，循环的次数为进程的数量；所有进程计算完毕，将结果发送给跟进程，得到最后的结果；

#include <stdio.h>
#include <mpi.h>
#include <malloc.h>
#include <stdlib.h>
#include <string.h>


int main(int argc, char**argv){
  int size, rank, i, j, k;
  int *A, *B, *C, *localA, *localB, *localC, *temp;
  int N, len, circle = 0, index, cnt;//circle为循环的次数
  int source, dest;
  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &size);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  if(rank == 0){
    printf("the Matrix of N:");
    scanf("%d", &N);
    A  = (int*)malloc(sizeof(int)*N*N);
    B  = (int*)malloc(sizeof(int)*N*N);
    C  = (int*)malloc(sizeof(int)*N*N);
    for(i = 0; i < N*N; i++){
      srand((unsigned int)i);
      A[i] = rand() % 10;
      srand((unsigned int)i+1);
      B[i] = rand() % 10;
    }
    /*
    for(i = 0; i < N; i++){
      for(j = 0; j<=i; j++){
        k = B[N*i+j];
          B[N*i+j] = B[N*j+i];
          B[N*j+i] = k;
      }
    }
    */
  }
  MPI_Bcast(&N, 1, MPI_INT, 0, MPI_COMM_WORLD);
  len = N / size;//每个进程分配的矩阵A的行数
  localA = (int*)malloc(sizeof(int)*len*N);
  localB = (int*)malloc(sizeof(int)*len*N);
  localC = (int*)malloc(sizeof(int)*len*N);
  temp   = (int*)malloc(sizeof(int)*len*N);
  MPI_Scatter(A, len*N, MPI_INT, localA, len*N, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Scatter(B, len*N, MPI_INT, localB, len*N, MPI_INT, 0, MPI_COMM_WORLD);
  index = rank;//指明当前进程实际所在进程号
  //cnt = 0;
  dest = rank==0 ? size-1 : rank-1;//发送的目标进程号
  source = (rank+1) % size;//接收的目标进程号
  do{
    for(i = 0; i < len; i++){
      for(j = 0; j < len; j++){
        cnt = j*N + (index*len+i);
        localC[cnt] = 0;
        for(k = 0; k <N; k++)
          localC[cnt] += localA[N*i+k] * localB[N*j+k];
      }
      //cnt++;
    }//计算矩阵C中的部分元素

    //一次传递矩阵A中的数据
    if(rank == 0){
      MPI_Ssend(localA, len*N, MPI_INT, dest, 99, MPI_COMM_WORLD);
      MPI_Recv(localA, len*N, MPI_INT, source, 99, MPI_COMM_WORLD,                    MPI_STATUS_IGNORE);
      /*printf("Process %d:\n", rank);
      for(i = 0; i < len*N; i++)
        printf("%d%s", localA[i], (i+1)%N == 0?"\n":" ");*/
      } else{
        memcpy(temp, localA, sizeof(int)*len*N);
        MPI_Recv(localA, len*N, MPI_INT, source, 99, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
      MPI_Ssend(temp, len*N, MPI_INT, dest, 99, MPI_COMM_WORLD);
      /*printf("Process %d:\n", rank);
      for(i = 0; i < len*N; i++)
        printf("%d%s", localA[i], (i+1)%N == 0?"\n":" ");*/
    }
    circle++;
    index = (index+1) % size;
  } while(circle < size);
  MPI_Gather(localC, len*N, MPI_INT, C, len*N, MPI_INT, 0, MPI_COMM_WORLD);
  if(rank == 0){
    printf("the Matrix of A:\n");
    for(i = 0; i < N*N; i++)
      printf("%3d%s", A[i], (i+1) % N == 0 ? "\n" : " ");
    printf("the Matirx of B:\n");
    for(i = 0; i < N; i++)
      for(j = 0; j < N; j++)
        printf("%3d%s", B[N*j+i], (j+1) == N ? "\n" : " ");//矩阵B按列存
    printf("the Matrix of C:\n");
    for(i = 0; i < N; i++)
      for(j = 0; j < N; j++)
        printf("%3d%s", C[N*j+i], (j+1) == N ? "\n" : " "); 
  }
  MPI_Finalize();
  return 0;
}