按照老师所说，可以把矩阵的每一行都存其列号(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，结果正确。
扫一扫
1899
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
