examples/smpi/MM/2.5D_MM.c

   1 /* Copyright (c) 2012-2014. The SimGrid Team.
   2  * All rights reserved.                                                     */
   3
   4 /* This program is free software; you can redistribute it and/or modify it
   5  * under the terms of the license (GNU LGPL) which comes with this package. */
   6
   7 /*!
   8  * 2.5D Block Matrix Multiplication example
   9  *
  10  */
  11
  12 #include "Matrix_init.h"
  13 #include "Summa.h"
  14 #include "2.5D_MM.h"
  15 #include <stdlib.h>
  16 #include <stdio.h>
  17 #include "xbt/log.h"
  18 #define CHECK_25D 1
  19
  20  XBT_LOG_NEW_DEFAULT_CATEGORY(25D_MM,
  21                              "Messages specific for this msg example");
  22
  23 double two_dot_five(
  24                     size_t m, size_t k, size_t n,
  25                     size_t Block_size, size_t group, size_t key,
  26                     size_t size_row, size_t size_col, size_t NB_groups ){
  27   double *a,  *b,  *c;
  28   double *res = NULL;
  29   /* Split the communicator into groups */
  30
  31   /* Find out my identity in the default communicator */
  32   int myrank;
  33   int NB_proc;
  34   int err;
  35   int useless = 0;
  36
  37
  38   double time, communication_time = 0;
  39   double start_time, end_time; //time mesure
  40   double end_time_intern; //time mesure
  41   double start_time_reduce, end_time_reduce; //time mesure
  42
  43   MPI_Comm my_world;
  44
  45   if ( group >= NB_groups ){
  46     XBT_DEBUG(
  47                 "Not enough group NB_groups : %zu my group id : %zu\n",
  48                 NB_groups, group);
  49     MPI_Comm_split(MPI_COMM_WORLD, 0, key, &my_world);
  50     return -1;
  51   }else{
  52     MPI_Comm_split(MPI_COMM_WORLD, 1, key, &my_world);
  53   }
  54
  55   MPI_Comm_size (my_world, &NB_proc);
  56
  57   if ( NB_proc < (int)(size_row*size_col*NB_groups) ){
  58     XBT_INFO(
  59                 "Not enough processors NB_proc : %d required : %zu\n",
  60                 NB_proc, size_row*size_col*NB_groups);
  61     return -1;
  62   }
  63
  64
  65
  66   MPI_Comm group_comm;
  67   MPI_Comm_split(my_world, group, key, &group_comm);
  68
  69   MPI_Comm_rank(group_comm, &myrank);
  70   MPI_Comm_size (group_comm, &NB_proc);
  71   /* for each group start the execution of his */
  72
  73
  74   NB_proc=size_row*size_col;
  75   size_t row = myrank / size_row;
  76   size_t col = myrank % size_row;
  77
  78
  79
  80
  81   /*-------------------------Check some mandatory conditions------------------*/
  82   size_t NB_Block = k / Block_size;
  83   if ( k % Block_size != 0 ){
  84     XBT_INFO(
  85                 "The matrix size has to be proportionnal to the number\
  86                 of blocks: %zu\n", NB_Block);
  87     return -1;
  88   }
  89
  90   if ( size_row > NB_Block || size_col > NB_Block ){
  91     XBT_INFO(
  92                 "Number of blocks is too small compare to the number of"
  93                 " processors (%zu,%zu) in a row or a col (%zu)\n",
  94                 size_col, size_row, NB_Block);
  95     return -1;
  96   }
  97
  98   if ( NB_Block % size_row != 0 || NB_Block % size_col != 0){
  99     XBT_INFO("The number of Block by processor is not an %s",
 100                 "integer\n");
 101     return -1;
 102   }
 103
 104
 105
 106   if(row >= size_col || col >= size_row){
 107     XBT_INFO( "I'm useless bye!!! col: %zu row: %zu, "
 108                 "size_col: %zu , size_row: %zu \n",
 109                 col,row,size_col,size_row);
 110     useless = 1;
 111   }
 112
 113
 114   if(useless == 1){
 115     /*----------------------Prepare the Communication Layer-------------------*/
 116     /* add useless processor on a new color to execute the matrix
 117      * multiplication with the other processors*/
 118
 119     /* Split comm size_to row and column comms */
 120     MPI_Comm row_comm, col_comm, group_line;
 121     MPI_Comm_split(my_world, myrank, MPI_UNDEFINED, &group_line);
 122     /* color by row, rank by column */
 123     MPI_Comm_split(group_comm, size_row, MPI_UNDEFINED, &row_comm);
 124     /* color by column, rank by row */
 125     MPI_Comm_split(group_comm, size_col, MPI_UNDEFINED, &col_comm);
 126     /*------------------------Communication Layer can be used-----------------*/
 127
 128     return 0;
 129   }
 130   XBT_DEBUG( "I'm initialized col: %zu row: %zu, "
 131               "size_col: %zu , size_row: %zu, my rank: %d \n"
 132               ,col,row,size_col,size_row, myrank);
 133
 134
 135
 136   /*------------------------Initialize the matrices---------------------------*/
 137
 138   /* think about a common interface
 139    *  int pdgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
 140    *                 m, n, k, alpha, a, ia, ja, lda, b, ib, jb, ldb,
 141    *                 beta, c, ldc, Comm, rank );
 142    */
 143
 144   /*------------------------Prepare the Communication Layer-------------------*/
 145   /* Split comm size_to row and column comms */
 146   MPI_Comm row_comm, col_comm, group_line;
 147   MPI_Comm_split(my_world, myrank, group, &group_line);
 148   /* color by row, rank by column */
 149   MPI_Comm_split(group_comm, row, col, &row_comm);
 150   /* color by column, rank by row */
 151   MPI_Comm_split(group_comm, col, row, &col_comm);
 152   /*-------------------------Communication Layer can be used------------------*/
 153
 154   // matrix sizes
 155   m   = m   / size_col;
 156   n   = n   / size_row;
 157   size_t k_a = k / size_row;
 158   size_t k_b = k / size_col;
 159
 160
 161
 162   /*only on the group 0*/
 163   if( group == 0 ) {
 164     matrices_initialisation(&a, &b, &c, m, k_a, k_b, n, row, col);
 165     if( NB_groups > 1 ) res = malloc( sizeof(double)*m*n );
 166   } else matrices_allocation(&a, &b, &c, m, k_a, k_b, n);
 167
 168
 169   /*-------------------Configuration for Summa algorihtm--------------------*/
 170   /*--------------------Allocation of matrices block-------------------------*/
 171   double *a_Summa, *b_Summa;
 172   blocks_initialisation(&a_Summa, &b_Summa, m, Block_size, n);
 173
 174   /*--------------------Communication types for MPI--------------------------*/
 175   MPI_Datatype Block_a;
 176   MPI_Datatype Block_a_local;
 177   MPI_Datatype Block_b;
 178   MPI_Type_vector(m , Block_size, k_a, MPI_DOUBLE, &Block_a);
 179   MPI_Type_vector(m , Block_size, Block_size, MPI_DOUBLE, &Block_a_local);
 180   MPI_Type_vector(Block_size, n, n, MPI_DOUBLE, &Block_b);
 181   MPI_Type_commit(&Block_a);
 182   MPI_Type_commit(&Block_a_local);
 183   MPI_Type_commit(&Block_b);
 184   /*-------------Communication types for MPI are configured------------------*/
 185
 186
 187
 188
 189   MPI_Barrier(my_world);
 190   start_time = MPI_Wtime();
 191   if( NB_groups > 1 ) {
 192     err = MPI_Bcast(a, m*k_a, MPI_DOUBLE, 0, group_line);
 193     if (err != MPI_SUCCESS) {
 194       perror("Error Bcast A\n");
 195       return -1;
 196     }
 197     err = MPI_Bcast(b, n*k_b, MPI_DOUBLE, 0, group_line);
 198     if (err != MPI_SUCCESS) {
 199       perror("Error Bcast B\n");
 200       return -1;
 201     }
 202     MPI_Barrier(my_world);
 203   }
 204   end_time_intern = MPI_Wtime();
 205   communication_time += end_time_intern - start_time;
 206
 207   XBT_INFO( "group %zu NB_block: %zu, NB_groups %zu\n"
 208               ,group,NB_Block, NB_groups);
 209   XBT_INFO(
 210               "m %zu,  k_a %zu,  k_b %zu,  n %zu,  Block_size %zu, "
 211               "group*NB_Block/NB_groups %zu, "
 212               "(group+1)*NB_Block/NB_groups %zu, row %zu,  col %zu,"
 213               "size_row %zu,  size_col %zu\n" ,
 214               m, k_a, k_b, n, Block_size,
 215               group*NB_Block/NB_groups, (group+1)*NB_Block/NB_groups,
 216               row, col, size_row, size_col);
 217
 218
 219   Summa( a, b, c, k_a, n, n, m, k_a, k_b, n, Block_size,
 220           group*NB_Block/NB_groups, (group+1)*NB_Block/NB_groups,
 221           row, col, size_row, size_col, a_Summa, b_Summa, Block_a,
 222           Block_a_local, Block_b, row_comm, col_comm, 0);
 223
 224   /*-------------------------End Summa algorihtm-----------------------------*/
 225
 226   MPI_Comm_rank(group_line, &myrank);
 227
 228   MPI_Barrier(my_world);
 229   start_time_reduce = MPI_Wtime();
 230   if( NB_groups > 1 ) {
 231     // a gather is better?
 232     err = MPI_Reduce(c, res, m*n, MPI_DOUBLE, MPI_SUM, 0, group_line);
 233     if (err != MPI_SUCCESS) {
 234       perror("Error Bcast A\n");
 235       return -1;
 236     }
 237   }else{
 238     double *swap= c;
 239     c = res;
 240     res=swap;
 241   }
 242   MPI_Barrier(my_world);
 243   end_time_reduce = MPI_Wtime();
 244
 245   MPI_Barrier(my_world);
 246   end_time = MPI_Wtime();
 247   time = end_time - start_time;
 248   double reduce_time = end_time_reduce - start_time_reduce;
 249   printf("communication time: %e reduce time: %e seconds, "
 250          "total time: %e seconds\n",communication_time,reduce_time,time);
 251   MPI_Barrier(my_world);
 252
 253 #if CHECK_25D
 254   if(myrank == 0)
 255     check_result(res, a, b, m, n, k_a, k_b, row, col, size_row, size_col);
 256 #endif
 257
 258   // close properly the pragram
 259
 260   MPI_Type_free(&Block_a);
 261   MPI_Type_free(&Block_a_local);
 262   MPI_Type_free(&Block_b);
 263
 264   free(a_Summa);
 265   free(b_Summa);
 266
 267
 268   free( a );
 269   free( b );
 270   if( NB_groups > 1 ) {
 271     free( c );
 272   }
 273   free(res);
 274
 275   MPI_Barrier(MPI_COMM_WORLD);
 276   MPI_Comm_free(&my_world);
 277   MPI_Comm_free(&group_comm);
 278   MPI_Comm_free(&group_line);
 279   MPI_Comm_free(&row_comm);
 280   MPI_Comm_free(&col_comm);
 281   return 0;
 282 }
 283
 284
 285