46 #define mpiTestError(ok,cad) \
47 if ( ok != MPI_SUCCESS ) { \
48 fprintf(stderr,"RunTime-Error, Rank %d: %s - %d\n", mpiRank, cad, ok);\
67 result = (
int *)calloc(2,
sizeof(
int));
69 fprintf(stderr,
"RunTime-Error, Rank %d: Allocating factors2D result\n",
mpiRank);
72 low = high = (int)sqrt((
double)numProc);
76 while ( product != numProc ) {
77 if (product < numProc) {
79 product = product + low;
83 if (low==1) { high=numProc; product=numProc; }
84 else product = product - high;
98 int main(
int argc,
char *argv[]) {
102 setbuf(stderr, NULL);
103 setbuf(stdout, NULL);
107 int rows, columns, numIter;
110 fprintf(stderr,
"\nUsage: %s <numRows> <numColumns> <numIterations>\n", argv[0]);
113 rows = atoi( argv[1] );
114 columns = atoi( argv[2] );
115 numIter = atoi( argv[3] );
119 ok = MPI_Init( &argc, &argv);
123 MPI_Comm_rank(MPI_COMM_WORLD, &
mpiRank);
124 MPI_Comm_size(MPI_COMM_WORLD, &
mpiNProc);
129 matrix = (
double *)calloc( (
size_t)(rows * columns),
sizeof(
double) );
130 matrixCopy = (
double *)calloc( (
size_t)(rows * columns),
sizeof(
double) );
132 if ( matrix==NULL || matrixCopy==NULL) {
133 fprintf(stderr,
"RunTime-Error, Rank %d: Allocating data structures\n",
mpiRank);
137 #define element(mat, idx1, idx2) (mat[(idx1)*columns+(idx2)])
140 for (i=0; i<rows; i++) {
141 for (j=0; j<columns; j++) {
146 for (i=0; i<rows; i++) {
148 element(matrix, i, columns-1) = 4.0;
150 element(matrixCopy, i, 0) = 3.0;
151 element(matrixCopy, i, columns-1) = 4.0;
153 for (i=0; i<columns; i++) {
155 element(matrix, rows-1, i) = 2.0;
157 element(matrixCopy, 0, i) = 1.0;
158 element(matrixCopy, rows-1, i) = 2.0;
162 MPI_Barrier( MPI_COMM_WORLD );
166 MPI_Comm cartesianComm;
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");
175 int periods[2] = { 0, 0 };
176 ok = MPI_Cart_create(MPI_COMM_WORLD, 2, procs, periods, 0, &cartesianComm);
180 ok = MPI_Cart_coords(cartesianComm,
mpiRank, 2, myCoords);
183 int rankUp=-1000, rankDown=-1000, rankLeft=-1000, rankRight=-1000;
186 if ( myCoords[0] != 0 ) {
187 ok = MPI_Cart_shift(cartesianComm, 0, -1, &foo, &rankUp);
188 mpiTestError( ok,
"Shifting cartesian coordinates up" );
191 if ( myCoords[0] != procs[0]-1 ) {
192 ok = MPI_Cart_shift(cartesianComm, 0, +1, &foo, &rankDown);
193 mpiTestError( ok,
"Shifting cartesian coordinates down" );
196 if ( myCoords[1] != 0 ) {
197 ok = MPI_Cart_shift(cartesianComm, 1, -1, &foo, &rankLeft);
198 mpiTestError( ok,
"Shifting cartesian coordinates left" );
201 if ( myCoords[1] != procs[1]-1 ) {
202 ok = MPI_Cart_shift(cartesianComm, 1, +1, &foo, &rankRight);
203 mpiTestError( ok,
"Shifting cartesian coordinates right" );
207 fprintf(stderr,
"CTRL RANKS: %d (%d,%d,%d,%d)\n",
mpiRank, rankUp, rankDown, rankLeft, rankRight);
212 procRatio[0] = (double)(rows-2)/procs[0];
213 procRatio[1] = (double)(columns-2)/procs[1];
219 fprintf(stderr,
"CTRL %d, myCoords %d,%d\n",
mpiRank, myCoords[0], myCoords[1]);
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) );
228 fprintf(stderr,
"CTRL %d, matrix[%d:%d][%d:%d]\n",
mpiRank, begin[0], end[0], begin[1], end[1]);
232 MPI_Datatype verticalBorderType;
233 ok = MPI_Type_hvector( end[0]-begin[0]+1, 1, (
int)
sizeof(
double)*columns, MPI_DOUBLE, &verticalBorderType );
235 ok = MPI_Type_commit( &verticalBorderType );
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];
243 for (loopIndex = 0; loopIndex < numIter-1; loopIndex++) {
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;
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] );
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] );
261 if ( myCoords[1] != 0 ) {
262 ok = MPI_Isend( &
element(matrix,begin[0], begin[1]), 1, verticalBorderType, rankLeft, 2, MPI_COMM_WORLD, &request[2] );
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] );
271 for (i=begin[0]; i<=end[0]; i++) {
272 for (j=begin[1]; j<=end[1]; j++) {
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] );
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] );
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] );
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] );
296 MPI_Waitall( 8, request, statusV );
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;
309 ok = MPI_Type_free( &verticalBorderType );
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));
324 for ( rank=1; rank<
mpiNProc; rank++ ) {
327 ok = MPI_Cart_coords(cartesianComm, rank, 2, remoteCoords);
328 mpiTestError( ok,
"Getting remote cartesian coordinates in end comm." );
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) );
339 fprintf(stderr,
"CTRL RemoteMatrix %d, [%d:%d][%d:%d]\n", rank, remoteBegin[0], remoteEnd[0], remoteBegin[1], remoteEnd[1]);
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 );
345 ok = MPI_Type_commit( &remoteType );
346 mpiTestError( ok,
"Commiting the remote type in end comm." );
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." );
351 ok = MPI_Type_free( &remoteType );
358 fprintf(stderr,
"CTRL LocalMatrix %d, [%d:%d][%d:%d]\n",
mpiRank, begin[0], end[0], begin[1], end[1]);
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 );
363 ok = MPI_Type_commit( &localType );
364 mpiTestError( ok,
"Commiting the local type in end comm." );
366 ok = MPI_Send( &
element( matrix,begin[0],begin[1] ), 1, localType, 0,
mpiRank, MPI_COMM_WORLD );
369 ok = MPI_Type_free( &localType );
375 MPI_Reduce( &
clock, &
clockReduce, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
382 if ((f=fopen(
"Result.out.dtxt",
"w")) == NULL) {
383 fprintf(stderr,
"Error: Imposible to open output file\n");
386 for (i=0; i<rows; i++) {
387 for (j=0; j<rows; j++) {
388 fprintf(f,
"%014.4lf",
element(matrix,i,j) );
400 MPI_Comm_free( &cartesianComm );
#define element(mat, idx1, idx2)
int main(int argc, char *argv[])
#define mpiTestError(ok, cad)
int * factors2D(int numProc)