connon.txt

connon并行矩阵运算

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


float **A, **B, **C;              
float *a, *b, *c, *tmp_a, *tmp_b;
int dg, dl, dl2,p, sp;           
int my_rank, my_row, my_col; 
MPI_Status status;



int get_index(int row, int col, int sp)
{
   return ((row+sp)%sp)*sp + (col+sp)%sp;
}



/* 函数名：scatter_A_B
 * 功能：rank为0的处理器向其他处理器发送A、B矩阵的相关块
 */
void scatter_A_B()
{
   int i,j,k,l;
   int p_imin,p_imax,p_jmin,p_jmax;

   for(k=0; k<p; k++)

   {
	  
	  p_jmin = (k % sp    ) * dl;
  	  p_jmax = (k % sp + 1) * dl-1;
	  p_imin = (k - (k % sp))/sp * dl;
	  p_imax = ((k - (k % sp))/sp +1) *dl -1;
      l = 0;

     
      for(i=p_imin; i<=p_imax; i++)
      {
      	  for(j=p_jmin; j<=p_jmax; j++)
      	  {
              tmp_a[l] = A[i][j];
	      tmp_b[l] = B[i][j];
	      l++;
          }
      }

     
      if(k==0)
      {
         memcpy(a, tmp_a, dl2 * sizeof(float));
	 memcpy(b, tmp_b, dl2 * sizeof(float));
      } else   
      {
          MPI_Send(tmp_a, dl2, MPI_FLOAT, k, 1, MPI_COMM_WORLD);
	  MPI_Send(tmp_b, dl2, MPI_FLOAT, k, 2, MPI_COMM_WORLD);
      }
   }
}

/*
 *函数名:init_alignment
 *功能:矩阵A和B初始对准
 */
void init_alignment()
{
   /*将A中坐标为(i,j)的分块A(i,j)向左循环移动i步*/
   MPI_Sendrecv(a, dl2, MPI_FLOAT, get_index(my_row,my_col-my_row,sp), 1,
            tmp_a, dl2, MPI_FLOAT, get_index(my_row,my_col+my_row,sp), 1, MPI_COMM_WORLD, &status);
   memcpy(a, tmp_a, dl2 * sizeof(float) );

   /*将B中坐标为(i,j)的分块B(i,j)向上循环移动j步*/
   MPI_Sendrecv(b, dl2, MPI_FLOAT, get_index(my_row-my_col,my_col,sp), 1,
            tmp_b, dl2, MPI_FLOAT, get_index(my_row+my_col,my_col,sp), 1, MPI_COMM_WORLD, &status);
   memcpy(b, tmp_b, dl2 * sizeof(float) );
}

/*
 *函数名：main_shift
 *功能：分块矩阵左移和上移，并计算分块c
 */
void main_shift()
{
   int i,j,k,l;

   for(l=0; l<sp; l++)
   {
     /*矩阵块相乘，c+=a*b */
     for(i=0; i<dl; i++)
       for(j=0; j<dl; j++)
         for(k=0; k<dl; k++)
           c[i*dl+j] += a[i*dl+k]*b[k*dl+j];

 

      /* 将分块b上移1位 */
      MPI_Send(b , dl2, MPI_FLOAT, get_index(my_row-1, my_col, sp), 1, MPI_COMM_WORLD);
      MPI_Recv(b , dl2, MPI_FLOAT, get_index(my_row+1, my_col, sp), 1, MPI_COMM_WORLD, &status);
   }
}

/*
 *函数名：collect_c
 *功能：rank为0的处理器从其余处理器收集分块矩阵c
 */
