matrixvector.c
上传用户:jsljixie
上传日期:2013-08-15
资源大小:827k
文件大小:9k
源码类别:

并行计算

开发平台:

Visual C++

  1. /*
  2.  * 程序功能:该程序并行计算矩阵和向量的乘积,并将结果向量输出,采用带状划分的算法
  3.  * 变量说明:matrix[width,width]   存储从文件中读取的用于计算的矩阵的元素,矩阵元素为单精度浮点数
  4.  *           vec[width]    存储从文件中读取的用于计算的向量的元素,向量元素为单精度浮点数
  5.  *           width         存储从文件中读取的未知矩阵的宽度,和向量的长度,用来动态给
  6.  *                                  矩阵头指针matrix,和向量头指针vec分配空间,      整型变量
  7.  *           col            默认的值应该为width+1,如果不是报错f           整型变量
  8.  *           starttime,finishprepare,finishcompute,finishall
  9.  *                         上述四个全局变量用来记录四个时间点                        双精度浮点数
  10.  *           status         MPI接受消息时用到的变量      MPI_statuc 变量
  11.  *           rom_per_pro    记录每一个进程实际需要计算的行数               整形变量
  12.  *           fin            输入文件指针
  13.  *           i,j            控制循环,并用作消息传递时的tag       整形变量
  14.  *           my_rank        存储当前进程的标识                    整型变量
  15.  *           group_size     存储当前总的进程的总数                整型变量
  16.  *           vecin[width]   在每一个进程中用来暂存vec的变量       元素为单精度浮点数
  17.  *           ma[rom_per_pro,width]
  18.  *                          每一个进程中用来暂存被分配的那部分的矩阵的元素  元素为单精度浮点数
  19.  *           result[rom_per_pro]
  20.  *                          每一个进程中用来记录计算结果的变量    元素为单精度浮点数
  21.  *           sum            计算过程中用到的中间变量              单精度浮点数
  22.  *           p              处理器的总数                           整型变量
  23.  *           realp          实际用到的处理器的数目                 整型变量
  24.  */
  25. #include "stdio.h"
  26. #include "stdlib.h"
  27. #include "mpi.h"
  28. #include "math.h"
  29. #define matrix(x,y) matrix[x*width+y]
  30. #define vec(x) vec[x]
  31. #define vecin(x) vecin[x]
  32. #define result(x) result[x]
  33. #define ma(x,y) ma[x*width+y]
  34. #define intsize sizeof(int)
  35. #define floatsize sizeof(float)
  36. int realp,p;
  37. int width,col,row_per_pro;
  38. float *matrix;
  39. float *vec;
  40. double starttime,finishprepare,finishcompute,finishall;
  41. MPI_Status status;
  42. FILE *fin;
  43. /*
  44.  函数名: free_all
  45.  功能:释放每个进程中动态分配的三个变量空间
  46.  输入:三个双精度浮点数指针变量ma,result,vecin
  47.  输出:返回为空
  48. */
  49. void free_all(float *ma,float *result,float *vecin)
  50. {
  51.     free(ma);
  52.     free(result);
  53.     free(vecin);
  54. }
  55. /*
  56.  函数名: main
  57.  功能:从文件中读入一个矩阵和一个向量,并行他们的乘积,在屏幕上输出结果变量
  58.  输入:argc为命令行参数个数
  59.        argv为每个命令行参数组成的字符串数组。
  60.  输出:返回0代表程序正常结束
  61. */
  62. int main(int argc, char **argv)
  63. {
  64.     int i,j,my_rank,group_size;
  65.     float sum;
  66.     float *ma;
  67.     float *result;
  68.     float *vecin;
  69.     MPI_Init(&argc,&argv);
  70.     MPI_Comm_size(MPI_COMM_WORLD,&group_size);
  71.     MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
  72.     p=group_size;
  73.     /*
  74.      由进程0完成matrix,vec的动态空间分配
  75.      同时将向量和矩阵读入到vec和matrix的空间中
  76.     */
  77.     if(my_rank==0)
  78.     {
  79.         /*
  80.          记下初始运行的时间点
  81.          打开文件dataIn.txt并读出width和col值
  82.          行列数不对,报错并中断返回
  83.         */
  84.         starttime=MPI_Wtime();
  85.         fin=fopen("dataIn.txt","r");
  86.         fscanf(fin,"%d %d", &width, &col);
  87.         if (width!=(col-1))
  88.         {
  89.             printf("The input is wrong,check the input matrix width and the total columnsnnnn");
  90.             exit(1);
  91.         }
  92.         /*
  93.          在进程0中为matrix和vec分配空间从dataIn.txt中读出相应的值
  94.          matrix 分配了width*width大小的浮点数空间
  95.          vec分配了width长的浮点数空间
  96.         */
  97.         matrix=(float *)malloc(floatsize*width*width);
  98.         vec=(float *)malloc(floatsize*width);
  99.         for(i = 0; i < width; i++)
  100.         {
  101.             for(j = 0; j < width; j++)
  102.             {
  103.                 fscanf(fin,"%f", matrix+i*width+j);
  104.             }
  105.             fscanf(fin,"%f", vec+i);
  106.         }
  107.         /* 将从文件dataIn.txt中读入的结果显示在屏幕上 */
  108.         fclose(fin);
  109.         printf("The matrix and the vector has been read from file "dataIn.txt" is: %d widthn ",width);
  110.         printf("The element of the matrix is as follows:n");
  111.         for(i=0;i<width;i++)
  112.         {
  113.             printf("row %dt",i);
  114.             for(j=0;j<width;j++)
  115.                 printf("%ft",matrix(i,j));
  116.             printf("n");
  117.         }
  118.         printf("the element of the vector is as follows:n");
  119.         for(i=0;i<width;i++)
  120.         {
  121.             printf("row %d is %fn",i,vec(i));
  122.         }
  123.     }
  124.     /* 进程0已经得到width值,所以它将这个结果广播给其他进程 */
  125.     MPI_Bcast(&width,1,MPI_INT,0,MPI_COMM_WORLD);
  126.     /*
  127.      计算每一个进程要处理的行数,注意不能整除时要使得row_per_pro的值加1
  128.      计算实际需要的最多的进程数,对于用不上的进程结束其执行
  129.      将实际用到的进程数和每个进程需要处理的行数输出到屏幕上
  130.     */
  131.     row_per_pro=width/p;
  132.     if (width%p!=0)
  133.         row_per_pro++;
  134.     realp=0;
  135.     while(row_per_pro*realp<width)
  136.         realp++;
  137.     if(my_rank==0)
  138.         printf("now there is %d processors are working on it and %d rows per processorsn",realp,row_per_pro);
  139.     /* 现在为每一个进程中的三个浮点数指针变量分配空间 */
  140.     vecin=(float *)malloc(floatsize*width);
  141.     ma=(float *)malloc(floatsize*row_per_pro*width);
  142.     result=(float *)malloc(floatsize*row_per_pro);
  143.     /* 将进程0中的vec向量的值广播到各个进程的vecin变量中保存 */
  144.     if (my_rank==0)
  145.     {
  146.         for(i=0;i<width;i++)
  147.             vecin(i)=vec(i);
  148.         free(vec);
  149.     }
  150.     MPI_Bcast(vecin,width,MPI_FLOAT,0,MPI_COMM_WORLD);
  151.     /* 将进程0中的matrix的各个元素按各个进程的分工进行传播 */
  152.     if (my_rank==0)
  153.     {
  154.         for(i=0;i<row_per_pro;i++)
  155.             for(j=0;j<width;j++)
  156.                 ma(i,j)=matrix(i,j);
  157.     }
  158.     if (my_rank==0)
  159.     {
  160.         if(width%realp==0)
  161.         {
  162.             for(i=1;i<realp;i++)
  163.                 MPI_Send(&(matrix(row_per_pro*i,0)),row_per_pro*width,MPI_FLOAT,i,i,MPI_COMM_WORLD);
  164.         }
  165.         else
  166.         {
  167.             for(i=1;i<(realp-1);i++)
  168.                 MPI_Send(&(matrix(row_per_pro*i,0)),row_per_pro*width,MPI_FLOAT,i,i,MPI_COMM_WORLD);
  169.             MPI_Send(&(matrix(row_per_pro*(realp-1),0)),
  170.                      (width-row_per_pro*(realp-1))*width,
  171.                      MPI_FLOAT,(realp-1),(realp-1),MPI_COMM_WORLD);
  172.         }
  173.         free(matrix);
  174.     }
  175.     else
  176.     {
  177.         if (my_rank<realp)
  178.         {
  179.             if(width%realp==0)
  180.                 MPI_Recv(ma,row_per_pro*width,MPI_FLOAT,0,my_rank,MPI_COMM_WORLD,&status);
  181.             else
  182.             {
  183.                 if(my_rank!=realp-1)
  184.                     MPI_Recv(ma,row_per_pro*width,MPI_FLOAT,0,my_rank,MPI_COMM_WORLD,&status);
  185.                 else
  186.                     MPI_Recv(ma,(width-row_per_pro*(realp-1))*width,MPI_FLOAT,0,realp-1,MPI_COMM_WORLD,&status);
  187.             }
  188.         }
  189.     }
  190.     /* 所有的数据都已经传送完毕,同步各个进程并,记录现在的时间点 */
  191.     MPI_Barrier(MPI_COMM_WORLD);
  192.     finishprepare=MPI_Wtime();
  193.     /* 各个进程作自己那一部分的计算工作 */
  194.     if (my_rank==0)
  195.     {
  196.         for(i=0;i<row_per_pro;i++)
  197.         {
  198.             sum=0.0;
  199.             for(j=0;j<width;j++)
  200.                 sum=sum+ma(i,j)*vecin(j);
  201.             result(i)=sum;
  202.         }
  203.     }
  204.     else
  205.     {
  206.         if(my_rank<(realp-1))
  207.         {
  208.             for(i=0;i<row_per_pro;i++)
  209.             {
  210.                 sum=0.0;
  211.                 for(j=0;j<width;j++)
  212.                     sum=sum+ma(i,j)*vecin(j);
  213.                 result(i)=sum;
  214.             }
  215.         }
  216.         else
  217.         {
  218.             if (my_rank<realp)
  219.             {
  220.                 for(i=0;i<(width-row_per_pro*(realp-1));i++)
  221.                 {
  222.                     sum=0.0;
  223.                     for(j=0;j<width;j++)
  224.                         sum=sum+ma(i,j)*vecin(j);
  225.                     result(i)=sum;
  226.                 }
  227.             }
  228.         }
  229.     }
  230.     /* 每一个进程都计算完毕,记下此时的时间点 */
  231.     MPI_Barrier(MPI_COMM_WORLD);
  232.     finishcompute=MPI_Wtime();
  233.     /* 将结果传回进程0 */
  234.     if (my_rank==0)
  235.     {
  236.         for(i=0;i<row_per_pro;i++)
  237.             vecin(i)=result(i);
  238.         for(i=1;i<(realp-1);i++)
  239.             MPI_Recv(&(vecin(i*row_per_pro)),row_per_pro,MPI_FLOAT,i,i,MPI_COMM_WORLD,&status);
  240.         if(realp!=1)
  241.         {
  242.             MPI_Recv(&(vecin(row_per_pro*(realp-1))),
  243.                      (width-row_per_pro*(realp-1)),
  244.                      MPI_FLOAT,(realp-1),(realp-1),
  245.                      MPI_COMM_WORLD,&status);
  246.         }
  247.     }
  248.     else
  249.     {
  250.         if (my_rank<realp)
  251.         {
  252.             if (my_rank!=(realp-1))
  253.                 MPI_Send(result,row_per_pro,MPI_FLOAT,0,my_rank,MPI_COMM_WORLD);
  254.             else
  255.                 MPI_Send(result,(width-row_per_pro*(realp-1)),MPI_FLOAT,0,(realp-1),MPI_COMM_WORLD);
  256.         }
  257.     }
  258.     /* 将进程0内的vecin内保存的结果输出 */
  259.     if (my_rank==0)
  260.     {
  261.         printf("the result is as follows:n");
  262.         for(i=0;i<width;i++)
  263.             printf("row%d %fn",i,vecin(i));
  264.     }
  265.     MPI_Barrier(MPI_COMM_WORLD);
  266.     finishall=MPI_Wtime();
  267.     /* 最后将几个时间输出 */
  268.     if (my_rank==0)
  269.     {
  270.         printf("n");
  271.         printf("Whole running time    = %f secondsn",finishall-starttime);
  272.         printf("prepare time  = %f secondsn",finishprepare-starttime);
  273.         printf("Parallel compute time = %f secondsn",finishcompute-finishprepare);
  274.     }
  275.     MPI_Barrier(MPI_COMM_WORLD);
  276.     MPI_Finalize();
  277.     free_all(ma,result,vecin);
  278.     return (0);
  279. }