并行计算openMPI程序设计（二）

并行计算编程学习笔记（名词简介）

2.1设计通讯需要考虑的因素

在设计程序任务之间的通讯时，有大量的重要因素需要考虑：

• 通讯开销： 1）任务间通讯几乎总是意味着开销。2）而可以用于计算的机器周期以及资源会转而用于对数据的封装和传输。3）频繁的通讯也需要任务之间的同步，这有可能会导致任务花费时间等待而不是执行。4）竞争通讯流量可能使可用的网络带宽饱和，从而进一步加剧性能问题。

• 延迟 vs. 带宽： 1）延迟 指的是从A点到B点发送最小量的信息所需要花费的时间，通常以毫秒计。2）带宽 指的是单位时间内可以传输的数据总量，通常以M/S或者G/S来计。3）发送大量的短消息可能会导致延迟成为通讯的主要开销。通常情况下将大量小信息封装成为大消息会更加有效，从而提高通讯带宽的利用效率。

• 通讯可见性： 1）在消息传递模型中，通讯往往是显式和可见的，并且在编程者的控制之下。2）在数据并行模型中，通讯对编程者来说往往是透明的，尤其是在分布式内存架构中。编程者往往甚至不能明确知道任务之间的通讯是如何完成的。

• 同步 vs. 异步通讯： 1) 同步通讯需要共享数据的任务之间某种意义上的“握手”。这既可以由编程者显式地指定，也可以在底层被隐式地实现而不为编程者所知。2）同步通讯业常常被称为“阻塞通讯”，因为一些任务必须等待直到它们和其它任务之间的通讯完成。3）异步通讯允许任务之间独立地传输数据。例如任务1可以准备并且发送消息给任务2，然后立即开始做其它工作，它并不关心任务2什么时候真正受到数据。4）异步通讯也常常被称为“非阻塞通讯”，因为在通讯发生的过程中，任务还可以完成其它工作。5）在计算和通讯自由转换是异步通讯的最大优势所在。

• 通讯的范围： 明确哪些任务之间需要通讯在设计并行代码的过程中是非常关键的。下面两种通讯范围既可以被设计为同步的，也可以被设计为异步的：1）点对点通讯： 涉及到两个任务，其中一个扮演消息发送者/生产者的角色，另外一个扮演消息接受者/消费者的角色。2）广播通讯： 涉及到多于两个任务之间的数据共享。这些任务通常处于一个组或者集合中。

• 通讯的效率： 通常编程者具有影响通讯性能的选择，这里列举其中一些：1）对于一个给定的模型，究竟应该采用哪一种实现？例如对于消息传递模型而言，一种MPI的实现可能在某个给定的硬件下比其它实现要快。2）什么采用什么类型的通讯操作？正如前面所提到的，异步通讯操作往往可以提高程序的整体性能。3）网络结构（network fabric）：某些平台可能会提供多于一个的网络结构。那么究竟哪一个最好？

• 死锁与活锁请见：传送门

2.2 同步 (Synchronization)

管理工作的顺序和执行它的任务是大多数并行程序设计的关键，它也可能是提升程序性能的关键，通常也需要对某些程序进行“串行化”。

• 屏障： 1）这通常意味着会涉及到所有任务；2）每个任务都执行自身的工作，直到它遇到屏障，然后它们就停止，或者“阻塞”；3）当最后一个任务到达屏障时，所有任务得以同步；4）接下来可能发生的事情就有所变化了。通常会执行一段串行代码，或者所有的任务在这里都结束了。

• 锁/信号量： 1）可以涉及任意多个任务；2）通常用于对全局数据或者某段代码的存取串行化（保护），在任一时刻，只有一个任务可以使用锁/信号量；3）第一个任务会获得一个锁，然后该任务就可以安全地对该保护数据进行存取；4）其它任务可以尝试去获得锁，但必须等到当前拥有该锁的任务释放锁才行；5）可以是阻塞的也可以是非阻塞的。

• 同步通讯操作： 1）仅仅涉及到执行数据通讯操作的任务；2）当一个任务执行数据通讯操作时，通常需要在参与通讯的任务之间建立某种协调机制。例如，在一个任务发送消息时，它必须收到接受任务的确认，以明确当前是可以发送消息的；3）在消息通讯一节中也已经说明。

2.3粒度 (Granularity)

• 计算通讯比 (computation / Communication Ratio)： 在并行计算中，粒度是对计算与通讯的比例的定性度量。计算周期通常通过同步时间与通讯周期分离。

