Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
refMPIGaussSeidel2D.v1.c
Go to the documentation of this file.
1 /*
2 * refMPIGaussSeidel2D.c
3 * MPI reference code manually developed and optimized
4 * Stencil code: Gauss-Seidel 2D for the heat equation.
5 *
6 * v1.0
7 * (c) 2013, Arturo González Escribano
8 */
9 
10 /*
11  * <license>
12  *
13  * Hitmap v1.2
14  *
15  * This software is provided to enhance knowledge and encourage progress in the scientific
16  * community. It should be used only for research and educational purposes. Any reproduction
17  * or use for commercial purpose, public redistribution, in source or binary forms, with or
18  * without modifications, is NOT ALLOWED without the previous authorization of the copyright
19  * holder. The origin of this software must not be misrepresented; you must not claim that you
20  * wrote the original software. If you use this software for any purpose (e.g. publication),
21  * a reference to the software package and the authors must be included.
22  *
23  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS "AS IS" AND ANY
24  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
25  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
26  * THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
27  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
28  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
29  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
30  * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32  *
33  * Copyright (c) 2007-2015, Trasgo Group, Universidad de Valladolid.
34  * All rights reserved.
35  *
36  * More information on http://trasgo.infor.uva.es/
37  *
38  * </license>
39 */
40 
41 #include<stdio.h>
42 #include<math.h>
43 #include<stdlib.h>
44 #include<mpi.h>
45 
46 
47 #define mpiTestError(ok,cad) \
48  if ( ok != MPI_SUCCESS ) { \
49  fprintf(stderr,"RunTime-Error, Rank %d: %s - %d\n", mpiRank, cad, ok);\
50  exit(-1); \
51  }
52 
53 
54 int mpiRank;
57 
58 
59 /*
60  * UTILITY: DECOMPOSE A NUMBER IN 2 PRIME FACTORS AS NEAR AS POSSIBLE TO THE SQUARE ROOT
61  */
62 int * factors2D(int numProc) {
63  int *result;
64  int low, high;
65  int product;
66 
67  /* 1. ALLOCATE MEMORY FOR result */
68  result = (int *)calloc(2, sizeof(int));
69  if ( result == NULL )
70  fprintf(stderr,"RunTime-Error, Rank %d: Allocating factors2D result\n", mpiRank);
71 
72  /* 2. START AS NEAR AS POSSIBLE TO THE SQUARE ROOT */
73  low = high = (int)sqrt((double)numProc);
74  product = low*high;
75 
76  /* 3. MOVE UP/DOWN TO SEARCH THE NEAREST FACTORS */
77  while ( product != numProc ) {
78  if (product < numProc) {
79  high++;
80  product = product + low;
81  }
82  else {
83  low--;
84  if (low==1) { high=numProc; product=numProc; }
85  else product = product - high;
86  }
87  }
88 
89  /* 4. RETURN RESULT */
90  result[0] = high;
91  result[1] = low;
92  return result;
93 }
94 
95 
96 /*
97  * MAIN PROGRAM
98  */
99 int main(int argc, char *argv[]) {
100  int ok;
101  int i,j;
102 #ifdef DEBUG
103  setbuf(stderr, NULL);
104  setbuf(stdout, NULL);
105 #endif
106 
107  /* 0. READ PARAMETERS */
108  int rows, columns, numIter;
109 
110  if ( argc != 4 ) {
111  fprintf(stderr, "\nUsage: %s <numRows> <numColumns> <numIterations>\n", argv[0]);
112  exit(EXIT_FAILURE);
113  }
114  rows = atoi( argv[1] );
115  columns = atoi( argv[2] );
116  numIter = atoi( argv[3] );
117 
118 
119  /* 1. INITIALIZE MPI */
120  ok = MPI_Init( &argc, &argv);
121  mpiTestError( ok, "Initializing MPI" );
122 
123  /* 2. GET BASIC GROUP PARAMETERS */
124  MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank);
125  MPI_Comm_size(MPI_COMM_WORLD, &mpiNProc);
126 
127  /* 3. DEFINE MAIN DATA STRUCTURE */
128  double *matrix;
129  matrix = (double *)calloc( (size_t)(rows * columns), sizeof(double) );
130 
131  if ( matrix==NULL ) {
132  fprintf(stderr,"RunTime-Error, Rank %d: Allocating data structures\n", mpiRank);
133  exit(2);
134  }
135 
136 #define element(mat, idx1, idx2) (mat[(idx1)*columns+(idx2)])
137 
138  /* 4. INITIAL MATRIX VALUES */
139  for (i=0; i<rows; i++) {
140  for (j=0; j<columns; j++) {
141  element(matrix, i, j) = 0.0;
142  }
143  }
144  /* 5. FIXED BOUNDARY VALUES */
145  for (i=0; i<rows; i++) {
146  element(matrix, i, 0) = 3.0;
147  element(matrix, i, columns-1) = 4.0;
148  }
149  for (i=0; i<columns; i++) {
150  element(matrix, 0, i) = 1.0;
151  element(matrix, rows-1, i) = 2.0;
152  }
153 
154  /* 6. START CLOCK */
155  MPI_Barrier( MPI_COMM_WORLD );
156  clock = MPI_Wtime();
157 
158  /* 7. CARTESIAN TOPOLOGY */
159  MPI_Comm cartesianComm;
160  int *procs = factors2D( mpiNProc );
161 
162  int periods[2] = { 0, 0 };
163  ok = MPI_Cart_create(MPI_COMM_WORLD, 2, procs, periods, 0, &cartesianComm);
164  mpiTestError( ok, "Creating cartesian toology" );
165 
166  int myCoords[2];
167  ok = MPI_Cart_coords(cartesianComm, mpiRank, 2, myCoords);
168  mpiTestError( ok, "Getting cartesian coordinates" );
169 
170  int rankUp=-1000, rankDown=-1000, rankLeft=-1000, rankRight=-1000;
171  int foo;
172 
173  if ( myCoords[0] != 0 ) {
174  ok = MPI_Cart_shift(cartesianComm, 0, -1, &foo, &rankUp);
175  mpiTestError( ok, "Shifting cartesian coordinates up" );
176  }
177 
178  if ( myCoords[0] != procs[0]-1 ) {
179  ok = MPI_Cart_shift(cartesianComm, 0, +1, &foo, &rankDown);
180  mpiTestError( ok, "Shifting cartesian coordinates down" );
181  }
182 
183  if ( myCoords[1] != 0 ) {
184  ok = MPI_Cart_shift(cartesianComm, 1, -1, &foo, &rankLeft);
185  mpiTestError( ok, "Shifting cartesian coordinates left" );
186  }
187 
188  if ( myCoords[1] != procs[1]-1 ) {
189  ok = MPI_Cart_shift(cartesianComm, 1, +1, &foo, &rankRight);
190  mpiTestError( ok, "Shifting cartesian coordinates right" );
191  }
192 
193  /* 8. SANITY CHECK (NO MORE PROCs THAN DATA ELEMENTS IN BOTH DIMENSIONS) */
194  if ( procs[0] > rows-2 || procs[1] > columns-2 ) {
195  if ( mpiRank == 0 )
196  fprintf(stderr, "Error: This code does not support partitioning less rows/columns than processors in an axis\n");
197  exit(2);
198  }
199 
200 #ifdef DEBUG
201 fprintf(stderr,"CTRL RANKS: %d (%d,%d,%d,%d)\n", mpiRank, rankUp, rankDown, rankLeft, rankRight);
202 #endif
203 
204  /* 9. COMPUTE WHICH IS THE LOCAL PART OF THE MATRIX */
205  double procRatio[2];
206  procRatio[0] = (double)(rows-2)/procs[0];
207  procRatio[1] = (double)(columns-2)/procs[1];
208 
209  int begin[2];
210  int end[2];
211 
212 #ifdef DEBUG
213 fprintf(stderr,"CTRL %d, myCoords %d,%d\n", mpiRank, myCoords[0], myCoords[1]);
214 #endif
215 
216  begin[0] = 1 + (int)floor( myCoords[0]*procRatio[0] );
217  end[0] = 1 + (int)floor( (myCoords[0]+((rows > procs[0]) ? 1 : 0))*procRatio[0] - ((rows > procs[0]) ? 1 : 0) );
218  begin[1] = 1 + (int)floor( myCoords[1]*procRatio[1] );
219  end[1] = 1 + (int)floor( (myCoords[1]+((columns > procs[1]) ? 1 : 0))*procRatio[1] - ((columns > procs[1]) ? 1 : 0) );
220 
221 #ifdef DEBUG
222 fprintf(stderr,"CTRL %d, matrix[%d:%d][%d:%d]\n", mpiRank, begin[0], end[0], begin[1], end[1]);
223 #endif
224 
225  /* 10. COMMIT MPI TYPES */
226  MPI_Datatype verticalBorderType;
227  ok = MPI_Type_hvector( end[0]-begin[0]+1, 1, (int)sizeof(double)*columns, MPI_DOUBLE, &verticalBorderType );
228  mpiTestError( ok, "Creating the columns type" );
229  ok = MPI_Type_commit( &verticalBorderType );
230  mpiTestError( ok, "Commiting the columns type" );
231 
232  MPI_Request request[4];
233  MPI_Status status;
234 
235  /* 11. LOOP */
236  int loopIndex;
237  for (loopIndex = 0; loopIndex < numIter; loopIndex++) {
238  /* 12. COMMUNICATE BORDERS: RECV IN MATRIX (Inner dependences) */
239  if ( myCoords[0] != 0 ) {
240  ok = MPI_Recv( &element(matrix, begin[0]-1, begin[1]), end[1]-begin[1]+1, MPI_DOUBLE, rankUp, 1, MPI_COMM_WORLD, &status );
241  mpiTestError( ok, "Recv from Up" );
242  }
243  if ( myCoords[1] != 0 ) {
244  ok = MPI_Recv( &element(matrix, begin[0], begin[1]-1), 1, verticalBorderType, rankLeft, 3, MPI_COMM_WORLD, &status );
245  mpiTestError( ok, "Recv from Left" );
246  }
247 
248  /* 13. COMPUTE LOCAL PART */
249  for (i=begin[0]; i<=end[0]; i++) {
250  for (j=begin[1]; j<=end[1]; j++) {
251  element( matrix,i,j ) = ( element( matrix,i-1,j ) + element( matrix,i+1,j ) + element( matrix,i,j-1 ) + element( matrix,i,j+1 ) ) / 4;
252  }
253  }
254 
255  /* 14. COMMUNICATE BORDERS: SEND MATRIX VALUES */
256  if ( myCoords[0] != 0 ) {
257  ok = MPI_Isend( &element(matrix, begin[0], begin[1]), end[1]-begin[1]+1, MPI_DOUBLE, rankUp, 0, MPI_COMM_WORLD, &request[0] );
258  mpiTestError( ok, "Send Up" );
259  }
260  if ( myCoords[0] != procs[0]-1 ) {
261  ok = MPI_Isend( &element(matrix, end[0], begin[1]), end[1]-begin[1]+1, MPI_DOUBLE, rankDown, 1, MPI_COMM_WORLD, &request[1] );
262  mpiTestError( ok, "Send Down" );
263  }
264  if ( myCoords[1] != 0 ) {
265  ok = MPI_Isend( &element(matrix,begin[0], begin[1]), 1, verticalBorderType, rankLeft, 2, MPI_COMM_WORLD, &request[2] );
266  mpiTestError( ok, "Send Left" );
267  }
268  if ( myCoords[1] != procs[1]-1 ) {
269  ok = MPI_Isend( &element(matrix,begin[0],end[1]), 1, verticalBorderType, rankRight, 3, MPI_COMM_WORLD, &request[3] );
270  mpiTestError( ok, "Send Right" );
271  }
272 
273  /* 15. COMMUNICATE BORDERS: RECV IN MATRIX (External loop carried dependences) */
274  if ( myCoords[0] != procs[0]-1 ) {
275  ok = MPI_Recv( &element(matrix, end[0]+1, begin[1]), end[1]-begin[1]+1, MPI_DOUBLE, rankDown, 0, MPI_COMM_WORLD, &status );
276  mpiTestError( ok, "Recv from Down" );
277  }
278  if ( myCoords[1] != procs[1]-1 ) {
279  ok = MPI_Recv( &element(matrix, begin[0], end[1]+1), 1, verticalBorderType, rankRight, 2, MPI_COMM_WORLD, &status );
280  mpiTestError( ok, "Recv from Right" );
281  }
282 
283  /* 16. END Isends */
284  if ( myCoords[0] != 0 ) {
285  MPI_Wait( &request[0], &status );
286  }
287  if ( myCoords[0] != procs[0]-1 ) {
288  MPI_Wait( &request[1], &status );
289  }
290  if ( myCoords[1] != 0 ) {
291  MPI_Wait( &request[2], &status );
292  }
293  if ( myCoords[1] != procs[1]-1 ) {
294  MPI_Wait( &request[3], &status );
295  }
296  } /* END LOOP */
297  ok = MPI_Type_free( &verticalBorderType );
298  mpiTestError( ok, "Freeing vertical border type" );
299 
300 #ifdef DEBUG
301 for (i=begin[0]; i<=end[0]; i++) {
302  for (j=begin[1]; j<=end[1]; j++) {
303  fprintf(stderr,"M[%d][%d] %lf\n", i,j, element(matrix,i,j));
304  }
305 }
306 #endif
307 
308  /* 17. COMMUNICATE ALL LOCAL PIECES */
309  if ( mpiRank == 0 ) {
310  /* 17.1. RANK 0 RECEIVES ALL */
311  int rank;
312  for ( rank=1; rank<mpiNProc; rank++ ) {
313  /* COMPUTE THE REMOTE LOCAL PART */
314  int remoteCoords[2];
315  ok = MPI_Cart_coords(cartesianComm, rank, 2, remoteCoords);
316  mpiTestError( ok, "Getting remote cartesian coordinates in end comm." );
317 
318  int remoteBegin[2];
319  int remoteEnd[2];
320 
321  remoteBegin[0] = 1 + (int)floor( remoteCoords[0]*procRatio[0] );
322  remoteEnd[0] = 1 + (int)floor( (remoteCoords[0]+((rows > procs[0]) ? 1 : 0))*procRatio[0] - ((rows > procs[0]) ? 1 : 0) );
323  remoteBegin[1] = 1 + (int)floor( remoteCoords[1]*procRatio[1] );
324  remoteEnd[1] = 1 + (int)floor( (remoteCoords[1]+((columns > procs[1]) ? 1 : 0))*procRatio[1] - ((columns > procs[1]) ? 1 : 0) );
325 
326 #ifdef DEBUG
327 fprintf(stderr,"CTRL RemoteMatrix %d, [%d:%d][%d:%d]\n", rank, remoteBegin[0], remoteEnd[0], remoteBegin[1], remoteEnd[1]);
328 #endif
329 
330  MPI_Datatype remoteType;
331  ok = MPI_Type_hvector( remoteEnd[0]-remoteBegin[0]+1, remoteEnd[1]-remoteBegin[1]+1, (int)sizeof(double)*columns, MPI_DOUBLE, &remoteType );
332  mpiTestError( ok, "Getting remote type in end comm." );
333  ok = MPI_Type_commit( &remoteType );
334  mpiTestError( ok, "Commiting the remote type in end comm." );
335 
336  ok = MPI_Recv( &element( matrix,remoteBegin[0],remoteBegin[1] ), 1, remoteType, rank, rank, MPI_COMM_WORLD, &status );
337  mpiTestError( ok, "Receiving remote parts in end comm." );
338 
339  ok = MPI_Type_free( &remoteType );
340  mpiTestError( ok, "Freeing remote type in end comm." );
341  }
342  }
343  else {
344  /* 17.2. OTHERS SEND THE LOCAL PART */
345 #ifdef DEBUG
346 fprintf(stderr,"CTRL LocalMatrix %d, [%d:%d][%d:%d]\n", mpiRank, begin[0], end[0], begin[1], end[1]);
347 #endif
348  MPI_Datatype localType;
349  ok = MPI_Type_hvector( end[0]-begin[0]+1, end[1]-begin[1]+1, (int)sizeof(double)*columns, MPI_DOUBLE, &localType );
350  mpiTestError( ok, "Getting local type in end comm." );
351  ok = MPI_Type_commit( &localType );
352  mpiTestError( ok, "Commiting the local type in end comm." );
353 
354  ok = MPI_Send( &element( matrix,begin[0],begin[1] ), 1, localType, 0, mpiRank, MPI_COMM_WORLD );
355  mpiTestError( ok, "Sending local part in end comm." );
356 
357  ok = MPI_Type_free( &localType );
358  mpiTestError( ok, "Freeing local type in end comm." );
359  }
360 
361  /* 18. STOP AND PRINT CLOCK */
362  clock = MPI_Wtime() - clock;
363  MPI_Reduce( &clock, &clockReduce, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
364  if (mpiRank==0) printf("Clock_%s: %14.6lf\n", "mainClock", clockReduce);
365 
366  /* 19. WRITE RESULT IN FILE */
367 #ifdef WRITE_RESULT
368  if ( mpiRank == 0 ) {
369  FILE *f;
370  if ((f=fopen("Result.out.dtxt","w")) == NULL) {
371  fprintf(stderr,"Error: Imposible to open output file\n");
372  exit(-1);
373  }
374  for (i=0; i<rows; i++) {
375  for (j=0; j<rows; j++) {
376  fprintf(f,"%014.4lf", element(matrix,i,j) );
377  fprintf(f,"\n");
378  }
379  }
380  fclose(f);
381  }
382 #endif
383 
384  /* 20. END PROGRAM */
385  MPI_Finalize();
386 
387  /* 21. RETURN CORRECT */
388  return 0;
389 }
390 
#define mpiTestError(ok, cad)
int mpiNProc
Definition: refMPICannon.c:69
int rank
Definition: SWpar_ref.c:181
double clockReduce
Definition: refMPICannon.c:70
#define element(mat, idx1, idx2)
int mpiRank
Definition: refMPICannon.c:68
int main(int argc, char *argv[])
Definition: cannonAsync.c:62
double clock
int * factors2D(int numProc)