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