KSC 2013 3번 문제 및 답안

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;
}