51 #define MATRIX_ROWS 10
52 #define MATRIX_COLUMNS 20
54 #define mpiTestError(ok,cad) \
55 if ( ok != MPI_SUCCESS ) { \
56 fprintf(stderr,"RunTime-Error, Rank %d: %s - %d\n", mpiRank, cad, ok);\
61 #define elemA(i,j) A[(i)*MATRIX_COLUMNS+(j)]
62 #define elemCopyA(i,j) copyA[(i)*MATRIX_COLUMNS+(j)]
63 #define elemB(i,j) B[(i)*MATRIX_ROWS+(j)]
64 #define elemCopyB(i,j) copyB[(i)*MATRIX_ROWS+(j)]
65 #define elemC(i,j) C[(i)*MATRIX_ROWS+(j)]
73 int main(
int argc,
char *argv[]) {
82 ok = MPI_Init( &argc, &argv);
86 MPI_Comm_rank(MPI_COMM_WORLD, &
mpiRank);
87 MPI_Comm_size(MPI_COMM_WORLD, &
mpiNProc);
91 if ( sqProc != floor(sqProc) ) {
93 fprintf(stderr,
"Error: This code does not support a number of procs. without a perfect square root: %d\n",
mpiNProc);
100 fprintf(stderr,
"Error: This code does not support matrix sizes which are not integer divisible by the square root of the number of processors\n");
112 if ( A==NULL || B==NULL || C==NULL ) {
113 fprintf(stderr,
"Error, Rank %d: Allocating data structures\n",
mpiRank);
129 elemB(i,j) = (i==j) ? 1.0 : 0.0;
149 MPI_Barrier( MPI_COMM_WORLD );
153 MPI_Comm cartesianComm;
157 int periods[2] = { 1, 1 };
158 ok = MPI_Cart_create(MPI_COMM_WORLD, 2, procs, periods, 0, &cartesianComm);
162 ok = MPI_Cart_coords(cartesianComm,
mpiRank, 2, myCoords);
165 int rankUp=-1000, rankDown=-1000, rankLeft=-1000, rankRight=-1000;
166 int rankSendIniA=-1000, rankSendIniB=-1000;
167 int rankRecvIniA=-1000, rankRecvIniB=-1000;
169 ok = MPI_Cart_shift(cartesianComm, 0, +1, &rankUp, &rankDown);
170 mpiTestError( ok,
"Shifting cartesian coordinates vertical" );
172 ok = MPI_Cart_shift(cartesianComm, 1, +1, &rankLeft, &rankRight);
173 mpiTestError( ok,
"Shifting cartesian coordinates horizontal" );
175 ok = MPI_Cart_shift(cartesianComm, 1, -myCoords[0], &rankRecvIniA, &rankSendIniA);
176 mpiTestError( ok,
"Shifting cartesian coordinates ini A" );
178 ok = MPI_Cart_shift(cartesianComm, 0, -myCoords[1], &rankRecvIniB, &rankSendIniB);
179 mpiTestError( ok,
"Shifting cartesian coordinates ini B" );
182 fprintf(stderr,
"CTRL %d RANKS: (%d,%d,%d,%d)\n",
mpiRank, rankUp, rankDown, rankLeft, rankRight);
183 fprintf(stderr,
"CTRL %d INI RANKS: A(send %d, recv %d) B(send %d, recv %d)\n",
mpiRank, rankSendIniA, rankRecvIniA, rankSendIniB, rankRecvIniB);
188 #define procRatioRows0 MATRIX_ROWS/procs[0]
189 #define procRatioRows1 MATRIX_ROWS/procs[1]
190 #define procRatioColumns0 MATRIX_COLUMNS/procs[0]
191 #define procRatioColumns1 MATRIX_COLUMNS/procs[1]
201 fprintf(stderr,
"CTRL %d, myCoords %d,%d\n",
mpiRank, myCoords[0], myCoords[1]);
205 endA[0] = (int)floor( ((
double)myCoords[0]+((MATRIX_ROWS > procs[0]) ? 1 : 0))*
procRatioRows0 - ((MATRIX_ROWS > procs[0]) ? 1 : 0) );
207 endA[1] = (int)floor( ((
double)myCoords[1]+((MATRIX_COLUMNS > procs[1]) ? 1 : 0))*
procRatioColumns1 - ((MATRIX_COLUMNS > procs[1]) ? 1 : 0) );
210 endB[0] = (int)floor( ((
double)myCoords[0]+((MATRIX_COLUMNS > procs[0]) ? 1 : 0))*
procRatioColumns0 - ((MATRIX_COLUMNS > procs[0]) ? 1 : 0) );
212 endB[1] = (int)floor( ((
double)myCoords[1]+((MATRIX_ROWS > procs[1]) ? 1 : 0))*
procRatioRows1 - ((MATRIX_ROWS > procs[1]) ? 1 : 0) );
215 endC[0] = (int)floor( ((
double)myCoords[0]+((MATRIX_ROWS > procs[0]) ? 1 : 0))*
procRatioRows0 - ((MATRIX_ROWS > procs[0]) ? 1 : 0) );
217 endC[1] = (int)floor( ((
double)myCoords[1]+((MATRIX_ROWS > procs[1]) ? 1 : 0))*
procRatioRows1 - ((MATRIX_ROWS > procs[1]) ? 1 : 0) );
220 fprintf(stderr,
"CTRL %d, A[%d:%d][%d:%d]\n",
mpiRank, beginA[0], endA[0], beginA[1], endA[1]);
221 fprintf(stderr,
"CTRL %d, B[%d:%d][%d:%d]\n",
mpiRank, beginB[0], endB[0], beginB[1], endB[1]);
222 fprintf(stderr,
"CTRL %d, C[%d:%d][%d:%d]\n",
mpiRank, beginC[0], endC[0], beginC[1], endC[1]);
226 MPI_Datatype dataBlockA;
227 MPI_Datatype dataBlockB;
228 MPI_Datatype dataBlockC;
230 ok = MPI_Type_hvector( endA[0]-beginA[0]+1, endA[1]-beginA[1]+1,
sizeof(
double)*(
size_t)MATRIX_COLUMNS, MPI_DOUBLE, &dataBlockA );
232 ok = MPI_Type_commit( &dataBlockA );
235 ok = MPI_Type_hvector( endB[0]-beginB[0]+1, endB[1]-beginB[1]+1,
sizeof(
double)*(
size_t)MATRIX_ROWS, MPI_DOUBLE, &dataBlockB );
237 ok = MPI_Type_commit( &dataBlockB );
240 ok = MPI_Type_hvector( endC[0]-beginC[0]+1, endC[1]-beginC[1]+1,
sizeof(
double)*(
size_t)MATRIX_ROWS, MPI_DOUBLE, &dataBlockC );
242 ok = MPI_Type_commit( &dataBlockC );
245 MPI_Request request[2];
251 for (i=beginA[0]; i<=endA[0]; i++) {
252 for (j=beginA[1]; j<=endA[1]; j++) {
253 fprintf(stderr,
"CTRL %d Before shift - A[%d][%d]=%lf\n",
mpiRank, i,j,
elemA(i,j) );
256 for (i=beginB[0]; i<=endB[0]; i++) {
257 for (j=beginB[1]; j<=endB[1]; j++) {
258 fprintf(stderr,
"CTRL %d Before shift - B[%d][%d]=%lf\n",
mpiRank, i,j,
elemB(i,j) );
265 ok = MPI_Isend( &(
elemA(beginA[0],beginA[1])), 1, dataBlockA, rankSendIniA, 0, MPI_COMM_WORLD, &request[0] );
268 ok = MPI_Isend( &(
elemB(beginB[0],beginB[1])), 1, dataBlockB, rankSendIniB, 1, MPI_COMM_WORLD, &request[1] );
271 ok = MPI_Recv( &(
elemCopyA(beginA[0],beginA[1])), 1, dataBlockA, rankRecvIniA, 0, MPI_COMM_WORLD, &status );
274 ok = MPI_Recv( &(
elemCopyB(beginB[0],beginB[1])), 1, dataBlockB, rankRecvIniB, 1, MPI_COMM_WORLD, &status );
277 MPI_Wait( &request[0], &status );
278 MPI_Wait( &request[1], &status );
280 for (i=beginA[0]; i<=endA[0]; i++) {
281 for (j=beginA[1]; j<=endA[1]; j++) {
285 for (i=beginB[0]; i<=endB[0]; i++) {
286 for (j=beginB[1]; j<=endB[1]; j++) {
295 for (i=beginA[0]; i<=endA[0]; i++) {
296 for (j=beginA[1]; j<=endA[1]; j++) {
297 fprintf(stderr,
"CTRL %d After shift - A[%d][%d]=%lf\n",
mpiRank, i,j,
elemA(i,j) );
300 for (i=beginB[0]; i<=endB[0]; i++) {
301 for (j=beginB[1]; j<=endB[1]; j++) {
302 fprintf(stderr,
"CTRL %d After shift - B[%d][%d]=%lf\n",
mpiRank, i,j,
elemB(i,j) );
311 for (loopIndex = 0; loopIndex < procs[0]-1; loopIndex++) {
314 for (i=beginC[0]; i<=endC[0]; i++) {
315 for (j=beginC[1]; j<=endC[1]; j++) {
316 for (k=0; k<=endA[1]-beginA[1]; k++) {
318 fprintf(stderr,
"CTRL %d, i,j,k:(%d,%d,%d) %lf x %lf Stage 1\n",
mpiRank, i,j,k,
elemA(i,k+begin[1]),
elemB(k+beginB[0],j));
326 ok = MPI_Isend( &
elemA(beginA[0],beginA[1]), 1, dataBlockA, rankLeft, 0, MPI_COMM_WORLD, &request[0] );
329 ok = MPI_Isend( &
elemB(beginB[0],beginB[1]), 1, dataBlockB, rankUp, 1, MPI_COMM_WORLD, &request[1] );
332 ok = MPI_Recv( &
elemCopyA(beginA[0],beginA[1]), 1, dataBlockA, rankRight, 0, MPI_COMM_WORLD, &status );
335 ok = MPI_Recv( &
elemCopyB(beginB[0],beginB[1]), 1, dataBlockB, rankDown, 1, MPI_COMM_WORLD, &status );
339 for (i=beginA[0]; i<=endA[0]; i++) {
340 for (j=beginA[1]; j<=endA[1]; j++) {
341 fprintf(stderr,
"CTRL %d After repeated comm. A[%d][%d] %lf\n",
mpiRank, i,j,
elemA(i,j));
346 MPI_Wait( &request[0], &status );
347 MPI_Wait( &request[1], &status );
349 for (i=beginA[0]; i<=endA[0]; i++) {
350 for (j=beginA[1]; j<=endA[1]; j++) {
354 for (i=beginB[0]; i<=endB[0]; i++) {
355 for (j=beginB[1]; j<=endB[1]; j++) {
365 ok = MPI_Type_free( &dataBlockA );
367 ok = MPI_Type_free( &dataBlockB );
371 for (i=beginC[0]; i<=endC[0]; i++) {
372 for (j=beginC[1]; j<=endC[1]; j++) {
373 for (k=0; k<=endA[1]-beginA[1]; k++) {
375 fprintf(stderr,
"CTRL %d, i,j,k:(%d,%d,%d) %lf x %lf Stage 1\n",
mpiRank, i,j,k,
elemA(i,k+begin[1]),
elemB(k+beginB[0],j));
383 for (i=beginC[0]; i<=endC[0]; i++) {
384 for (j=beginC[1]; j<=endC[1]; j++) {
385 fprintf(stderr,
"CTRL %d Final C[%d][%d] %lf\n",
mpiRank, i,j,
elemC(i,j));
394 for ( rank=1; rank<
mpiNProc; rank++ ) {
397 ok = MPI_Cart_coords(cartesianComm, rank, 2, remoteCoords);
398 mpiTestError( ok,
"Getting remote cartesian coordinates in end comm." );
406 remoteEndC[0] = (int)floor( (remoteCoords[0]+((MATRIX_ROWS > procs[0]) ? 1 : 0))*procRatioRows0 - ((MATRIX_ROWS > procs[0]) ? 1 : 0) );
407 remoteEndC[1] = (int)floor( (remoteCoords[1]+((MATRIX_ROWS > procs[1]) ? 1 : 0))*procRatioRows1 - ((MATRIX_ROWS > procs[1]) ? 1 : 0) );
409 fprintf(stderr,
"CTRL RemoteC %d, [%d:%d][%d:%d]\n", rank, remoteBeginC[0], remoteEndC[0], remoteBeginC[1], remoteEndC[1]);
412 ok = MPI_Recv( &(
elemC(remoteBeginC[0],remoteBeginC[1])), 1, dataBlockC, rank, rank, MPI_COMM_WORLD, &status );
413 mpiTestError( ok,
"Receiving remote parts in end comm." );
419 fprintf(stderr,
"CTRL LocalC %d, [%d:%d][%d:%d]\n",
mpiRank, beginC[0], endC[0], beginC[1], endC[1]);
421 ok = MPI_Send( &(
elemC(beginC[0],beginC[1])), 1, dataBlockC, 0,
mpiRank, MPI_COMM_WORLD );
425 ok = MPI_Type_free( &dataBlockC );
437 if ((f=fopen(
"C.dtxt",
"w")) == NULL) {
438 fprintf(stderr,
"Error: Imposible to open output file\n");
443 fprintf(f,
"%014.4lf",
elemC(i,j));
#define mpiTestError(ok, cad)
int main(int argc, char *argv[])
#define procRatioColumns1
#define procRatioColumns0