From 2487fe8b9b110435e7dbccaba42962883cde2179 Mon Sep 17 00:00:00 2001
From: David Goz <david.goz@inaf.it>
Date: Fri, 21 Jun 2024 14:44:00 +0200
Subject: [PATCH] mpi/omp_comm round-off error spotting
---
jacobi/mpi/comp_comm/include/allvars.h | 2 +-
jacobi/mpi/comp_comm/make_mpi_path | 2 +-
.../comp_comm/src/jacobi_2D_mpi_comp_comm.c | 140 ++++++++++--------
3 files changed, 81 insertions(+), 63 deletions(-)
diff --git a/jacobi/mpi/comp_comm/include/allvars.h b/jacobi/mpi/comp_comm/include/allvars.h
index f6f93c1..c065228 100644
--- a/jacobi/mpi/comp_comm/include/allvars.h
+++ b/jacobi/mpi/comp_comm/include/allvars.h
@@ -16,5 +16,5 @@ typedef float MyData;
#else
/* double precision is assumed by default */
typedef double MyData;
-#define MPI_MyDatatype MPI_DOUBLE_PRECISION
+#define MPI_MyDatatype MPI_DOUBLE_PRECISION
#endif /* SINGLE_PRECISION */
diff --git a/jacobi/mpi/comp_comm/make_mpi_path b/jacobi/mpi/comp_comm/make_mpi_path
index f2c3de9..75146ce 100644
--- a/jacobi/mpi/comp_comm/make_mpi_path
+++ b/jacobi/mpi/comp_comm/make_mpi_path
@@ -1,6 +1,6 @@
# set the MPI install path
-MPI_INSTALL_PATH = /home/gozzilla/software/openmpi/openmpi-5.0.3
+MPI_INSTALL_PATH = /home/darkenergy/software/openmpi/openmpi-5.0.3
diff --git a/jacobi/mpi/comp_comm/src/jacobi_2D_mpi_comp_comm.c b/jacobi/mpi/comp_comm/src/jacobi_2D_mpi_comp_comm.c
index 0a9dc72..b63dca1 100644
--- a/jacobi/mpi/comp_comm/src/jacobi_2D_mpi_comp_comm.c
+++ b/jacobi/mpi/comp_comm/src/jacobi_2D_mpi_comp_comm.c
@@ -51,7 +51,7 @@ void JacobiAlgorithm(MyData **const restrict, MyData **const restrict, const MyD
void Jacobi_Communication(MyData **const restrict Phi,
MyData **const restrict Phi0,
- MyData *const restrict error,
+ MyData *const restrict error,
const MyData * restrict delta,
Task *const restrict ThisTask);
@@ -217,22 +217,8 @@ int main(int argc, char **argv)
1. Set grid indices
-------------------------------------------------------- */
- const int ibeg = NGHOST;
- /* const int iend = ibeg + NX_GLOB - 1; */
- /* const int nx = iend - ibeg + 1; */
- /* const int nx_tot = nx + 2 * NGHOST; */
-
- const int jbeg = NGHOST;
- /* const int jend = jbeg + NY_GLOB - 1; */
- /* const int ny = jend - jbeg + 1; */
- /* const int ny_tot = ny + 2 * NGHOST; */
-
- /* if (rank == MASTER) */
- /* { */
- /* printf("\n\t Grid indices:"); */
- /* printf("\n\t\t ibeg, iend = %d, %d; nx_tot = %d" ,ibeg, iend, nx_tot); */
- /* printf("\n\t\t jbeg, jend = %d, %d; ny_tot = %d\n\n",jbeg, jend, ny_tot); */
- /* } */
+ const int ibeg = NGHOST;
+ const int jbeg = NGHOST;
/* --------------------------------------------------------
2. Generate grids, allocate memory
@@ -425,14 +411,35 @@ void JacobiAlgorithm(MyData **const restrict Phi,
{
Phi[j][i] = 0.25 * (Phi0[j][i-1] + Phi0[j][i+1] +
Phi0[j-1][i] + Phi0[j+1][i]);
-
- *error += delta[X] * delta[Y] * fabs(Phi[j][i] - Phi0[j][i]);
+
+ *error += (delta[X] * delta[Y] * fabs(Phi[j][i] - Phi0[j][i]));
} /* loop over columns */
} /* loop over rows */
return;
}
+void get_error(MyData **const restrict Phi,
+ MyData **const restrict Phi0,
+ const MyData *restrict delta,
+ const int jbeg,
+ const int jend,
+ const int ibeg,
+ const int iend,
+ MyData *const restrict error)
+{
+ *error = 0.0;
+ for (int j=jbeg ; j<=jend ; j++)
+ {
+ for (int i=ibeg ; i<=iend ; i++)
+ {
+ *error += (delta[X] * delta[Y] * fabs(Phi[j][i] - Phi0[j][i]));
+ } /* loop over columns */
+ } /* loop over rows */
+
+ return;
+}
+
void Jacobi_Communication(MyData **const restrict Phi,
MyData **const restrict Phi0,
MyData *const restrict error,
@@ -448,21 +455,25 @@ void Jacobi_Communication(MyData **const restrict Phi,
const int data_row_size = ThisTask->domain.dim[Y];
/* First task: issue the communication */
- MPI_Request request[8];
+ MPI_Request request[4];
MyData **const restrict buffer = Phi0;
- /* 4 nonblocking MPI_Irecv */
- MPI_Irecv(&buffer[ThisTask->domain.local_start[X] - 1][ThisTask->domain.local_start[Y] ], data_row_size, MPI_MyDatatype, ThisTask->nbrbottom, 0, ThisTask->comm2d, &request[0]);
- MPI_Irecv(&buffer[ThisTask->domain.local_end[X] + 1][ThisTask->domain.local_start[Y] ], data_row_size, MPI_MyDatatype, ThisTask->nbrtop, 1, ThisTask->comm2d, &request[1]);
- MPI_Irecv(&buffer[ThisTask->domain.local_start[X] ][ThisTask->domain.local_start[Y] - 1], 1, column, ThisTask->nbrleft, 2, ThisTask->comm2d, &request[2]);
- MPI_Irecv(&buffer[ThisTask->domain.local_start[X] ][ThisTask->domain.local_end[Y] + 1], 1, column, ThisTask->nbrright, 3, ThisTask->comm2d, &request[3]);
-
- /* 4 nonblocking MPI_Isend */
- MPI_Isend(&buffer[ThisTask->domain.local_end[X] ][ThisTask->domain.local_start[Y]], data_row_size, MPI_MyDatatype, ThisTask->nbrtop, 0, ThisTask->comm2d, &request[4]);
- MPI_Isend(&buffer[ThisTask->domain.local_start[X]][ThisTask->domain.local_start[Y]], data_row_size, MPI_MyDatatype, ThisTask->nbrbottom, 1, ThisTask->comm2d, &request[5]);
- MPI_Isend(&buffer[ThisTask->domain.local_start[X]][ThisTask->domain.local_end[Y] ], 1, column, ThisTask->nbrright, 2, ThisTask->comm2d, &request[6]);
- MPI_Isend(&buffer[ThisTask->domain.local_start[X]][ThisTask->domain.local_start[Y]], 1, column, ThisTask->nbrleft, 3, ThisTask->comm2d, &request[7]);
+ MPI_Isendrecv(&buffer[ThisTask->domain.local_end[X] ][ThisTask->domain.local_start[Y] ], data_row_size, MPI_MyDatatype, ThisTask->nbrtop, 0,
+ &buffer[ThisTask->domain.local_start[X] - 1][ThisTask->domain.local_start[Y] ], data_row_size, MPI_MyDatatype, ThisTask->nbrbottom, 0,
+ ThisTask->comm2d, &request[0]);
+
+ MPI_Isendrecv(&buffer[ThisTask->domain.local_start[X] ][ThisTask->domain.local_start[Y] ], data_row_size, MPI_MyDatatype, ThisTask->nbrbottom, 1,
+ &buffer[ThisTask->domain.local_end[X] + 1][ThisTask->domain.local_start[Y] ], data_row_size, MPI_MyDatatype, ThisTask->nbrtop, 1,
+ ThisTask->comm2d, &request[1]);
+
+ MPI_Isendrecv(&buffer[ThisTask->domain.local_start[X] ][ThisTask->domain.local_end[Y] ], 1, column, ThisTask->nbrright, 2,
+ &buffer[ThisTask->domain.local_start[X] ][ThisTask->domain.local_start[Y] - 1], 1, column, ThisTask->nbrleft, 2,
+ ThisTask->comm2d, &request[2]);
+
+ MPI_Isendrecv(&buffer[ThisTask->domain.local_start[X] ][ThisTask->domain.local_start[Y] ], 1, column, ThisTask->nbrleft, 3,
+ &buffer[ThisTask->domain.local_start[X] ][ThisTask->domain.local_end[Y] + 1], 1, column, ThisTask->nbrright, 3,
+ ThisTask->comm2d, &request[3]);
/**************************************** computation ****************************************/
/* perform the computation with the local data, (i.e. ghost cells are not required) */
@@ -472,40 +483,47 @@ void Jacobi_Communication(MyData **const restrict Phi,
const int jend = ThisTask->domain.local_end[X] - 1;
const int ibeg = ThisTask->domain.local_start[Y] + 1;
const int iend = ThisTask->domain.local_end[Y] - 1;
-
- /* wait the data on the boundaries */
- MPI_Waitall(8, request, MPI_STATUSES_IGNORE);
*error = 0.0;
- /* JacobiAlgorithm(Phi, Phi0, delta, jbeg, jend, ibeg, iend, error); */
- JacobiAlgorithm(Phi, Phi0, delta, jbeg-1, jend+1, ibeg-1, iend+1, error);
+ JacobiAlgorithm(Phi, Phi0, delta, jbeg, jend, ibeg, iend, error);
- /* /\* wait the data on the boundaries *\/ */
- /* MPI_Waitall(8, request, MPI_STATUSES_IGNORE); */
+ /*********************************************************************************************/
+
+ /* wait the data on the boundaries */
+ MPI_Waitall(4, request, MPI_STATUSES_IGNORE);
- /* /\* nbrbottom *\/ */
- /* JacobiAlgorithm(Phi, Phi0, delta, */
- /* ThisTask->domain.local_start[X], ThisTask->domain.local_start[X], */
- /* ThisTask->domain.local_start[Y], ThisTask->domain.local_end[Y], */
- /* error); */
-
- /* /\* nbrtop *\/ */
- /* JacobiAlgorithm(Phi, Phi0, delta, */
- /* ThisTask->domain.local_end[X], ThisTask->domain.local_end[X], */
- /* ThisTask->domain.local_start[Y], ThisTask->domain.local_end[Y], */
- /* error); */
-
- /* /\* nbrleft *\/ */
- /* JacobiAlgorithm(Phi, Phi0, delta, */
- /* ThisTask->domain.local_start[X], ThisTask->domain.local_end[X], */
- /* ThisTask->domain.local_start[Y], ThisTask->domain.local_start[Y], */
- /* error); */
-
- /* /\* nbrright *\/ */
- /* JacobiAlgorithm(Phi, Phi0, delta, */
- /* ThisTask->domain.local_start[X], ThisTask->domain.local_end[X], */
- /* ThisTask->domain.local_end[Y], ThisTask->domain.local_end[Y], */
- /* error); */
+ /* nbrbottom */
+ JacobiAlgorithm(Phi, Phi0, delta,
+ ThisTask->domain.local_start[X], ThisTask->domain.local_start[X],
+ ThisTask->domain.local_start[Y], ThisTask->domain.local_end[Y],
+ error);
+
+ /* nbrtop */
+ JacobiAlgorithm(Phi, Phi0, delta,
+ ThisTask->domain.local_end[X], ThisTask->domain.local_end[X],
+ ThisTask->domain.local_start[Y], ThisTask->domain.local_end[Y],
+ error);
+
+ /* nbrleft */
+ JacobiAlgorithm(Phi, Phi0, delta,
+ ThisTask->domain.local_start[X], ThisTask->domain.local_end[X],
+ ThisTask->domain.local_start[Y], ThisTask->domain.local_start[Y],
+ error);
+
+ /* nbrright */
+ JacobiAlgorithm(Phi, Phi0, delta,
+ ThisTask->domain.local_start[X], ThisTask->domain.local_end[X],
+ ThisTask->domain.local_end[Y], ThisTask->domain.local_end[Y],
+ error);
+
+ /* Round-off error fixing ??? */
+ {
+ *error = 0.0;
+ get_error(Phi, Phi0, delta,
+ ThisTask->domain.local_start[X], ThisTask->domain.local_end[X],
+ ThisTask->domain.local_start[Y], ThisTask->domain.local_end[Y],
+ error);
+ }
MPI_Type_free(&column);
--
GitLab