examples/smpi/MM/Summa.c

   1 /*!
   2  * Classical Block Matrix Multiplication example
   3  *
   4  */
   5
   6 #include "Matrix_init.h"
   7 #include "Summa.h"
   8 #include "xbt/log.h"
   9 #include <stdio.h>
  10
  11  XBT_LOG_NEW_DEFAULT_CATEGORY(MM_Summa,
  12                              "Messages specific for this msg example");
  13
  14 double Summa(double *a, double *b, double *c,
  15              size_t lda, size_t ldb, size_t ldc,
  16              size_t m, size_t k_a, size_t k_b, size_t n,
  17              size_t Block_size, size_t start, size_t end,
  18              size_t row, size_t col, size_t size_row, size_t size_col,
  19              double *a_local, double *b_local,
  20              MPI_Datatype Block_a, MPI_Datatype Block_a_local,
  21              MPI_Datatype Block_b,
  22              MPI_Comm row_comm, MPI_Comm col_comm, int subs)
  23 {
  24   double *B_a     , *B_b     ; //matrix blocks
  25   size_t err;
  26   //double alpha = 1, beta = 1;  //C := alpha * a * b + beta * c
  27   size_t B_proc_row; // Number of bloc(row or col) on one processor
  28 #ifndef CYCLIC
  29   size_t B_proc_col;
  30   B_proc_col =  k_b / Block_size;  // Number of block on one processor
  31 #endif
  32   B_proc_row = k_a / Block_size; // Number of block on one processor
  33
  34   //size_t lda = k_a, ldb = n, ldc = n;
  35   size_t lda_local = lda;
  36   size_t ldb_local = ldb;
  37
  38
  39   double time, computation_time = 0, communication_time = 0;
  40   double start_time, end_time; //time mesure
  41   double start_time_intern, end_time_intern; //time mesure
  42
  43
  44
  45
  46   start_time = MPI_Wtime();
  47
  48   /*-------------Distributed Matrix Multiplication algorithm-----------------*/
  49   size_t iter;
  50   for( iter = start; iter < end; iter++ ){
  51     size_t pivot_row, pivot_col, pos_a, pos_b;
  52 #ifdef CYCLIC
  53     // pivot row on processor layer
  54     pivot_row = (iter % size_col);
  55     pivot_col = (iter % size_row);
  56     //position of the block
  57     if(subs == 1){
  58       pos_a = (size_t)((iter - start) / size_row) * Block_size;
  59       pos_b = (size_t)((iter - start) / size_col) * ldb * Block_size;
  60     }else{
  61       pos_a = (size_t)(iter / size_row) * Block_size;
  62       pos_b = (size_t)(iter / size_col) * ldb * Block_size;
  63     }
  64 #else
  65     // pivot row on processor layer
  66     pivot_row = (size_t)(iter / B_proc_col) % size_col;
  67     pivot_col = (size_t)(iter / B_proc_row) % size_row;
  68     //position of the block
  69     if(subs == 1){
  70       pos_a = ((iter - start) % B_proc_row) * Block_size;
  71       pos_b = ((iter - start) % B_proc_col) * ldb * Block_size;
  72     }else{
  73       pos_a = (iter % B_proc_row) * Block_size;
  74       pos_b = (iter % B_proc_col) * ldb * Block_size;
  75     }
  76 #endif
  77     XBT_DEBUG( "pivot: %zu, iter: %zu, B_proc_col: %zu, "
  78                 "size_col:%zu, size_row: %zu\n",
  79                 pivot_row, iter, B_proc_row,size_col,size_row);
  80 /*    MPI_Barrier(row_comm);*/
  81 /*    MPI_Barrier(col_comm);*/
  82
  83     start_time_intern = MPI_Wtime();
  84     //Broadcast the row
  85     if(size_row > 1){
  86       MPI_Datatype * Block;
  87       if( pivot_col != col ){
  88         B_a = a_local;
  89         lda_local = Block_size;
  90         XBT_DEBUG("recieve B_a %zu,%zu \n",m , Block_size);
  91         Block = &Block_a_local;
  92       }else{
  93         B_a = a + pos_a;
  94         lda_local = lda;
  95         XBT_DEBUG("sent B_a %zu,%zu \n",m , Block_size);
  96         Block = &Block_a;
  97       }
  98       err = MPI_Bcast(B_a, 1, *Block, pivot_col, row_comm);
  99       if (err != MPI_SUCCESS) {
 100         perror("Error Bcast A\n");
 101         return -1;
 102       }
 103     }else{
 104       B_a = a + pos_a;
 105       XBT_DEBUG("position of B_a: %zu \n", pos_a);
 106     }
 107
 108     //Broadcast the col
 109     if(size_col > 1){
 110       if( pivot_row == row ){
 111         B_b = b + pos_b;
 112         XBT_DEBUG("sent B_b Block_size: %zu, pos:%zu \n",
 113                     ldb, pos_b);
 114       }else{
 115         B_b = b_local;
 116         XBT_DEBUG("recieve B_b %zu,%zu \n", Block_size,n);
 117       }
 118       err = MPI_Bcast(B_b, 1, Block_b, pivot_row, col_comm );
 119       if (err != MPI_SUCCESS) {
 120         perror("Error Bcast B\n");
 121         MPI_Finalize();
 122         exit(-1);
 123       }
 124     }else{
 125       B_b = b + pos_b;
 126       XBT_DEBUG("position of B_b: %zu \n", pos_b);
 127     }
 128     end_time_intern = MPI_Wtime();
 129     communication_time += end_time_intern - start_time_intern;
 130
 131 /*    MPI_Barrier(row_comm);*/
 132 /*    MPI_Barrier(col_comm);*/
 133     start_time_intern = MPI_Wtime();
 134     XBT_DEBUG("execute Gemm number: %zu\n", iter);
 135     //We have recieved a line of block and a colomn
 136    //              cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
 137    //               m, n, Block_size, alpha, B_a, lda_local, B_b, ldb_local,
 138    //               beta, c, ldc );
 139     int i, j, k;
 140     for(i = 0; i < m; i++)
 141       for(j = 0; j < n; j++)
 142         for(k = 0; k < Block_size; k++)
 143           c[i*ldc+j] += B_a[i*lda_local+k]*B_b[k*ldb_local+j];
 144
 145     end_time_intern = MPI_Wtime();
 146     computation_time += end_time_intern - start_time_intern;
 147
 148   }
 149   MPI_Barrier(row_comm);
 150   MPI_Barrier(col_comm);
 151
 152   end_time = MPI_Wtime();
 153   time = end_time - start_time ;
 154   printf("communication time: %e seconds, "
 155          "computation time: %e seconds\n",
 156          communication_time, computation_time);
 157
 158
 159   return time;
 160 }