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