Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
refMPIJacobi2D.v2.c
Go to the documentation of this file.
1 /*
2 * refMPIJacobi2D.v2.c
3 * Simple 2D Cellular automata. MPI reference code manually developed and optimized
4 *
5 * v2.0
6 * (c) 2007-2014, 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[8] = { MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL };
239  MPI_Status statusV[8];
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_Irecv( &element(matrixCopy, end[0]+1, begin[1]), end[1]-begin[1]+1, MPI_DOUBLE, rankDown, 0, MPI_COMM_WORLD, &request[4] );
280  mpiTestError( ok, "Recv from Down" );
281  }
282  if ( myCoords[0] != 0 ) {
283  ok = MPI_Irecv( &element(matrixCopy, begin[0]-1, begin[1]), end[1]-begin[1]+1, MPI_DOUBLE, rankUp, 1, MPI_COMM_WORLD, &request[5] );
284  mpiTestError( ok, "Recv from Up" );
285  }
286  if ( myCoords[1] != procs[1]-1 ) {
287  ok = MPI_Irecv( &element(matrixCopy, begin[0], end[1]+1), 1, verticalBorderType, rankRight, 2, MPI_COMM_WORLD, &request[6] );
288  mpiTestError( ok, "Recv from Right" );
289  }
290  if ( myCoords[1] != 0 ) {
291  ok = MPI_Irecv( &element(matrixCopy, begin[0], begin[1]-1), 1, verticalBorderType, rankLeft, 3, MPI_COMM_WORLD, &request[7] );
292  mpiTestError( ok, "Recv from Left" );
293  }
294 
295  /* 15. END Communications */
296  MPI_Waitall( 8, request, statusV );
297 
298  } /* END LOOP */
299 
300  // UNROLLED LOOP: LAST ITERATION WITHOUT COMMUNICATION
301  /* 16. COMPUTE LOCAL PART */
302  for (i=begin[0]; i<=end[0]; i++) {
303  for (j=begin[1]; j<=end[1]; j++) {
304  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;
305  }
306  }
307 
308  /* 17. FREE COMMUNICATION PATTERN RESOURCES */
309  ok = MPI_Type_free( &verticalBorderType );
310  mpiTestError( ok, "Freeing vertical border type" );
311 
312 #ifdef DEBUG
313 for (i=begin[0]; i<=end[0]; i++) {
314  for (j=begin[1]; j<=end[1]; j++) {
315  fprintf(stderr,"M[%d][%d] %lf\n", i,j, element(matrix,i,j));
316  }
317 }
318 #endif
319 
320  /* 18. COMMUNICATE ALL LOCAL PIECES */
321  if ( mpiRank == 0 ) {
322  /* 18.1. RANK 0 RECEIVES ALL */
323  int rank;
324  for ( rank=1; rank<mpiNProc; rank++ ) {
325  /* COMPUTE THE POSITION OF THE REMOTE LOCAL PART TO BE RECEIVED */
326  int remoteCoords[2];
327  ok = MPI_Cart_coords(cartesianComm, rank, 2, remoteCoords);
328  mpiTestError( ok, "Getting remote cartesian coordinates in end comm." );
329 
330  int remoteBegin[2];
331  int remoteEnd[2];
332 
333  remoteBegin[0] = 1 + (int)floor( remoteCoords[0]*procRatio[0] );
334  remoteEnd[0] = 1 + (int)floor( (remoteCoords[0]+((rows > procs[0]) ? 1 : 0))*procRatio[0] - ((rows > procs[0]) ? 1 : 0) );
335  remoteBegin[1] = 1 + (int)floor( remoteCoords[1]*procRatio[1] );
336  remoteEnd[1] = 1 + (int)floor( (remoteCoords[1]+((columns > procs[1]) ? 1 : 0))*procRatio[1] - ((columns > procs[1]) ? 1 : 0) );
337 
338 #ifdef DEBUG
339 fprintf(stderr,"CTRL RemoteMatrix %d, [%d:%d][%d:%d]\n", rank, remoteBegin[0], remoteEnd[0], remoteBegin[1], remoteEnd[1]);
340 #endif
341  MPI_Status status;
342  MPI_Datatype remoteType;
343  ok = MPI_Type_hvector( remoteEnd[0]-remoteBegin[0]+1, remoteEnd[1]-remoteBegin[1]+1, (MPI_Aint)((int)sizeof(double)*columns), MPI_DOUBLE, &remoteType );
344  mpiTestError( ok, "Getting remote type in end comm." );
345  ok = MPI_Type_commit( &remoteType );
346  mpiTestError( ok, "Commiting the remote type in end comm." );
347 
348  ok = MPI_Recv( &element( matrix,remoteBegin[0],remoteBegin[1] ), 1, remoteType, rank, rank, MPI_COMM_WORLD, &status );
349  mpiTestError( ok, "Receiving remote parts in end comm." );
350 
351  ok = MPI_Type_free( &remoteType );
352  mpiTestError( ok, "Freeing remote type in end comm." );
353  }
354  }
355  else {
356  /* 18.2. OTHERS SEND THE LOCAL PART */
357 #ifdef DEBUG
358 fprintf(stderr,"CTRL LocalMatrix %d, [%d:%d][%d:%d]\n", mpiRank, begin[0], end[0], begin[1], end[1]);
359 #endif
360  MPI_Datatype localType;
361  ok = MPI_Type_hvector( end[0]-begin[0]+1, end[1]-begin[1]+1, (MPI_Aint)((int)sizeof(double)*columns), MPI_DOUBLE, &localType );
362  mpiTestError( ok, "Getting local type in end comm." );
363  ok = MPI_Type_commit( &localType );
364  mpiTestError( ok, "Commiting the local type in end comm." );
365 
366  ok = MPI_Send( &element( matrix,begin[0],begin[1] ), 1, localType, 0, mpiRank, MPI_COMM_WORLD );
367  mpiTestError( ok, "Sending local part in end comm." );
368 
369  ok = MPI_Type_free( &localType );
370  mpiTestError( ok, "Freeing local type in end comm." );
371  }
372 
373  /* 19. STOP AND PRINT CLOCK */
374  clock = MPI_Wtime() - clock;
375  MPI_Reduce( &clock, &clockReduce, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
376  if (mpiRank==0) printf("Clock_%s: %14.6lf\n", "mainClock", clockReduce);
377 
378  /* 20. WRITE RESULT IN FILE */
379 #ifdef WRITE_RESULT
380  if ( mpiRank == 0 ) {
381  FILE *f;
382  if ((f=fopen("Result.out.dtxt","w")) == NULL) {
383  fprintf(stderr,"Error: Imposible to open output file\n");
384  exit(-1);
385  }
386  for (i=0; i<rows; i++) {
387  for (j=0; j<rows; j++) {
388  fprintf(f,"%014.4lf", element(matrix,i,j) );
389  fprintf(f,"\n");
390  }
391  }
392  fclose(f);
393  }
394 #endif
395 
396  /* 21. END PROGRAM */
397  free( procs );
398  free( matrix );
399  free( matrixCopy );
400  MPI_Comm_free( &cartesianComm );
401  MPI_Finalize();
402 
403  /* 22. RETURN CORRECT */
404  return 0;
405 }
406 
int mpiNProc
Definition: refMPICannon.c:69
#define element(mat, idx1, idx2)
int rank
Definition: SWpar_ref.c:181
double clockReduce
Definition: refMPICannon.c:70
int mpiRank
Definition: refMPICannon.c:68
int main(int argc, char *argv[])
Definition: cannonAsync.c:62
#define mpiTestError(ok, cad)
double clock
int * factors2D(int numProc)