50 #define mpiTestError(ok,cad) \
51 if ( ok != MPI_SUCCESS ) { \
52 fprintf(stderr,"RunTime-Error, Rank %d: %s - %d\n", mpiRank, cad, ok);\
71 result = (
int *)calloc(2,
sizeof(
int));
73 fprintf(stderr,
"RunTime-Error, Rank %d: Allocating factors2D result\n",
mpiRank);
76 low = high = (int)sqrt((
double)numProc);
80 while ( product != numProc ) {
81 if (product < numProc) {
83 product = product + low;
87 if (low==1) { high=numProc; product=numProc; }
88 else product = product - high;
102 int main(
int argc,
char *argv[]) {
106 setbuf(stderr, NULL);
107 setbuf(stdout, NULL);
111 int rows, columns, numIter;
114 fprintf(stderr,
"\nUsage: %s <numRows> <numColumns> <numIterations>\n", argv[0]);
117 rows = atoi( argv[1] );
118 columns = atoi( argv[2] );
119 numIter = atoi( argv[3] );
123 ok = MPI_Init( &argc, &argv);
127 MPI_Comm_rank(MPI_COMM_WORLD, &
mpiRank);
128 MPI_Comm_size(MPI_COMM_WORLD, &
mpiNProc);
131 MPI_Barrier( MPI_COMM_WORLD );
135 MPI_Comm cartesianComm;
139 if ( procs[0] > rows-2 || procs[1] > columns-2 ) {
140 if (
mpiRank == 0 ) fprintf(stderr,
"Error: This code does not support partitioning less rows/columns than processors in that axis\n");
144 int periods[2] = { 0, 0 };
145 ok = MPI_Cart_create(MPI_COMM_WORLD, 2, procs, periods, 0, &cartesianComm);
149 ok = MPI_Cart_coords(cartesianComm,
mpiRank, 2, myCoords);
152 int rankUp=-1000, rankDown=-1000, rankLeft=-1000, rankRight=-1000;
155 if ( myCoords[0] != 0 ) {
156 ok = MPI_Cart_shift(cartesianComm, 0, -1, &foo, &rankUp);
157 mpiTestError( ok,
"Shifting cartesian coordinates up" );
160 if ( myCoords[0] != procs[0]-1 ) {
161 ok = MPI_Cart_shift(cartesianComm, 0, +1, &foo, &rankDown);
162 mpiTestError( ok,
"Shifting cartesian coordinates down" );
165 if ( myCoords[1] != 0 ) {
166 ok = MPI_Cart_shift(cartesianComm, 1, -1, &foo, &rankLeft);
167 mpiTestError( ok,
"Shifting cartesian coordinates left" );
170 if ( myCoords[1] != procs[1]-1 ) {
171 ok = MPI_Cart_shift(cartesianComm, 1, +1, &foo, &rankRight);
172 mpiTestError( ok,
"Shifting cartesian coordinates right" );
176 fprintf(stderr,
"CTRL RANKS: %d (%d,%d,%d,%d)\n",
mpiRank, rankUp, rankDown, rankLeft, rankRight);
181 procRatio[0] = (double)(rows-2)/procs[0];
182 procRatio[1] = (double)(columns-2)/procs[1];
188 fprintf(stderr,
"CTRL %d, myCoords %d,%d\n",
mpiRank, myCoords[0], myCoords[1]);
191 begin[0] = 1 + (int)floor( myCoords[0]*procRatio[0] );
192 end[0] = 1 + (int)floor( (myCoords[0]+((rows > procs[0]) ? 1 : 0))*procRatio[0] - ((rows > procs[0]) ? 1 : 0) );
193 begin[1] = 1 + (int)floor( myCoords[1]*procRatio[1] );
194 end[1] = 1 + (int)floor( (myCoords[1]+((columns > procs[1]) ? 1 : 0))*procRatio[1] - ((columns > procs[1]) ? 1 : 0) );
197 fprintf(stderr,
"CTRL %d, matrix[%d:%d][%d:%d]\n",
mpiRank, begin[0], end[0], begin[1], end[1]);
201 int size[2] = { end[0]-begin[0]+1+2, end[1]-begin[1]+1+2 };
204 matrix = (
double *)calloc( (
size_t)(size[0] * size[1]),
sizeof(
double) );
205 matrixCopy = (
double *)calloc( (
size_t)(size[0] * size[1]),
sizeof(
double) );
207 if ( matrix==NULL || matrixCopy==NULL) {
208 fprintf(stderr,
"RunTime-Error, Rank %d: Allocating data structures\n",
mpiRank);
212 #define elementGlobal(mat, idx1, idx2) (mat[(idx1)*columns+(idx2)])
213 #define element(mat, idx1, idx2) (mat[(idx1)*size[1]+(idx2)])
216 for (i=1; i<size[0]; i++) {
217 for (j=1; j<size[1]; j++) {
222 if ( myCoords[1] == 0 ) {
223 for (i=0; i<size[0]; i++) {
225 element(matrixCopy, i, 0) = 3.0;
228 if ( myCoords[1] == procs[1]-1 ) {
229 for (i=0; i<size[0]; i++) {
230 element(matrix, i, size[1]-1) = 4.0;
231 element(matrixCopy, i, size[1]-1) = 4.0;
234 if ( myCoords[0] == 0 ) {
235 for (i=0; i<size[1]; i++) {
237 element(matrixCopy, 0, i) = 1.0;
240 if ( myCoords[0] == procs[0]-1 ) {
241 for (i=0; i<size[1]; i++) {
242 element(matrix, size[0]-1, i) = 2.0;
243 element(matrixCopy, size[0]-1, i) = 2.0;
248 MPI_Datatype verticalBorderType;
249 ok = MPI_Type_hvector( size[0]-2, 1, (
int)
sizeof(
double)*size[1], MPI_DOUBLE, &verticalBorderType );
251 ok = MPI_Type_commit( &verticalBorderType );
254 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 };
255 MPI_Status statusV[8];
259 for (loopIndex = 0; loopIndex < numIter-1; loopIndex++) {
262 for (i=1; i<size[0]-1; i++) {
263 for (j=1; j<size[1]-1; j++) {
264 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;
269 if ( myCoords[0] != 0 ) {
270 ok = MPI_Isend( &
element(matrix, 1, 1), size[1]-2, MPI_DOUBLE, rankUp, 0, MPI_COMM_WORLD, &request[0] );
273 if ( myCoords[0] != procs[0]-1 ) {
274 ok = MPI_Isend( &
element(matrix, size[0]-2, 1), size[1]-2, MPI_DOUBLE, rankDown, 1, MPI_COMM_WORLD, &request[1] );
277 if ( myCoords[1] != 0 ) {
278 ok = MPI_Isend( &
element(matrix, 1, 1), 1, verticalBorderType, rankLeft, 2, MPI_COMM_WORLD, &request[2] );
281 if ( myCoords[1] != procs[1]-1 ) {
282 ok = MPI_Isend( &
element(matrix, 1, size[1]-2), 1, verticalBorderType, rankRight, 3, MPI_COMM_WORLD, &request[3] );
287 for (i=1; i<size[0]-1; i++) {
288 for (j=1; j<size[1]-1; j++) {
294 if ( myCoords[0] != procs[0]-1 ) {
295 ok = MPI_Irecv( &
element(matrixCopy, size[0]-1, 1), size[1]-2, MPI_DOUBLE, rankDown, 0, MPI_COMM_WORLD, &request[4] );
298 if ( myCoords[0] != 0 ) {
299 ok = MPI_Irecv( &
element(matrixCopy, 0, 1), size[1]-2, MPI_DOUBLE, rankUp, 1, MPI_COMM_WORLD, &request[5] );
302 if ( myCoords[1] != procs[1]-1 ) {
303 ok = MPI_Irecv( &
element(matrixCopy, 1, size[1]-1), 1, verticalBorderType, rankRight, 2, MPI_COMM_WORLD, &request[6] );
306 if ( myCoords[1] != 0 ) {
307 ok = MPI_Irecv( &
element(matrixCopy, 1, 0), 1, verticalBorderType, rankLeft, 3, MPI_COMM_WORLD, &request[7] );
312 MPI_Waitall( 8, request, statusV );
318 for (i=1; i<size[0]-1; i++) {
319 for (j=1; j<size[1]-1; j++) {
320 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;
325 ok = MPI_Type_free( &verticalBorderType );
329 for (i=1; i<size[0]-1; i++) {
330 for (j=1; j<size[1]-1; j++) {
331 fprintf(stderr,
"M[%d][%d] %lf\n", i,j,
element(matrix,i,j));
336 double *matrixResult = NULL;
340 matrixResult = (
double *)calloc( (
size_t)(rows * columns),
sizeof(
double) );
343 for (i=0; i<size[0]; i++) {
344 for (j=0; j<size[1]; j++) {
351 for ( rank=1; rank<
mpiNProc; rank++ ) {
354 ok = MPI_Cart_coords(cartesianComm, rank, 2, remoteCoords);
355 mpiTestError( ok,
"Getting remote cartesian coordinates in end comm." );
360 remoteBegin[0] = 1 + (int)floor( remoteCoords[0]*procRatio[0] );
361 remoteEnd[0] = 1 + (int)floor( (remoteCoords[0]+((rows > procs[0]) ? 1 : 0))*procRatio[0] - ((rows > procs[0]) ? 1 : 0) );
362 remoteBegin[1] = 1 + (int)floor( remoteCoords[1]*procRatio[1] );
363 remoteEnd[1] = 1 + (int)floor( (remoteCoords[1]+((columns > procs[1]) ? 1 : 0))*procRatio[1] - ((columns > procs[1]) ? 1 : 0) );
366 if ( remoteBegin[0] == 1 ) remoteBegin[0]--;
367 if ( remoteBegin[1] == 1 ) remoteBegin[1]--;
368 if ( remoteEnd[0] == rows-2 ) remoteEnd[0]++;
369 if ( remoteEnd[1] == columns-2 ) remoteEnd[1]++;
372 fprintf(stderr,
"CTRL RemoteMatrix %d, [%d:%d][%d:%d]\n", rank, remoteBegin[0], remoteEnd[0], remoteBegin[1], remoteEnd[1]);
375 MPI_Datatype remoteType;
376 ok = MPI_Type_hvector( remoteEnd[0]-remoteBegin[0]+1, remoteEnd[1]-remoteBegin[1]+1, (MPI_Aint)((
int)
sizeof(
double)*columns), MPI_DOUBLE, &remoteType );
378 ok = MPI_Type_commit( &remoteType );
379 mpiTestError( ok,
"Commiting the remote type in end comm." );
381 ok = MPI_Recv( &
elementGlobal( matrixResult,remoteBegin[0],remoteBegin[1] ), 1, remoteType, rank, rank, MPI_COMM_WORLD, &status );
382 mpiTestError( ok,
"Receiving remote parts in end comm." );
384 ok = MPI_Type_free( &remoteType );
391 fprintf(stderr,
"CTRL LocalMatrix %d, [%d:%d][%d:%d]\n",
mpiRank, begin[0], end[0], begin[1], end[1]);
397 localBegin[0] = ( begin[0] == 1 ) ? 0 : 1;
398 localBegin[1] = ( begin[1] == 1 ) ? 0 : 1;
399 localEnd[0] = ( end[0] == rows-2 ) ? size[0]-1 : size[0]-2;
400 localEnd[1] = ( end[1] == columns-2 ) ? size[1]-1 : size[1]-2;
402 MPI_Datatype localType;
403 ok = MPI_Type_hvector( localEnd[0]-localBegin[0]+1, localEnd[1]-localBegin[1]+1, (
int)
sizeof(
double)*size[1], MPI_DOUBLE, &localType );
405 ok = MPI_Type_commit( &localType );
406 mpiTestError( ok,
"Commiting the local type in end comm." );
408 ok = MPI_Send( &
element( matrix,localBegin[0],localBegin[1] ), 1, localType, 0,
mpiRank, MPI_COMM_WORLD );
411 ok = MPI_Type_free( &localType );
417 MPI_Reduce( &
clock, &
clockReduce, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
424 if ((f=fopen(
"Result.out.dtxt",
"w")) == NULL) {
425 fprintf(stderr,
"Error: Imposible to open output file\n");
428 for (i=0; i<rows; i++) {
429 for (j=0; j<rows; j++) {
442 if (
mpiRank == 0 ) free( matrixResult );
443 MPI_Comm_free( &cartesianComm );
#define element(mat, idx1, idx2)
#define elementGlobal(mat, idx1, idx2)
int main(int argc, char *argv[])
int * factors2D(int numProc)
#define mpiTestError(ok, cad)