void collect_C()
{
   int i,j,i2,j2,k;
   int p_imin,p_imax,p_jmin,p_jmax; 

   
   for (i=0;i<dl;i++)
	 for(j=0;j<dl;j++)
	   C[i][j]=c[i*dl+j];

   for (k=1;k<p;k++)
   {
      
       MPI_Recv(c, dl2, MPI_FLOAT, k, 1, MPI_COMM_WORLD, &status);

       p_jmin = (k % sp    ) *dl;
       p_jmax = (k % sp + 1) *dl-1;
       p_imin =  (k - (k % sp))/sp     *dl;
       p_imax = ((k - (k % sp))/sp +1) *dl -1;

       i2=0;
      
       for(i=p_imin; i<=p_imax; i++)
       {
           j2=0;
           for(j=p_jmin; j<=p_jmax; j++)
           {
               C[i][j]=c[i2*dl+j2];
               j2++;
           }
           i2++;
       }
   }
}

/*函数名：print
 *功能：打印矩阵
 *输入：指向矩阵指针的指针，字符串
 */
void print(float **m,char *str)
{
   int i,j;
   printf("%s",str);
   /*打印矩阵m*/
   for(i=0;i<dg;i++)
   {
       for(j=0;j<dg;j++)
           printf("%15.0f    ",m[i][j]);
       printf("\n");
   }
   printf("\n");
}

/*
 *函数名：main
 *功能：主过程，Cannon算法，矩阵相乘
 *输入：argc为命令行参数个数，argv为每个命令行参数组成的字符串数组
 */
int main(int argc, char *argv[])
{
   int i;

   MPI_Init(&argc, &argv);                  
   MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); 

   sp = sqrt(p);

   /* 确保处理器个数是完全平方数，否则打印错误信息，程序退出 */
   if (sp*sp != p)
   {
      if (my_rank == 0)
	  printf("Number of processors is not a quadratic number!\n");
      MPI_Finalize();
      exit(1);
   }

   if (argc != 2)
   {
      if (my_rank == 0)
          printf("usage: mpirun -np ProcNum cannon MatrixDimension\n");
      MPI_Finalize();
      exit(1);
   }

   dg  = atoi(argv[1]);    /* 总矩阵维数 */
   dl  = dg / sp;          /* 计算分块矩阵维数 */
   dl2 = dl * dl;

   /* 计算处理器在逻辑阵列中的坐标 */
   my_col =  my_rank % sp ;
   my_row = (my_rank-my_col) / sp ;

   /* 为a、b、c分配空间 */
   a = (float *)malloc( dl2 * sizeof(float) );
   b = (float *)malloc( dl2 * sizeof(float) );
   c = (float *)malloc( dl2 * sizeof(float) );

   /* 初始化c */
   for(i=0; i<dl2 ; i++)
     c[i] = 0.0;

   /* 为tmp_a、tmp_b分配空间 */
   tmp_a = (float *)malloc( dl2 * sizeof(float) );
   tmp_b = (float *)malloc( dl2 * sizeof(float) );

   if (my_rank == 0)
   {
      /* rank为0的处理器为A、B、C分配空间 */
      A = (float **)malloc( dg * sizeof(float*) );
      B = (float **)malloc( dg * sizeof(float*) );
      C = (float **)malloc( dg * sizeof(float*) );

      for(i=0; i<dg; i++)
      {
         A[i] = (float *)malloc( dg * sizeof(float) );
         B[i] = (float *)malloc( dg * sizeof(float) );
         C[i] = (float *)malloc( dg * sizeof(float) );
      }
      random_A_B();    
      scatter_A_B();   
   } else               
   {
       MPI_Recv(a, dl2, MPI_FLOAT, 0 , 1, MPI_COMM_WORLD, &status);
       MPI_Recv(b, dl2, MPI_FLOAT, 0 , 2, MPI_COMM_WORLD, &status);
   }

   
   {
      MPI_Send(c,dl2,MPI_FLOAT,0,1,MPI_COMM_WORLD); 
   }

   MPI_Barrier(MPI_COMM_WORLD);        /* 同步所有处理器 */
   MPI_Finalize();                     /* 结束MPI计算 */

   return 0;
}