Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
shunt a warning, quite dirtily
[simgrid.git] / examples / smpi / MM / Summa.c
1 /*!
2  * Classical Block Matrix Multiplication example
3  *
4  * Authors: Quintin Jean-Noël
5  */
6 #include "Matrix_init.h"
7 #include "Summa.h"
8 #include "timer.h"
9 #include "xbt/log.h"
10  XBT_LOG_NEW_DEFAULT_CATEGORY(MM_Summa,
11                              "Messages specific for this msg example");
12
13 inline double Summa(
14                      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   struct timespec start_time, end_time; //time mesure
41   struct timespec start_time_intern, end_time_intern; //time mesure
42
43
44
45
46   get_time(&start_time);
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     get_time(&start_time_intern);
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     get_time(&end_time_intern);
129     communication_time += get_timediff(&start_time_intern,&end_time_intern);
130
131     MPI_Barrier(row_comm);
132     MPI_Barrier(col_comm);
133     get_time(&start_time_intern);
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     get_time(&end_time_intern);
146     computation_time += get_timediff(&start_time_intern,&end_time_intern);
147
148   }
149   MPI_Barrier(row_comm);
150   MPI_Barrier(col_comm);
151
152   get_time(&end_time);
153   time = get_timediff(&start_time,&end_time);
154   printf("communication time: %le nanoseconds, "
155          "computation time: %le nanoseconds\n",
156          communication_time, computation_time);
157
158
159   return time;
160 }