diff --git a/PoissonSolver/2D-mpi-sq/src/main.c b/PoissonSolver/2D-mpi-sq/src/main.c index 89e92d1..7d73456 100644 --- a/PoissonSolver/2D-mpi-sq/src/main.c +++ b/PoissonSolver/2D-mpi-sq/src/main.c @@ -14,15 +14,11 @@ int main(int argc, char** argv) { MPI_Init(&argc, &argv); - /******* Debug Params *********/ - int rank, size; - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - MPI_Comm_size(MPI_COMM_WORLD, &size); - printf("Process %d of %d \n", rank, size); - + double startTime, endTime; Parameter params; Solver solver; + initParameter(¶ms); if (argc < 2) { @@ -36,9 +32,7 @@ int main(int argc, char** argv) solveRB(&solver); endTime = getTimeStamp(); writeResult(&solver, "p.dat"); - // printf(" %.2fs\n", endTime - startTime); - // writeResult(&solver, "p.dat"); - + if(solver.rank == 0) printf(" %.2fs\n", endTime - startTime); MPI_Finalize(); return EXIT_SUCCESS; diff --git a/PoissonSolver/2D-mpi-sq/src/solver.c b/PoissonSolver/2D-mpi-sq/src/solver.c index d2091be..6562cb2 100644 --- a/PoissonSolver/2D-mpi-sq/src/solver.c +++ b/PoissonSolver/2D-mpi-sq/src/solver.c @@ -11,6 +11,7 @@ #include "solver.h" #include #include +#include #include #define PI 3.14159265358979323846 @@ -31,8 +32,36 @@ int get_dim_start(int dim_rank, int dim_comm_size, int dim_size) return dim_start; } +void write_matrix_to_file( + const char* _Nonnull filename, double* matrix, size_t rows, size_t cols) +{ + FILE* fp; + fp = fopen(filename, "w"); + + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + fprintf(fp, "%f ", matrix[i * (cols) + j]); + } + fprintf(fp, "\n"); + } + + fclose(fp); +} + void initSolver(Solver* solver, Parameter* params, int problem) { + int x_low, x_high, y_low, y_high; + int rank, size; + MPI_Datatype contiguos_boundary, vector_boundary; + + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); + solver->imax = params->imax; solver->jmax = params->jmax; solver->dx = params->xlength / params->imax; @@ -41,17 +70,13 @@ void initSolver(Solver* solver, Parameter* params, int problem) solver->omega = params->omg; solver->rho = params->rho; solver->itermax = params->itermax; - - int imax = solver->imax; - int jmax = solver->jmax; - - int rank, size; solver->dims[0] = 0; solver->dims[1] = 0; - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - MPI_Comm_size(MPI_COMM_WORLD, &size); - solver->rank = rank; - solver->size = size; + solver->rank = 0; + solver->size = 0; + solver->rank = rank; + solver->size = size; + MPI_Dims_create(size, 2, solver->dims); MPI_Cart_create(MPI_COMM_WORLD, 2, @@ -61,9 +86,9 @@ void initSolver(Solver* solver, Parameter* params, int problem) &solver->cart_comm); MPI_Cart_coords(solver->cart_comm, rank, 2, solver->coords); solver->boundary = NB; - int x_low, x_high, y_low, y_high; MPI_Cart_shift(solver->cart_comm, 0, 1, &x_low, &x_high); MPI_Cart_shift(solver->cart_comm, 1, 1, &y_low, &y_high); + if (x_low == MPI_PROC_NULL) solver->boundary += BOTTOM; if (y_low == MPI_PROC_NULL) solver->boundary += LEFT; if (x_high == MPI_PROC_NULL) solver->boundary += TOP; @@ -72,6 +97,15 @@ void initSolver(Solver* solver, Parameter* params, int problem) solver->limax = distribute_dim(solver->coords[0], solver->dims[0], solver->imax); solver->ljmax = distribute_dim(solver->coords[1], solver->dims[1], solver->jmax); + MPI_Type_contiguous(solver->ljmax, MPI_DOUBLE, &contiguos_boundary); + MPI_Type_vector(solver->limax, 1, solver->ljmax + 2, MPI_DOUBLE, &vector_boundary); + MPI_Type_commit(&contiguos_boundary); + MPI_Type_commit(&vector_boundary); + solver->types[0] = contiguos_boundary; + solver->types[1] = contiguos_boundary; + solver->types[2] = vector_boundary; + solver->types[3] = vector_boundary; + solver->snd_displacements[0] = ((solver->ljmax + 2) + 1) * sizeof(double); solver->snd_displacements[1] = ((solver->ljmax + 2) * (solver->limax) + 1) * sizeof(double); @@ -85,53 +119,39 @@ void initSolver(Solver* solver, Parameter* params, int problem) solver->rcv_displacements[3] = ((solver->ljmax + 2) + solver->ljmax + 1) * sizeof(double); - MPI_Datatype contiguos_boundary; - MPI_Datatype vector_boundary; - - MPI_Type_contiguous(solver->ljmax, MPI_DOUBLE, &contiguos_boundary); - MPI_Type_vector(solver->limax, 1, solver->ljmax + 2, MPI_DOUBLE, &vector_boundary); - MPI_Type_commit(&contiguos_boundary); - MPI_Type_commit(&vector_boundary); - solver->types[0] = contiguos_boundary; - solver->types[1] = contiguos_boundary; - solver->types[2] = vector_boundary; - solver->types[3] = vector_boundary; - size_t bytesize = (solver->limax + 2) * (solver->ljmax + 2) * sizeof(double); - solver->p = malloc(bytesize); - solver->rhs = malloc(bytesize); + solver->p = allocate(64, bytesize); + solver->rhs = allocate(64, bytesize); double dx = solver->dx; double dy = solver->dy; double* p = solver->p; double* rhs = solver->rhs; - int i_start = get_dim_start(solver->coords[0], solver->dims[0], solver->imax); - int j_start = get_dim_start(solver->coords[1], solver->dims[1], solver->jmax); + solver->i_start = get_dim_start(solver->coords[0], solver->dims[0], solver->imax); + solver->j_start = get_dim_start(solver->coords[1], solver->dims[1], solver->jmax); +#ifdef DEBUG printf("Prc %d; {%d,%d}, {%d,%d}, {%d,%d}\n", rank, solver->coords[0], solver->coords[1], solver->limax, solver->ljmax, - i_start, - j_start); + solver->i_start, + solver->j_start); fflush(stdout); - // for (int i = 0; i < solver->limax + 2; i++) { - // for (int j = 0; j < solver->ljmax + 2; j++) { - // P(i, j) = (double)solver->boundary; - // } - // } +#endif /* ifdef DEBUG */ + for (int i = 0; i < solver->limax + 2; i++) { for (int j = 0; j < solver->ljmax + 2; j++) { - P(i, j) = sin(2.0 * PI * (i_start + i) * dx * 2.0) + - sin(2.0 * PI * (j_start + j) * dy * 2.0); + P(i, j) = sin(2.0 * PI * (solver->i_start + i) * dx * 2.0) + + sin(2.0 * PI * (solver->j_start + j) * dy * 2.0); } } if (problem == 2) { for (int i = 0; i < solver->limax + 2; i++) { for (int j = 0; j < solver->ljmax + 2; j++) { - RHS(i, j) = sin(2.0 * PI * (i_start + i) * dx); + RHS(i, j) = sin(2.0 * PI * (solver->i_start + i) * dx); } } } else { @@ -160,25 +180,25 @@ void solveRB(Solver* solver) int it = 0; double res = 1.0; int pass, jsw, isw; + while ((res >= epssq) && (it < itermax)) { - res = 0.0; - jsw = 1; - int count[] = { 1, 1, 1, 1 }; + res = 0.0; + jsw = solver->j_start % 2 ? 1 : 2; + MPI_Neighbor_alltoallw(solver->p, - count, + (int[4]) { 1, 1, 1, 1 }, solver->snd_displacements, solver->types, solver->p, - count, + (int[4]) { 1, 1, 1, 1 }, solver->rcv_displacements, solver->types, solver->cart_comm); - for (pass = 0; pass < 2; pass++) { + + for (pass = 0; pass < 2; pass++) { // why ?? isw = jsw; - - for (int j = 1; j < solver->ljmax + 1; j++) { - for (int i = isw; i < solver->limax + 1; i += 2) { - + for (int i = 1; i < solver->limax + 1; i++) { + for (int j = isw; j < solver->ljmax + 1; j += 2) { double r = RHS(i, j) - ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 + (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2); @@ -190,7 +210,9 @@ void solveRB(Solver* solver) } jsw = 3 - jsw; } + MPI_Allreduce(MPI_IN_PLACE, &res, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + if (solver->boundary & BOTTOM) for (int j = 1; j < solver->ljmax + 1; j++) P(0, j) = P(1, j); @@ -205,31 +227,40 @@ void solveRB(Solver* solver) P(i, solver->ljmax + 1) = P(i, solver->ljmax); res = res / (double)(imax * jmax); + #ifdef DEBUG printf("%d Residuum: %e\n", it, res); #endif it++; } - if(solver->rank==0) printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); + if (solver->rank == 0) + printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); } void writeResult(Solver* solver, char* filename) { - int imax = solver->imax; - int jmax = solver->jmax; - double* p = solver->p; - int* rcv_count = NULL; - int* displacements = NULL; - double* p_total = NULL; + int imax = solver->imax; + int jmax = solver->jmax; + double* p = solver->p; + int* rcv_count = NULL; + int* displacements = NULL; + double* p_total = NULL; + double* send_buffer = NULL; + int local_size = solver->limax * solver->ljmax; + if (solver->rank == 0) { - rcv_count = malloc(solver->size * sizeof(int)); - displacements = calloc(solver->size, sizeof(int)); - p_total = malloc(sizeof(double) * solver->imax * solver->jmax); - double* p_total_ordered = malloc(sizeof(double) * solver->imax * solver->jmax); - int local_size = solver->limax * solver->ljmax; + + rcv_count = allocate(64, solver->size * sizeof(int)); + displacements = allocate(64, solver->size * sizeof(int)); + p_total = allocate(64, sizeof(double) * solver->imax * solver->jmax); + double* p_total_ordered = allocate(64, + sizeof(double) * solver->imax * solver->jmax); + MPI_Gather(&local_size, 1, MPI_INT, rcv_count, 1, MPI_INT, 0, MPI_COMM_WORLD); + displacements[0] = 0; + for (int i = 0; i < solver->limax; i++) for (int j = 0; j < solver->ljmax; j++) p_total_ordered[i * (solver->ljmax) + j] = @@ -237,14 +268,7 @@ void writeResult(Solver* solver, char* filename) for (int i = 1; i < solver->size; i++) displacements[i] = displacements[i - 1] + rcv_count[i - 1]; - printf("{ "); - for (int i = 0; i < solver->size; i++) - printf("%d ", rcv_count[i]); - printf("}\n"); - printf("{ "); - for (int i = 0; i < solver->size; i++) - printf("%d ", displacements[i]); - printf("}\n"); + MPI_Gatherv(p_total_ordered, local_size, MPI_DOUBLE, @@ -254,43 +278,37 @@ void writeResult(Solver* solver, char* filename) MPI_DOUBLE, 0, MPI_COMM_WORLD); + for (int process = 0; process < solver->size; process++) { + int coords[2] = { 0 }; + MPI_Cart_coords(solver->cart_comm, process, 2, coords); + int i_max_local = distribute_dim(coords[0], solver->dims[0], solver->imax); int j_max_local = distribute_dim(coords[1], solver->dims[1], solver->jmax); int i_start_local = get_dim_start(coords[0], solver->dims[0], solver->imax); int j_start_local = get_dim_start(coords[1], solver->dims[1], solver->jmax); + for (int i = 0; i < i_max_local; i++) for (int j = 0; j < j_max_local; j++) p_total_ordered[(i + i_start_local) * (solver->jmax) + (j + j_start_local)] = p_total[displacements[process] + i * (j_max_local) + j]; } - FILE* fp; - fp = fopen(filename, "w"); - if (fp == NULL) { - printf("Error!\n"); - exit(EXIT_FAILURE); - } + write_matrix_to_file(filename, p_total_ordered, solver->imax, solver->jmax); - for (int j = 0; j < jmax; j++) { - for (int i = 0; i < imax; i++) { - fprintf(fp, "%f ", p_total_ordered[i * (solver->jmax) + j]); - } - fprintf(fp, "\n"); - } - - fclose(fp); } else { - int local_size = solver->limax * solver->ljmax; - double* send_buffer = malloc(local_size * sizeof(double)); + + send_buffer = allocate(64, local_size * sizeof(double)); for (int i = 0; i < solver->limax; i++) for (int j = 0; j < solver->ljmax; j++) send_buffer[i * (solver->ljmax) + j] = solver->p[(i + 1) * (solver->ljmax + 2) + (j + 1)]; + MPI_Gather(&local_size, 1, MPI_INT, rcv_count, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Gatherv(send_buffer, local_size, MPI_DOUBLE, diff --git a/PoissonSolver/2D-mpi-sq/src/solver.h b/PoissonSolver/2D-mpi-sq/src/solver.h index 0b87fa6..407fc18 100644 --- a/PoissonSolver/2D-mpi-sq/src/solver.h +++ b/PoissonSolver/2D-mpi-sq/src/solver.h @@ -23,6 +23,7 @@ typedef struct { double dx, dy; int imax, jmax; int limax, ljmax; + int i_start, j_start; double *p, *rhs; double eps, omega, rho; int itermax;