sequential.c
#include<stdio.h> #include<stdlib.h> #include<math.h> void advance_field(int nx, int ny, double *f, double *g, double *c) { int i, j, idx; for (i=1; i<nx-1; i++) { for (j=1; j<ny-1; j++) { idx = i*ny + j; f[idx] = c[idx]*(g[idx+ny] + g[idx-ny] + g[idx+1] + g[idx-1] - 4*g[idx]) + 2*g[idx] - f[idx]; } } } void periodic_y(int nx, int ny, double *f) { int i, j; for (i=0; i<nx; i++) { f[i*ny+(ny-1)] = f[i*ny+1]; f[i*ny+0] = f[i*ny+(ny-2)]; } } int main() { int nx=2000, ny=2000; int tmax=1500; int width=120, thick=60, gap=1200; // slit parameters int distance=800; double c0=0.49; double *f, *g, *c; int i, j, tstep; //------------------------------------------------------------------------- // initialize coefficient and fields //------------------------------------------------------------------------- f = (double*)malloc(nx*ny*sizeof(double)); g = (double*)malloc(nx*ny*sizeof(double)); c = (double*)malloc(nx*ny*sizeof(double)); for (i=0; i<nx*ny; i++) { f[i] = 0; g[i] = 0; c[i] = c0; } //------------------------------------------------------------------------- // slit structure //------------------------------------------------------------------------- for (i=nx/2-1; i<nx/2+thick; i++) { for (j=0; j<ny; j++) c[i*ny+j] = 0; } for (i=nx/2-1; i<nx/2+thick; i++) { for (j=ny/2-(gap+width)/2-1; j<ny/2+(gap+width)/2; j++) c[i*ny+j] = c0; } for (i=nx/2-1; i<nx/2+thick; i++) { for (j=ny/2-(gap-width)/2-1; j<ny/2+(gap-width)/2; j++) c[i*ny+j] = 0; } //------------------------------------------------------------------------- // main loop for the time evolution //------------------------------------------------------------------------- for (tstep=1; tstep<=tmax; tstep++) { //------------------------------------ // point source //------------------------------------ //g[nx/2-1][ny/2-1] = sin(0.02*tstep); //------------------------------------ // line source //------------------------------------ for (j=0; j<ny; j++) g[(nx/2-distance-1)*ny+j] = sin(0.02*tstep); advance_field(nx, ny, f, g, c); periodic_y(nx, ny, f); advance_field(nx, ny, g, f, c); periodic_y(nx, ny, g); } //------------------------------------------------------------------------- // gather fields and save as binary files //------------------------------------------------------------------------- FILE *fout; fout = fopen("field.bin", "wb"); fwrite(f, sizeof(double), nx*ny, fout); fclose(fout); return 0; }
answer.c
#include<stdio.h> #include<stdlib.h> #include<math.h> #include<mpi.h> void advance_field(int nx, int ny, double *f, double *g, double *c) { int i, j, idx; for (i=1; i<nx-1; i++) { for (j=1; j<ny-1; j++) { idx = i*ny + j; f[idx] = c[idx]*(g[idx+ny] + g[idx-ny] + g[idx+1] + g[idx-1] - 4*g[idx]) + 2*g[idx] - f[idx]; } } } void periodic_y(int nx, int ny, double *f) { int i, j; for (i=0; i<nx; i++) { f[i*ny+(ny-1)] = f[i*ny+1]; f[i*ny+0] = f[i*ny+(ny-2)]; } } void exchange_boundary(int nx, int ny, double *f) { int i, j; int nprocs, myrank; MPI_Request req1, req2, req3, req4; MPI_Status stat; MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &myrank); if (myrank < nprocs-1) { MPI_Isend(&f[(nx-2)*ny], ny, MPI_DOUBLE, myrank+1, 10, MPI_COMM_WORLD, &req1); MPI_Irecv(&f[(nx-1)*ny], ny, MPI_DOUBLE, myrank+1, 20, MPI_COMM_WORLD, &req2); } if (myrank > 0 ) { MPI_Isend(&f[ny], ny, MPI_DOUBLE, myrank-1, 20, MPI_COMM_WORLD, &req3); MPI_Irecv(&f[0], ny, MPI_DOUBLE, myrank-1, 10, MPI_COMM_WORLD, &req4); } if (myrank < nprocs-1) { MPI_Wait(&req1, &stat); MPI_Wait(&req2, &stat); } if (myrank > 0 ) { MPI_Wait(&req3, &stat); MPI_Wait(&req4, &stat); } } int main(int argc, char *argv[]) { int tnx=8000, tny=8000; int tmax=3000; int width=120, thick=60, gap=1200; // slit parameters int distance=1500; double c0=0.49; double *f, *g, *c, *ftot; int nx, ny; int i, j, idx, tstep; int nprocs, myrank; //------------------------------------------------------------------------- // initialize the MPI environment //------------------------------------------------------------------------- MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &myrank); //------------------------------------------------------------------------- // initialize coefficient and fields //------------------------------------------------------------------------- nx = tnx/nprocs + 2; ny = tny; f = (double*)malloc(nx*ny*sizeof(double)); g = (double*)malloc(nx*ny*sizeof(double)); c = (double*)malloc(nx*ny*sizeof(double)); for (i=0; i<nx*ny; i++) { f[i] = 0; g[i] = 0; c[i] = c0; } //------------------------------------------------------------------------- // slit structure //------------------------------------------------------------------------- // We assume the nprocs is multiple of two. if (myrank == nprocs/2) { for (i=1; i<thick+2; i++) { for (j=0; j<ny; j++) c[i*ny+j] = 0; } for (i=1; i<thick+2; i++) { for (j=ny/2-(gap+width)/2-1; j<ny/2+(gap+width)/2; j++) c[i*ny+j] = c0; } for (i=1; i<thick+2; i++) { for (j=ny/2-(gap-width)/2-1; j<ny/2+(gap-width)/2; j++) c[i*ny+j] = 0; } } //------------------------------------------------------------------------- // main loop for the time evolution //------------------------------------------------------------------------- for (tstep=1; tstep<=tmax; tstep++) { //------------------------------------ // point source //------------------------------------ //if (myrank == 0) // g[nx/2-1][ny/2-1] = sin(0.02*tstep); //------------------------------------ // line source //------------------------------------ // We assume the nprocs is multiple of two. if (myrank == nprocs/2 - distance/nx - 1) for (j=0; j<ny; j++) g[(nx-(distance-distance/nx*nx))*ny+j] = sin(0.02*tstep); advance_field(nx, ny, f, g, c); periodic_y(nx, ny, f); exchange_boundary(nx, ny, f); advance_field(nx, ny, g, f, c); periodic_y(nx, ny, g); exchange_boundary(nx, ny, g); } //------------------------------------------------------------------------- // gather fields and save as binary files //------------------------------------------------------------------------- ftot = (double*)malloc(tnx*tny*sizeof(double)); MPI_Gather(&f[n], (nx-2)*ny, MPI_DOUBLE, ftot, (nx-2)*ny, MPI_DOUBLE, 0, MPI_COMM_WORLD); if (myrank == 0) { FILE *fout; fout = fopen("field.bin", "wb"); fwrite(ftot, sizeof(double), tnx*tny, fout); fclose(fout); } //------------------------------------------------------------------------- // finalize the MPI environment //------------------------------------------------------------------------- MPI_Finalize(); return 0; }
myparallel.c
#include<stdio.h> #include<stdlib.h> #include<math.h> #include<mpi.h> void advance_field(int nx, int ny, double *f, double *g, double *c) { int i, j, idx; for (i=1; i<nx-1; i++) { for (j=1; j<ny-1; j++) { idx = i*ny + j; f[idx] = c[idx]*(g[idx+ny] + g[idx-ny] + g[idx+1] + g[idx-1] - 4*g[idx]) + 2*g[idx] - f[idx]; } } } void periodic_y(int nx, int ny, double *f) { int i, j; for (i=0; i<nx; i++) { f[i*ny+(ny-1)] = f[i*ny+1]; f[i*ny+0] = f[i*ny+(ny-2)]; } } void exchange_boundary(int ny, int asize, double *f, int nprocs, int myrank){ int i, j; MPI_Request req[4]; MPI_Status stat[4]; if(myrank < nprocs-1){ MPI_Isend(&f[(asize)*ny], ny, MPI_DOUBLE, myrank+1, 123, MPI_COMM_WORLD, &req[0]); MPI_Irecv(&f[(asize+1)*ny], ny, MPI_DOUBLE, myrank+1, 321, MPI_COMM_WORLD, &req[1]); } if(myrank > 0){ MPI_Irecv(&f[0], ny, MPI_DOUBLE, myrank-1, 123, MPI_COMM_WORLD, &req[2]); MPI_Isend(&f[ny], ny, MPI_DOUBLE, myrank-1, 321, MPI_COMM_WORLD, &req[3]); } if(myrank < nprocs-1){ MPI_Waitall(2, req, stat); } if(myrank > 0){ MPI_Waitall(2, req+2, stat+2); } } int main(int argc, char **argv) { int nx=8000, ny=8000; int tmax=3000; int width=120, thick=60, gap=1200; // slit parameters int distance=800; double c0=0.49; double *f, *g, *c; int i, j, tstep; int nprocs, myrank; int q, r, asize, xspos, xepos; double elapsed_time; //---------------------------- // MPI Initialization //---------------------------- MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &myrank); q = nx/nprocs; r = nx%nprocs; asize = myrank<r?(q+1):q; xspos = myrank<r?(q+1)*myrank:(q*myrank+r); xepos = xspos + asize - 1; //------------------------------------------------------------------------- // initialize coefficient and fields //------------------------------------------------------------------------- f = (double*)malloc((asize+2)*ny*sizeof(double)); g = (double*)malloc((asize+2)*ny*sizeof(double)); c = (double*)malloc((asize+2)*ny*sizeof(double)); for (i=ny; i<(asize+1)*ny; i++) { f[i] = 0; g[i] = 0; c[i] = c0; } //------------------------------------------------------------------------- // slit structure //------------------------------------------------------------------------- int sspot = nx/2-1, espot = nx/2+thick-1, start, end; if(xspos <= sspot && sspot <= xepos){ start = sspot - xspos; if(xspos <= espot && espot <= xepos){ end = espot - xspos; } else end = xepos - xspos; for (i=start; i<=end; i++) { for (j=0; j<ny; j++) c[(i+1)*ny+j] = 0; } for (i=start; i<=end; i++) { for (j=ny/2-(gap+width)/2-1; j<ny/2+(gap+width)/2; j++) c[(i+1)*ny+j] = c0; } for (i=start; i<=end; i++) { for (j=ny/2-(gap-width)/2-1; j<ny/2+(gap-width)/2; j++) c[(i+1)*ny+j] = 0; } } else if(xspos <= espot && espot <= xepos){ start = 0; end = espot - xspos; for (i=start; i<=end; i++) { for (j=0; j<ny; j++) c[(i+1)*ny+j] = 0; } for (i=start; i<=end; i++) { for (j=ny/2-(gap+width)/2-1; j<ny/2+(gap+width)/2; j++) c[(i+1)*ny+j] = c0; } for (i=start; i<=end; i++) { for (j=ny/2-(gap-width)/2-1; j<ny/2+(gap-width)/2; j++) c[(i+1)*ny+j] = 0; } } //------------------------------------------------------------------------- // main loop for the time evolution //------------------------------------------------------------------------- for (tstep=1; tstep<=tmax; tstep++) { //------------------------------------ // point source //------------------------------------ //g[nx/2-1][ny/2-1] = sin(0.02*tstep); //------------------------------------ // line source //------------------------------------ if(myrank == 0 && tstep % 100 == 0) printf("%2.1f percent proceeded\n", tstep*100/(double)tmax); int spot = nx/2-distance-1; if(xspos <= spot && spot <= xepos){ for (j=0; j<ny; j++) g[(spot-xspos+1)*ny+j] = sin(0.02*tstep); } exchange_boundary(ny, asize, g, nprocs, myrank); if(myrank == 0) advance_field(asize+1, ny, f+ny, g+ny, c+ny); else if(myrank == nprocs-1) advance_field(asize+1, ny, f, g, c); else advance_field(asize+2, ny, f, g, c); periodic_y(asize, ny, f+ny); exchange_boundary(ny, asize, f, nprocs, myrank); if(myrank == 0) advance_field(asize+1, ny, g+ny, f+ny, c+ny); else if(myrank == nprocs-1) advance_field(asize+1, ny, g, f, c); else advance_field(asize+2, ny, g, f, c); periodic_y(asize, ny, g+ny); } int *ircnt, *idisp; double *ftot; if (myrank == 0){ printf("nx = %d, nprocs = %d, q = %d, r = %d\n", nx, nprocs, q, r); ircnt = (int*)malloc(sizeof(int)*nprocs); idisp = (int*)malloc(sizeof(int)*nprocs); ftot = (double*)malloc(sizeof(double)*nx*ny); idisp[0] = 0; for(i = 0; i < nprocs; i++){ ircnt[i] = (i<r?(q+1):q)*ny; if(i>0) idisp[i] = idisp[i-1] + ircnt[i-1]; } } MPI_Gatherv(&f[ny], asize*ny, MPI_DOUBLE, ftot, ircnt, idisp, MPI_DOUBLE, 0, MPI_COMM_WORLD); //------------------------------------------------------------------------- // gather fields and save as binary files //------------------------------------------------------------------------- if(myrank == 0) { FILE *fout; fout = fopen("field3.bin", "wb"); fwrite(ftot, sizeof(double), nx*ny, fout); fclose(fout); elapsed_time += MPI_Wtime(); printf("nprocs = %d, elapsed time = %f\n", nprocs, elapsed_time); free(ftot); free(ircnt); free(idisp); } free(f); free(g); MPI_Finalize(); return 0; }
'Parallel Programming' 카테고리의 다른 글
KSC 2014 2번 문제 및 답안 (0) | 2016.10.01 |
---|---|
KSC 2014 1번 문제 및 답안 (0) | 2016.10.01 |
KSC 2013 2번 문제 및 답안 (0) | 2016.10.01 |
KSC 2013 1번 문제 및 답안 (0) | 2016.10.01 |
유한차분법(finite difference method)에서 boundary exchange하기 (0) | 2016.09.28 |
댓글