「并行程序设计」Fox矩阵乘法的MPI实现

已于 2024-07-04 22:41:37 修改
阅读量873 收藏 23
点赞数 24
分类专栏： 并行计算 文章标签： 矩阵 c++ 分布式
于 2024-07-04 22:38:42 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/lk32767/article/details/140190015
版权

目录

算法介绍

用 Fox 矩阵乘法计算 C = A × B C=A\times B C=A×B。

分块

设处理器数量为 p p p，矩阵A、B和C的维度为 n × n n\times n n×n。将A、B分成 p = p × p p=\sqrt p \times \sqrt p p=p ​×p ​ 个方块 A i j A_{ij} Aij​和 B i j B_{ij} Bij​，每个块的大小为 n p × n p \frac{n}{\sqrt p}\times \frac{n}{\sqrt p} p ​n​×p ​n​，并将它们分配给 p × p \sqrt p \times\sqrt p p ​×p ​个处理器。开始时，处理器 P i j P_{ij} Pij​存放块 A i j A_{ij} Aij​和 B i j B_{ij} Bij​，负责计算 C i j C_{ij} Cij​。

算法步骤

  1. 处理器 P i i P_{ii} Pii​将块 A i j A_{ij} Aij​向所在行的其他处理器进行一到多播送；
  2. 各处理器将收到的A块与原有的B块相乘，累加到本地C块中；
  3. 计算完成后，各处理器向上发送本地的B块，再从下方处理器接收新的B块；
  4. 若 A i k A_{ik} Aik​是上一次第 i i i行播送的A块，则本次选择 A i , ( k + 1 ) m o d p A_{i,(k+1)mod \sqrt p} Ai,(k+1)modp ​​向所在行的其他处理器进行播送；
  5. 转步骤2，执行 p − 1 \sqrt p -1 p ​−1次；
  6. 最后由 P 00 P_{00} P00​处理器收集各处理器C块的计算结果。

请添加图片描述请添加图片描述

MPI代码实现

MPI（Massage Passing Interface，消息传递接口）是一种并行编程库，支持C、Fortran等多种编程语言。

本实验在集群上运行MPI程序。

算法过程

  1. 由标识符为0的处理器随机生成 n × n n\times n n×n维待乘矩阵A、B；
  2. 由处理器0将A、B分块，用 MPI_Send() 函数将对应块发送给另外 p × p − 1 \sqrt p \times\sqrt p -1 p ​×p ​−1个处理器，用MPI_Recv()函数接收；
  3. 所有处理器执行 F o x ( ) Fox() Fox()函数，得到本地计算结果C
  4. 最后，所有处理器将计算结果发送给处理器0，由处理器0打印计算结果和计算用时，并与串行矩阵乘法的用时进行比较。

全部代码如下：

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

// mpicc -o fox fox.c -lm   //将fox.c编译为可执行文件fox
// mpirun -np 4 ./fox  //选择线程数为4，并执行fox文件

const int N = 16;//矩阵A, B的维度

void Print_Mat(int *A, int n){//打印矩阵
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            printf("%d  ",A[i * n + j]);
        }
        printf("\n");
    }
    printf("\n");
}

int int_sqrt(int p){//计算p的根号，返回整数结果
    for(int i = 1; i <= p; i++){
        if(i * i == p){
            return i;
        }
    }
    return -1;
}

//获得A的第i,j个n*n方块,储存在a中
void Get_Block(int *A, int *a, int i, int j, int n){
    for(int k = 0; k < n; k++){
        for(int l = 0; l < n; l++){
            a[k * n + l] = A[i * n * N + j * n + k * N + l];
        }
    }
}

//矩阵相乘的串行算法，计算矩阵乘法 a*b，计算结果储存在c中
void Multi_Mat(int *a, int *b, int *c, int n){
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            for(int k = 0; k < n; k++){
                c[i * n + j] += a[i * n + k] * b[k * n + j];
            }
        }
    }
}

