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