• 细粒度并行化 (Fine-grain Parallelism)： 1）在通讯事件之外进行相对较少的计算工作；2）计算通讯率较低；3）方便负载均衡；4）意味着较高的通讯开销以及较少的性能提升机会；5）如果粒度过细，任务之间的通讯和同步的开销可能需要比计算更长的时间。

• 粗粒度并行化 (Coarse-grain Parallelism)： 1）在通讯/同步事件之外需要较大量的计算工作；2）较高的计算/通讯比；3）意味着较大的性能提升机会；4）难以进行较好的负载均衡。

• 更多详细内容原文请见：并行计算术语简介

并行计算实现矩阵乘法代码分析

3.1 矩阵并行计算程序(原文章点击此)

（1） 将A和C按行分为np块，将B按列分为np块（B可以按列存储）；
（2） 进程号为 id 的进程读取 A 和 B 的第id个分块；
（3） 循环np次：

• 各个进程用各自的A、B分块求解C的分块；
• 轮换B的分块（例如：id 号进程发送自己当前的B的分块到 id+1号进程）

#include<iostream>//输入输出是由iostream库提供
#include<mpi.h>//mpi接口库
#include<math.h>//数学计算公式的库，包括取绝对值，取证取余和三角函数等的计算。
#include<stdlib.h>//stdlib.h声明的库函数可分为六类：类型转换、伪随机数、动态内存分配与回收管理、进程控制、搜索及排序、简单数学。stdlib.h中定义了物种类型：一些宏和通用工具函数。是C的函数库，对应的C++函数库为cstdlib，C++兼容C，而cstdlib包含stdlib.h所有功能而且以C++风格编写。


void initMatrixWithRV(float *A, int rows, int cols);//RV:Random Variable表示随机变量 init：初始化
void copyMatrix(float *A, float *A_copy, int rows, int cols);//singlet hread:单线程  矩阵C=A*B
// A: m*p, B: p*n  ！！！ note that B is stored by column first（注意B是列优先储存）
void matMultiplyWithTransposedB(float *A, float *B, float *matResult, int m, int n, int p);

int main(int argc, char** argv)
{
　 int m = atoi(argv[1]);//atoi()函数  atoi(): int atoi(const char *str );   功能：把字符串转换成整型数。
　 int n = atoi(argv[2]);
　 int p = atoi(argv[3]);

   float *A, *B, *C;//A B C代表三个矩阵，初始化三个矩阵的地址
   float *bA, *bB_send, *bB_recv, *bC, *bC_send;//初始化各个进程需要的数组的地址
　 int myrank, numprocs;//储存进程号和进程数

    MPI_Status status;//是MPI里的指针

    MPI_Init(&argc, &argv);  // 并行开始
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs); //获得进程数
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank); //获得进程编号
 　　
 　　//行数和列数除进程数向下取整获得的每个分块的行数和列数
　  int bm = m / numprocs;
    int bn = n / numprocs;

    //各个进程需要储存各自处理的矩阵的部分，要存储就要申请空间
    bA = new float[bm * p];
    bB_send = new float[bn * p];
    bB_recv = new float[bn * p];
    bC = new float[bm * bn];
    bC_send = new float[bm * n];

    //主进程存储全部数据
    if(myrank == 0){
        A = new float[m * p];
        B = new float[n * p];
        C = new float[m * n];

        //通过随机生成数的方式来生成需要计算的目标矩阵
        initMatrixWithRV(A, m, p);
        initMatrixWithRV(B, n, p);
    }
    
    MPI_Barrier(MPI_COMM_WORLD);//阻塞，让所有进程同步。只有各个进程申请空间，主进程生成目标矩阵之后才能进行分配
    MPI_Scatter(A, bm * p, MPI_FLOAT, bA, bm * p, MPI_FLOAT, 0, MPI_COMM_WORLD);//将A分配给各个进程（可能会有剩余，即A整除numprocs可能有余数）
    MPI_Scatter(B, bn * p, MPI_FLOAT, bB_recv, bn * p, MPI_FLOAT, 0, MPI_COMM_WORLD);//将B分配给各个进程（可能会有剩余，即B整除numprocs可能有余数）
    //这里sendto取余的原因是最后一位进程需要和第一位进程通信，所以取余可以得到第一位的进程号  recvFrom要加numprocs的原因是防止进程号为负数(myrank是从零开始的，myrank一般是主进程，主进程也跟其他进程一样做运算)
    int sendTo = (myrank + 1) % numprocs;//每个进程得到目标进程的进程号
    int recvFrom = (myrank - 1 + numprocs) % numprocs;//每个进程获得源数据（发送数据）的进程号
    
    int circle = 0;  //分块循环的计数器
    do{
        matMultiplyWithTransposedB(bA, bB_recv, bC, bm, bn, p);
        int blocks_col = (myrank - circle + numprocs) % numprocs;//确定bc在bC_send 中的位置
        for(int i=0; i<bm; i++){
            for(int j=0; j<bn; j++){
                bC_send[i*n + blocks_col*bn + j] = bC[i*bn + j];//bc在bC_send只用考虑将列号变换即可
            }
        } 

        if(myrank % 2 == 0){
            copyMatrix(bB_recv, bB_send, bn, p);
            MPI_Ssend(bB_send, bn*p, MPI_FLOAT, sendTo, circle, MPI_COMM_WORLD);//使用消息标签的意义在于,避免两个进程中传递多个消息时,进程发送消息先后到达出现的冲突.可将本次发送的消息与同一进程向同一目的进程发送的其他消息区分开来.
            MPI_Recv(bB_recv, bn*p, MPI_FLOAT, recvFrom, circle, MPI_COMM_WORLD, &status);
        }else{
            MPI_Recv(bB_recv, bn*p, MPI_FLOAT, recvFrom, circle, MPI_COMM_WORLD, &status); 
            MPI_Ssend(bB_send, bn*p, MPI_FLOAT, sendTo, circle, MPI_COMM_WORLD);
            copyMatrix(bB_recv, bB_send, bn, p);    
        }

        circle++;
    }while(circle < numprocs);

