【MPI学习笔记】2：并行化方阵和向量的乘积(按行分配)

按照老师所说，可以把矩阵的每一行都存其列号(0~N-1)，然后列向量全部设置为1，这样得到的结果列向量一定每一位的值都应当是(N-1+0)*N/2，可以用这种方式检查程序写的对不对。

每个进程一次读完自己的任务

简述

假设机器是共享内存的(如果是分布式内存的，那16台这么大内存的机器能处理的规模是现在的16倍)，那么必须为所有进程考虑这点可用的空间，拿N=40000试了一下(死循环扔进后台然后查看内存使用情况)：

因为最主要的就是存这个大方阵，它所开的空间是N*N*sizeof(double)就等于40000*40000*8个字节，差不多是12.8个GB，如上图。

按照这种计算方式，用阶数是50000的方阵时所需要的空间差不多是20个G，一定不够用。

不过这种死循环式的测试手段似乎有问题，没有考虑操作系统本身。操作系统依据程序局部性原理，把不断命中的while(1){}这部分语句已经保持在了内存甚至cache里，然而这部分语句完全没有用到我要测试的那些malloc出来的内存里的数据，所以没有使用的那部分数据或许被暂时换到外存，所以这种测试得到的结果是这样：

实际上在这个程序里，不是真的能跑这么大的，如果放到前台来执行，跑不出结果来：

只好Ctrl+Z停止然后kill掉，在kill之前可以看一下主存和交换分区都满了，我猜测是因为不停的缺页替换所以始终没法得到结果。

把自己能处理的jobs都kill掉，然后看一下到底有多少内存可用：

差不多19980000MB=1.998x10^10B，4万多阶的double方阵就是极限了。

程序

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

//方阵的维度
#define N 40000

int main()
{
    int *vec=NULL;//列向量
    double *mat=NULL;//自己进程的那部分矩阵
    int my_rank;//自己进程的进程号
    int comm_sz;//总的进程数目
    int my_row;//本进程处理的行数
    int i,j;//游标
    double *result=NULL;//用来存本进程计算出的结果向量
    double *all_rst=NULL;//只给0号进程存总的结果

    /**********************/
    MPI_Init(NULL,NULL);
    MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
    MPI_Comm_size(MPI_COMM_WORLD,&comm_sz);
    //本进程处理的行数就是总阶数/进程数
    my_row=N/comm_sz;

    //每个进程都申请空间
    vec=malloc(N*sizeof(int));
    mat=malloc(N*my_row*sizeof(double));
    result=malloc(my_row*sizeof(double));

    //结果向量的分量初始值置0
    for(i=0;i<my_row;i++)
        result[i]=0;

    //矩阵的值的设定
    for(i=0;i<my_row;i++)
        for(j=0;j<N;j++)
            mat[i*N+j]=j;
    //向量的值的设定
    for(j=0;j<N;j++)
        vec[j]=1;

    //进程同步,让所有进程的资源都准备好
    MPI_Barrier(MPI_COMM_WORLD);

    //复杂度O(N*my_row)=O(N*N/comm_sz)
    for(i=0;i<my_row;i++)
    {
        for(j=0;j<N;j++)
        {
            //计算并加进来
            result[i]+=mat[i*N+j]*vec[j];
        }
    }

    //while(0){}

    //进程同步,让所有进程的资源都准备好
    MPI_Barrier(MPI_COMM_WORLD);

    //聚集给0号进程
    if(my_rank==0)
    {
        //先开辟存储空间
        all_rst=malloc(N*sizeof(double));

        //聚集
        MPI_Gather
            (
            result, /*发送内容的地址*/
            my_row, /*发送的长度*/
            MPI_DOUBLE, /*发送的数据类型*/
            all_rst, /*接收内容的地址*/
            my_row, /*接收的长度*/
            MPI_DOUBLE, /*接收的数据类型*/
            0, /*接收至哪个进程*/
            MPI_COMM_WORLD /*通信域*/
            );
    }
    else
    {
        //聚集
        MPI_Gather
            (
            result, /*发送内容的地址*/
            my_row, /*发送的长度*/
            MPI_DOUBLE, /*发送的数据类型*/
            all_rst, /*接收内容的地址*/
            my_row, /*接收的长度*/
            MPI_DOUBLE, /*接收的数据类型*/
            0, /*接收至哪个进程*/
            MPI_COMM_WORLD /*通信域*/
            );
    }

    //进程同步,让所有进程的资源都准备好
    MPI_Barrier(MPI_COMM_WORLD);

    //0号进程负责输出
    if(my_rank==0)
    {
        printf("The Result is:\n");
        for(i=0;i<N;i+=1)
            printf("%f\n",all_rst[i]);
    }

    MPI_Finalize();
    /**********************/

    free(vec);
    free(mat);
    free(result);
    free(all_rst);

    return 0;
}

运行结果

方阵每行是从0到4万-1，用计算器计算一下4万*(4万-1)/2，结果正确。

每个进程按行读取自己的任务

简述

受群里ll同学的方法启发，他的方式是把每个进程自己的小矩阵，再分配成多行一组的行组，每次计算出这些行组对应的子向量后，把行组free掉重新申请，这样只需要行组这么多的空间，就可以轮流读入和计算更多行的子矩阵。

