记一次使用MPI完成矩阵乘法程序

最新推荐文章于 2024-07-09 21:51:56 发布

LT_real

最新推荐文章于 2024-07-09 21:51:56 发布

阅读量387

点赞数 5

文章标签：矩阵线性代数

本文链接：https://blog.csdn.net/qq_47814951/article/details/139034527

版权

记录使用MPI完成矩阵乘法程序的过程。

一、MPI安装

在linux系统中，输入命令：

sudo apt-get install mpich

如果要卸载，则输入：

sudo apt-get --purge remove mpich

二、矩阵乘法程序

1.初始化

由于我们要使用MPI传递二维数组，故这里我们创建三个N*N的一维数组A、B、C以便后续数据传输。

    double *A = new double [n*n];
    double *B = new double [n*n];
    double *partA = new double [n*n];
    double *partC = new double [n*n];


    if(myid==0){
        for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                if(i<N && j<N){
                    A[i*n+j] = rand()%10;
                    B[i*n+j] = rand()%10;
                }else{
                    A[i*n+j] = 0;
                    B[i*n+j] = 0;
                }
                
            }
        }
        startWtime = MPI_Wtime();
    }

2.矩阵乘法

普通矩阵乘法只需要3个for循环即可，但考虑到Cache存储，我们可以将一个矩阵转置过来，提高Cache命中率：

void matM(double *A,double *B,double *C,int divideN,int n){//cache
    double mid;
    for(int i=0;i<divideN;i++)for(int j=0;j<i;j++){
        mid=B[i*n+j];
        B[i*n+j]=B[j*n+i];
        B[j*n+i]=mid;
    }
    for(int i=0;i<divideN;i++){
            for(int j=0;j<n;j++){
                C[i*n+j]=0;
                for(int k=0;k<n;k++){
                    C[i*n+j]+=A[i*n+k] * B[j*n+k];
                }
            }
        }
    for(int i=0;i<divideN;i++)for(int j=0;j<i;j++){
        mid=B[i*n+j];
        B[i*n+j]=B[j*n+i];
        B[j*n+i]=mid;
    }
}

3.并行矩阵乘法程序

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <stdbool.h>
#include <iostream>
#include <mpi.h>
#define N 8000

void matM(double *A,double *B,double *C,int divideN,int n){//cache
    double mid;
    for(int i=0;i<divideN;i++)for(int j=0;j<i;j++){
        mid=B[i*n+j];
        B[i*n+j]=B[j*n+i];
        B[j*n+i]=mid;
    }
    for(int i=0;i<divideN;i++){
            for(int j=0;j<n;j++){
                C[i*n+j]=0;
                for(int k=0;k<n;k++){
                    C[i*n+j]+=A[i*n+k] * B[j*n+k];
                }
            }
        }
    for(int i=0;i<divideN;i++)for(int j=0;j<i;j++){
        mid=B[i*n+j];
        B[i*n+j]=B[j*n+i];
        B[j*n+i]=mid;
    }
}

int main(int argc,char *argv[])
{
    double startWtime=0.0,endWtime=0.0;
    int myid,numProcs,jud=0;
    MPI_Init(&argc, &argv );
    MPI_Comm_rank(MPI_COMM_WORLD, &myid );
    MPI_Comm_size(MPI_COMM_WORLD, &numProcs );

    int n=0;
    if (N % numProcs!=0){
        n=N-N%numProcs+numProcs;
    }else{
        n=N;
    }
    int divideN = n/numProcs;

    srand(1);//保持一致性
    double *A = new double [n*n];
    double *B = new double [n*n];
    double *partA = new double [n*n];
    double *partC = new double [n*n];


    if(myid==0){
        for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                if(i<N && j<N){
                    A[i*n+j] = rand()%10;
                    B[i*n+j] = rand()%10;
                }else{
                    A[i*n+j] = 0;
                    B[i*n+j] = 0;
                }
                
            }
        }
        startWtime = MPI_Wtime();
    }
    
        MPI_Scatter(&A[0*n+0],divideN *n,MPI_DOUBLE,&partA[0*n+0],divideN *n,MPI_DOUBLE,0,MPI_COMM_WORLD);
        MPI_Bcast(&B[0*n+0],n*n,MPI_DOUBLE,0,MPI_COMM_WORLD);

        matM(partA,B,partC,divideN,n);
        
    double *C = nullptr;
    if(myid==0) C = new double [n*n];
        MPI_Gather(&partC[0*n+0],divideN *n,MPI_DOUBLE,&C[0*n+0],divideN *n,MPI_DOUBLE,0,MPI_COMM_WORLD);
    if(myid==0){
        endWtime = MPI_Wtime();
        printf("用时: %.3f 秒\n", endWtime-startWtime);
    }

    delete []A;
    delete []B;
    delete []C;
    delete []partA;
    delete []partC;
    MPI_Finalize();
    
    return 0;
}