　  MPI_Barrier(MPI_COMM_WORLD);//同步一下，确保每个进程都计算完成
    MPI_Gather(bC_send, bm * n, MPI_FLOAT, C, bm * n, MPI_FLOAT, 0, MPI_COMM_WORLD);//收集数据
   //进行剩余计算
    if(myrank == 0){
        int remainAStartId = bm * numprocs;
        int remainBStartId = bn * numprocs;
   //计算C矩阵中剩余的从remainAStartId行到m行
        for(int i=remainAStartId; i<m; i++){
            for(int j=0; j<n; j++){
                float temp=0;
                for(int k=0; k<p; k++){
                    temp += A[i*p + k] * B[j*p +k];
                }
                C[i*p + j] = temp;
            }
        }
     // 计算C矩阵中从第1行到remainAStartId行，remainBStartId列到n列的数据
        for(int i=0; i<remainAStartId; i++){
            for(int j=remainBStartId; j<n; j++){
                float temp = 0;
                for(int k=0; k<p; k++){
                    temp += A[i*p + k] * B[j*p +k];
                }
                C[i*p + j] = temp;
            }
        }
    }
    //释放空间
    delete[] bA;
    delete[] bB_send;
    delete[] bB_recv;
    delete[] bC;
    delete[] bC_send;

    if(myrank == 0){
        delete[] A;
        delete[] B;
        delete[] C;
    }

    MPI_Finalize(); // 并行结束

    return 0;
}

void initMatrixWithRV(float *A, int rows, int cols)
{
    srand((unsigned)time(NULL));
    for(int i = 0; i < rows*cols; i++){
        A[i] = (float)rand() / RAND_MAX;//RAND_MAX:字符常量 一般数值为32767
    }
}
void copyMatrix(float *A, float *A_copy, int rows, int cols)
{
    for(int i=0; i<rows*cols; i++){
        A_copy[i] = A[i];
    }
}

void matMultiplyWithTransposedB(float *A, float *B, float *matResult, int m, int p, int n)
{
    for(int i=0; i<m; i++){
        for(int j=0; j<n; j++){
            float temp = 0;
            for(int k=0; k<p; k++){
                temp += A[i*p+k] * B[j*p+k];
            }
            matResult[i*n+j] = temp;
        }
    }
}

编译（C++）：

mpicxx   [文件名.后缀] -o  [文件地址]    

（C）：

mpicc [文件名.后缀] -o [文件地址]

运行

mpiexec -np [开启进程数(整数)]  [文件地址]  [程序参数列表]

这里最需要注意的地方就是B的轮换。 有两点需要注意：

(1) 防阻塞机制。这里采用奇偶原则：偶数号进程先发送，再接收；奇数号进程则相反。这样可以避免所有进程同时发送造成死锁的情况；

(2) 数据备份。发送和接收的信息存储在不同的矩阵中，这样保证原来的信息不会被覆盖。
这种方法的优点是显而易见的。对于足够牛的服务器/计算机集群，开启成百上千个进程来并行完全不是问题。