我让这种方式处在极限状态下，即每行本身是前面所说的一个行组。
申请好这一行的空间后，进程每次读入一行覆盖上一次读入的这行，即反复使用这一行的空间，计算并放入自己负责的子向量中去。
因为每次只是覆盖上一次的值，一行始终是这么多数，不需要重复free和*alloc申请空间。

如果每个进程本身在单独的一台机器上，这种方式理论上能计算的方阵的阶的极限是：用进程所在的机器几乎全部的内存表示到一行，所能表示的向量阶数。如果每台机器都是16G的内存，前面计算了能开的大矩阵也就是4万多阶，就按4万阶算，那么理论上能计算约40000*40000=16亿行。

程序

ll同学在群里跑到了96万，我也按96万试一下(实际上这两种方式在这台机器上能跑的阶都远不止96万，但是跑的时间实在太长了)。

在测试中，为了快速获得结果和正当结束程序，可以把对一行值的设定放到循环外面(反正测试的是都一样的值)，并把输出时跳跃的步长从1改成几万，这样很快就会采样输出完然后结束。

/*
15121856 刘知昊
第1次作业(方阵*向量,按行分配)
 */
#include<stdio.h>
#include<mpi.h>
#include<stdlib.h>

//方阵的维度
#define N 960000

int main()
{
    int *vec=NULL;//列向量
    //double *mat=NULL;//自己进程的那部分矩阵
    int my_rank;//自己进程的进程号
    int comm_sz;//总的进程数目
    int my_row;//本进程处理的行数
    int i,j;//通用游标
    double *result=NULL;//用来存本进程计算出的结果向量
    double *all_rst=NULL;//只给0号进程存总的结果
    double *row_now=NULL;//每个进程每次仅申请的一行

    /**********************/
    MPI_Init(NULL,NULL);
    MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
    MPI_Comm_size(MPI_COMM_WORLD,&comm_sz);
    //本进程处理的行数就是总阶数/进程数
    my_row=N/comm_sz;

    //为每个进程都申请空间
    vec=malloc(N*sizeof(int));
    //mat=malloc(N*my_row*sizeof(double));//这是个进程每次申请完
    row_now=malloc(N*sizeof(double));//这是每个进程每次仅申请一行
    result=calloc(sizeof(double),my_row);//为了初始值置0而使用calloc

    //向量的值的设定,使用memset会报错
    for(j=0;j<N;j++)
        vec[j]=1;

    /*
    //本行的值的设定
    //在实际使用时,这里是读入本进程应处理的下一行
    for(j=0;j<N;j++)
        row_now[j]=j;
    */

    //复杂度O(N*my_row)=O(N*N/comm_sz)
    for(i=0;i<my_row;i++)
    {
        //本行的值的设定
        //在实际使用时,这里是读入本进程应处理的下一行
        for(j=0;j<N;j++)
            row_now[j]=j;

        for(j=0;j<N;j++)
        {
            //计算并加入到本进程结果向量的对应位置上
            result[i]+=row_now[j]*vec[j];
        }

        //下次循环直接覆盖上次使用的值
        //因此不需free和重新申请row_now
    }

    //为了机器考虑,确定不再使用的空间立即free掉
    free(row_now);

    //while(0){}//测试用

    //此处的进程同步是必要之举
    MPI_Barrier(MPI_COMM_WORLD);

    //聚集给0号进程
    if(my_rank==0)
    {
        //先开辟存储空间
        all_rst=malloc(N*sizeof(double));

        //聚集
        MPI_Gather
            (
            result, /*发送内容的地址*/
            my_row, /*发送的长度*/
            MPI_DOUBLE, /*发送的数据类型*/
            all_rst, /*接收内容的地址*/
            my_row, /*接收的长度*/
            MPI_DOUBLE, /*接收的数据类型*/
            0, /*接收至哪个进程*/
            MPI_COMM_WORLD /*通信域*/
            );
    }
    else
    {
        //聚集
        MPI_Gather
            (
            result, /*发送内容的地址*/
            my_row, /*发送的长度*/
            MPI_DOUBLE, /*发送的数据类型*/
            all_rst, /*接收内容的地址*/
            my_row, /*接收的长度*/
            MPI_DOUBLE, /*接收的数据类型*/
            0, /*接收至哪个进程*/
            MPI_COMM_WORLD /*通信域*/
            );
    }

    //进程同步,让所有进程的资源都准备好,稳妥之举
    MPI_Barrier(MPI_COMM_WORLD);

    //为了机器考虑,确定不再使用的空间立即free掉
    free(vec);
    free(result);
    //free(mat);


    //0号进程负责输出
    if(my_rank==0)
    {
        printf("The Result is:\n");
        //改变跨度可以采样获取结果,快速结束I/O
        for(i=0;i<N;i+=1)
            printf("%f\n",all_rst[i]);
    }

    MPI_Finalize();
    /**********************/

    //最终,free应无遗漏
    free(all_rst);

    return 0;
}

运行结果

方阵每行是从0到96万-1，用计算器计算一下96万*(96万-1)/2，结果正确。