diff --git a/PoissonSolver/2D-mpi-sq/src/main.c b/PoissonSolver/2D-mpi-sq/src/main.c index af09d50..7c3f8b8 100644 --- a/PoissonSolver/2D-mpi-sq/src/main.c +++ b/PoissonSolver/2D-mpi-sq/src/main.c @@ -29,17 +29,16 @@ int main(int argc, char** argv) printf("Usage: %s \n", argv[0]); exit(EXIT_SUCCESS); } + readParameter(¶ms, argv[1]); initSolver(&solver, ¶ms, 2); - // writeResult(&solver, "p-0.dat"); // startTime = getTimeStamp(); // solveRB(&solver); // endTime = getTimeStamp(); + writeResult(&solver, "p.dat"); // printf(" %.2fs\n", endTime - startTime); // writeResult(&solver, "p.dat"); - if (rank == 0) { - } MPI_Finalize(); return EXIT_SUCCESS; diff --git a/PoissonSolver/2D-mpi-sq/src/solver.c b/PoissonSolver/2D-mpi-sq/src/solver.c index 3346091..0f979e1 100644 --- a/PoissonSolver/2D-mpi-sq/src/solver.c +++ b/PoissonSolver/2D-mpi-sq/src/solver.c @@ -9,11 +9,13 @@ #include "allocate.h" #include "parameter.h" #include "solver.h" +#include #include +#include #define PI 3.14159265358979323846 -#define P(i, j) p[(j) * (solver->limax + 2) + (i)] -#define RHS(i, j) rhs[(j) * (solver->limax + 2) + (i)] +#define P(i, j) p[(i) * (solver->ljmax + 2) + (j)] +#define RHS(i, j) rhs[(i) * (solver->ljmax + 2) + (j)] int distribute_dim(int dim_rank, int dim_comm_size, int dim_size) { @@ -48,6 +50,8 @@ void initSolver(Solver* solver, Parameter* params, int problem) solver->dims[1] = 0; MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &size); + solver->rank = rank; + solver->size = size; MPI_Dims_create(size, 2, solver->dims); MPI_Cart_create(MPI_COMM_WORLD, 2, @@ -59,51 +63,41 @@ 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); - solver->snd_displacements[0] = ((solver->limax + 2) + 1) * sizeof(double); - solver->snd_displacements[1] = ((solver->limax + 2) * (solver->ljmax) + 1) * + solver->snd_displacements[0] = ((solver->ljmax + 2) + 1) * sizeof(double); + solver->snd_displacements[1] = ((solver->ljmax + 2) * (solver->limax) + 1) * sizeof(double); - solver->snd_displacements[2] = ((solver->limax + 2) + 1) * sizeof(double); - solver->snd_displacements[3] = ((solver->limax + 2) + solver->limax) * sizeof(double); + solver->snd_displacements[2] = ((solver->ljmax + 2) + 1) * sizeof(double); + solver->snd_displacements[3] = ((solver->ljmax + 2) + solver->ljmax) * sizeof(double); solver->rcv_displacements[0] = (1) * sizeof(double); - solver->rcv_displacements[1] = ((solver->limax + 2) * (solver->ljmax + 1) + 1) * + solver->rcv_displacements[1] = ((solver->ljmax + 2) * (solver->limax + 1) + 1) * sizeof(double); - solver->rcv_displacements[2] = ((solver->limax + 2)) * sizeof(double); - solver->rcv_displacements[3] = ((solver->limax + 2) + solver->limax + 1) * + solver->rcv_displacements[2] = ((solver->ljmax + 2)) * sizeof(double); + solver->rcv_displacements[3] = ((solver->ljmax + 2) + solver->ljmax + 1) * sizeof(double); MPI_Datatype contiguos_boundary; MPI_Datatype vector_boundary; - MPI_Datatype internal; - - MPI_Type_contiguous(solver->limax, MPI_DOUBLE, &contiguos_boundary); - MPI_Type_vector(solver->ljmax, 1, solver->limax + 2, MPI_DOUBLE, &vector_boundary); - MPI_Type_create_subarray(2, - (int[2]){solver->limax+2,solver->ljmax+2}, - (int[2]){solver->limax,solver->ljmax}, - (int[2]){1,1}, - MPI_ORDER_C, - MPI_DOUBLE, - &internal); - MPI_Type_commit(&internal); + 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->types[4] = internal; size_t bytesize = (solver->limax + 2) * (solver->ljmax + 2) * sizeof(double); - solver->p = allocate(64, bytesize); - solver->rhs = allocate(64, bytesize); + solver->p = malloc(bytesize); + solver->rhs = malloc(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); + printf("Prc %d; {%d,%d}, {%d,%d}, {%d,%d}\n", rank, solver->coords[0], @@ -113,23 +107,27 @@ void initSolver(Solver* solver, Parameter* params, int problem) i_start, j_start); fflush(stdout); - - for (int j = 0; j < solver->ljmax + 2; j++) { - for (int i = 0; i < solver->limax + 2; i++) { + for (int i = 0; i < solver->limax + 2; i++) { + for (int j = 0; j < solver->ljmax + 2; j++) { + P(i, j) = rank; + } + } + 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); } } if (problem == 2) { - for (int j = 0; j < solver->ljmax + 2; j++) { - for (int i = 0; i < solver->limax + 2; i++) { + 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); } } } else { - for (int j = 0; j < solver->ljmax + 2; j++) { - for (int i = 0; i < solver->limax + 2; i++) { + for (int i = 0; i < solver->limax + 2; i++) { + for (int j = 0; j < solver->ljmax + 2; j++) { RHS(i, j) = 0.0; } } @@ -153,16 +151,24 @@ 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; - + res = 0.0; + jsw = 1; + int count[] = { 1, 1, 1, 1 }; + MPI_Neighbor_alltoallw(solver->p, + count, + solver->snd_displacements, + solver->types, + solver->p, + count, + solver->rcv_displacements, + solver->types, + solver->cart_comm); for (pass = 0; pass < 2; pass++) { isw = jsw; - for (int j = 1; j < jmax + 1; j++) { - for (int i = isw; i < imax + 1; i += 2) { + for (int j = 1; j < solver->ljmax + 1; j++) { + for (int i = isw; i < solver->limax + 1; i += 2) { double r = RHS(i, j) - ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 + @@ -175,16 +181,16 @@ void solveRB(Solver* solver) } jsw = 3 - jsw; } - - for (int i = 1; i < imax + 1; i++) { - P(i, 0) = P(i, 1); - P(i, jmax + 1) = P(i, jmax); - } - - for (int j = 1; j < jmax + 1; j++) { - P(0, j) = P(1, j); - P(imax + 1, j) = P(imax, j); - } + // + // for (int i = 1; i < imax + 1; i++) { + // P(i, 0) = P(i, 1); + // P(i, jmax + 1) = P(i, jmax); + // } + // + // for (int j = 1; j < jmax + 1; j++) { + // P(0, j) = P(1, j); + // P(imax + 1, j) = P(imax, j); + // } res = res / (double)(imax * jmax); #ifdef DEBUG @@ -199,24 +205,89 @@ void solveRB(Solver* solver) void writeResult(Solver* solver, char* filename) { - int imax = solver->imax; - int jmax = solver->jmax; - double* p = solver->p; + int imax = solver->imax; + int jmax = solver->jmax; + double* p = solver->p; + int* rcv_count = NULL; + int* displacements = NULL; + double* p_total = NULL; + 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; + 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] = + solver->p[(i + 1) * (solver->ljmax + 2) + (j + 1)]; - FILE* fp; - fp = fopen(filename, "w"); - - if (fp == NULL) { - printf("Error!\n"); - exit(EXIT_FAILURE); - } - - for (int j = 0; j < jmax + 2; j++) { - for (int i = 0; i < imax + 2; i++) { - fprintf(fp, "%f ", P(i, j)); + 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, + p_total, + rcv_count, + displacements, + 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]; } - fprintf(fp, "\n"); - } + FILE* fp; + fp = fopen(filename, "w"); - fclose(fp); + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + 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)); + 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, + p_total, + rcv_count, + displacements, + MPI_DOUBLE, + 0, + MPI_COMM_WORLD); + } } diff --git a/PoissonSolver/2D-mpi-sq/src/solver.h b/PoissonSolver/2D-mpi-sq/src/solver.h index 01961b3..42d3eac 100644 --- a/PoissonSolver/2D-mpi-sq/src/solver.h +++ b/PoissonSolver/2D-mpi-sq/src/solver.h @@ -15,11 +15,12 @@ typedef struct { double eps, omega, rho; int itermax; MPI_Comm cart_comm; - MPI_Datatype types[5]; + MPI_Datatype types[4]; MPI_Aint rcv_displacements[4]; MPI_Aint snd_displacements[4]; int coords[2]; int dims[2]; + int rank, size; } Solver; extern void initSolver(Solver*, Parameter*, int problem);