//fox矩阵乘法
void Fox(int *a, int *b, int *c, int sp, int n, int myrank){
    int *temp_a = malloc(n * n * sizeof(int));//用来接收a
    int *temp_b = malloc(n * n * sizeof(int));//用来接收b

    //将处理器记为P_ij
    int j = myrank % sp;
    int i = (myrank - j) / sp;

    //对方块c初始化
    for(int k = 0; k < n * n; k++){
        c[k]=0;
    }
    int senddest = 0, recvdest = 0;
    //fox循环
    for(int round = 0; round < sp; round++){
        if(i == (j + round) % sp){//选出本轮在i行广播的A_ij
            for(int k = 0; k < sp; k++){
                if(k != j){
                    MPI_Send(a, n * n, MPI_INT, i * sp + k, (round + 1) * sp, MPI_COMM_WORLD);
                }
                else{
                    for(int l = 0; l < n * n; l++){
                        temp_a[l]=a[l];
                    }
                }
            }
        }
        else{
            if(i - round < 0){
                recvdest = i * sp + i - round + sp;
            }else{
                recvdest = i * sp + i -round;
            }
            MPI_Recv(temp_a, n * n, MPI_INT, recvdest, (round + 1) * sp, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        }

        //做矩阵乘法
        Multi_Mat(temp_a, b, c, n);

        if(round == sp - 1)break;
        
        //将块B_ij上移
        senddest = i > 0 ? i - 1:i - 1 + sp;
        recvdest = i < sp - 1 ? i + 1:i + 1 - sp;
        MPI_Sendrecv(b, n * n, MPI_INT, senddest * sp + j, 0, 
                    temp_b, n * n, MPI_INT, recvdest * sp + j, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        for(int i = 0; i < n * n; i++){
            b[i] = temp_b[i];
        }
    }
}

//将C的第i,j个n*n方块用c赋值
void Copy_Block(int *c, int *C, int i, int j, int n){
    for(int k = 0; k < n; k++){
        for(int l = 0; l < n; l++){
            C[i * n * N + j * n + k * N + l] = c[k * n + l];
        }
    }
}

// processor0 接收其他处理器的计算结果
void Recv_Block(int *c, int *C, int sp, int n){
    for(int i = 0; i < sp; i++){
        for(int j = 0; j < sp; j++){
            if(i + j != 0) 
                MPI_Recv(c, n * n, MPI_INT, i * sp + j, 3, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            Copy_Block(c, C, i, j, n);
        }
    }
}

int main(int argc, char **argv)
{
    MPI_Init(&argc, &argv);

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

    int *A, *B, *C;
    if(myrank == 0){
        // 由processor0 随机生成N*N矩阵A、B，和空矩阵C
        srand(time(NULL));
        A = (int *)malloc(N * N * sizeof(int));
        B = (int *)malloc(N * N * sizeof(int));
        C = (int *)malloc(N * N * sizeof(int));
        for(int i = 0; i < N * N; i++){
            A[i] = rand() % 10;
            B[i] = rand() % 10;
        }
        //Print_Mat(A, N);
    }

    int sp = int_sqrt(p);
    if(sp == -1 && myrank == 0){
        printf("sp = %d, 处理器数不是完全平方数\n",sp);
        return(0);
    }
    int n = N / sp;//方块的维度
    if(sp * n != N && myrank == 0){
        printf("处理器数不能均分矩阵\n");
        return(0);
    }
    //printf("n = %d\n", n);
    int *a = malloc(n * n * sizeof(int));//a, b, c储存方块
    int *b = malloc(n * n * sizeof(int));
    int *c = malloc(n * n * sizeof(int));

    //计时开始
    MPI_Barrier(MPI_COMM_WORLD);
    double start_fox = MPI_Wtime();

    if(myrank == 0){
        //由processor0把A，B的n*n方块发送给另外sp*sp-1个处理器
        for(int i = 0; i < sp; i++){
            for(int j = 0; j < sp; j++){
                if(i == 0 && j == 0) continue;
                Get_Block(A, a, i, j, n);
                Get_Block(B, b, i, j, n);
                //Print_Mat(a, n);
                MPI_Send(a, n * n, MPI_INT, i * sp + j, 'a', MPI_COMM_WORLD);
                MPI_Send(b, n * n, MPI_INT, i * sp + j, 'b', MPI_COMM_WORLD);
            }
        }
        // processor0 储存方块A_00, B_00
        Get_Block(A, a, 0, 0, n);
        Get_Block(B, b, 0, 0, n);
    }
    else{// 其他处理器接收方块A_ij, B_ij
        MPI_Recv(a, n * n, MPI_INT, 0, 'a', MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        MPI_Recv(b, n * n, MPI_INT, 0, 'b', MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    }
    
    // Fox矩阵乘法
    Fox(a, b, c, sp, n, myrank);

    if(myrank != 0){//其他处理器向 processor0 发送计算结果
        MPI_Send(c, n * n, MPI_INT, 0, 3, MPI_COMM_WORLD);
    }
    else{//processor0 接收方块结果，并赋值给矩阵C
        Recv_Block(c, C, sp, n);
    }

    //计时结束
    double finish_fox = MPI_Wtime();

    // 由 processor0 打印计算结果
    if (myrank == 0)
    {   
        // 检验fox乘法结果正确性
        #if 0
        printf("Matrix A = \n");Print_Mat(A, N);
        printf("Matrix B = \n");Print_Mat(B, N);
        printf("fox计算结果 C = \n");Print_Mat(C, N);
        #endif
        printf("fox并行乘法耗时: %f\n", finish_fox - start_fox);
        
        //归零
        for(int i=0; i < N * N; i++){
            C[i]=0;
        }

        //串行乘法
        double start =  MPI_Wtime();
        //Multi_Mat(A, B, C, N);
        double finish = MPI_Wtime();

        #if 0
        printf("串行计算结果 C = \n");Print_Mat(C, N);
        #endif
        printf("串行乘法耗时: %f\n", finish - start);

        //计算加速比
        double rate = (finish - start)/(finish_fox - start_fox);
        printf("并行加速比为: %f\n", rate);
    }

    MPI_Finalize();
    return 0;
}

实验结果

  1. 取矩阵维度 n = 16 n=16 n=16，处理器数量 p = 4 p=4 p=4，打印计算结果。

请添加图片描述

  1. 可以改变矩阵维度和处理器数量，例如依次取 n = 16 , 128 , 256 , 512 , 1024 n = 16,128,256, 512, 1024 n=16,128,256,512,1024和 p = 4 , 16 , 64 p=4,16,64 p=4,16,64，比较并行加速比的变化。
怀帝阍而不见
关注
  • 24
    点赞
  • 踩
  • 23
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
MPI群通信与矩阵乘法的Fox算法实现
Jarod 的半个程序员
12-04 6119
原本以为，MPI天生只能在Linux上运行。但这次却发现了Intel MPI Library 这个好用的东西。基本不需要设置，安上之后，用自己能登录windows的帐号和密码注册就行了。虽然不是局域网上的机器，但也可以让我的双核CPU达到100%（平时开个Matlab什么的都才是50%，软件优化真是关键啊）。FOX算法有一些恶心的要求：输入的矩阵必须是方阵，而且进程必须为平方数，方阵必须能均匀
矩阵乘法 mpi c语言,C语言基于MPI并行计算矩阵的乘法
weixin_35997546的博客
05-18 1428
浏览量：1,023(在苹果系统下，如果文章中的图片不能正常显示，请升级Safari浏览器到最新版本，或者使用Chrome、Firefox浏览器打开。)矩阵的计算是并行计算里的一个很重要的问题。矩阵是一种网格化的数据，是一组同类型数值的集合，矩阵的出现，使得代数系统更完善，对各种实际问题的求解产生了巨大的作用。但是其庞大的计算量往往令人生畏，稍微大一点的矩阵，计算就变得非常繁琐，不仅仅是人，而且串行...
矩阵相乘FOX并行算法
12-10
矩阵相乘FOX并行算法，矩阵相乘FOX并行算法
FOX算法的MPI实现
cancumt的博客
01-11 4940
FOX算法的MPI实现算法描述如下： 将待相乘的矩阵A和B分成p个方块Ai,j和Bi,j(0≤j,i≤√p-1)，每块大小为（n/√p）×（n/√p），并将他们分配给√p×/√p，个处理器。开始时处理器Pi,j存放有Ai,j和Bi,j，并负责计算块Ci,j。然后Fox算法执行以下√p次迭代，即可完成： ①选中对角块Ai,j，并将其向所在行的√p-1个处理器进行一到多播送； ②各处理器将所收到
MPI通信子和网格拓扑结构的应用——fox矩阵相乘算法
qq_43571494的博客
12-10 382
/* File: mpi_mat_vect_mult.c * * Purpose: Implement parallel matrix-vector multiplication using * one-dimensional arrays to store the vectors and the * matrix. Vectors use block distributions and the * matrix is dis
矩阵相乘的FOX并行实现
11-27
矩阵相乘的FOX并行实现，这是矩阵相乘的另外一种实现方法。
mpi大规模矩阵乘法C语言,C语言基于MPI并行计算矩阵的乘法
weixin_39753260的博客
05-24 401
浏览量：1,030(在苹果系统下，如果文章中的图片不能正常显示，请升级Safari浏览器到最新版本，或者使用Chrome、Firefox浏览器打开。)矩阵的计算是并行计算里的一个很重要的问题。矩阵是一种网格化的数据，是一组同类型数值的集合，矩阵的出现，使得代数系统更完善，对各种实际问题的求解产生了巨大的作用。但是其庞大的计算量往往令人生畏，稍微大一点的矩阵，计算就变得非常繁琐，不仅仅是人，而且串行...
矩阵乘法mpi实现 并行运算
05-27
矩阵乘法mpi实现 并行运算 在linux下成功运行 使用mpicc -o 编译 使用mpirun命令运行
矩阵乘法MPI并行程序报告.pdf
06-06
并行程序设计部分，核心在于MPI的使用。`MPI_Init(&argc, &argv)`初始化MPI环境，`MPI_Finalize()`则在程序结束时关闭环境。`MPI_Comm_rank`获取当前进程的ID，`MPI_Comm_size`得到进程组的总进程数。例如，`MPI_...
MPI并行计算 矩阵乘法
05-31
MPI并行计算 矩阵乘法 如果是10000个线程的话可能会出问题
MPI实现矩阵乘法
12-24
通过分块利用MPI通讯实现矩阵乘法的并行计算。
matvec-mpi:使用 MPI 的稀疏矩阵向量并行乘法实现
06-19
##稀疏矩阵向量乘法与 MPI 并行###Design 使用 MPI 并行化稀疏矩阵向量乘法： 在步骤 1 中使用一维行分解读取文件并将数据分发到所有处理器，这需要 O(n) 然后 O(nnz) 其中 n 是行数，nnz 是矩阵。 矩阵 A 数据以 ...
矩阵相乘FOX算法的mpi实现
koastal的博客
01-20 2618
#include "mpi.h" #include #include #include const int root_id = 0; const int max_procs_size = 16; int main(int argc,char *argv[]) { double start_time, end_time, time; int procs_id, procs_
矩阵相乘
叨叨令
08-20 1279
import java.util.Scanner; public class Test19 { /* 题目描述  如果A是个x行y列的矩阵，B是个y行z列的矩阵，把A和B相乘，其结果将是另一个x行z列的 矩阵C。这个矩阵的每个元素是由下面的公式决定的：    原型： voidmatrix_multiply(int *m1,int *m2,int *r, int x, i...
矩阵乘法的并行算法
热门推荐
realxie的专栏
02-15 2万+
设两个矩阵A和B，大小分别为M * N 和 N * P， 如果C = A * B， 则C的大小为M * P。 矩阵算法的算法表示，伪代码如下： for (i = 0; i < M; ++i){ for (j = 0; j < P; ++j){ C[i][j] = 0; for (k = 0; k < N; ++k){ C[i][j] += A[i][k] * B[k][j];
BC131 矩阵相等判定（c语言）
最新发布
scy2429828663的博客
08-23 296
【代码】BC131 矩阵相等判定（c语言）
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值