KSC 2013 2번 문제 및 답안

문제 다운로드

문제는 2차 FDM으로 2차 배열을 쪼개어 계산하는 것이다. Exchange Boundary function이 핵심적이라고 볼 수 있다.

sequential.c

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
void advance_field(int nx, int ny, double *f, double *g) {
	int i, j, idx;
	for (i=1; i<nx-1; i++) {
		for (j=1; j<ny-1; j++) {
			idx = i*ny + j;

			f[idx] = 0.49*(g[idx+ny] + g[idx-ny] + g[idx+1] + g[idx-1] - 4*g[idx]) 
			         + 2*g[idx] - f[idx];
		}
	}
}

int main() {
	const int nx=2000, ny=2000;
	int tmax=500;
	double *f, *g;
	int i, j, idx, tstep;

    //-------------------------------------------------------------------------
    // initialize coefficient and fields
    //-------------------------------------------------------------------------
	f = (double*)malloc(nx*ny*sizeof(double));
	g = (double*)malloc(nx*ny*sizeof(double));

	for (i=0; i<nx; i++) {
		for (j=0; j<ny; j++) { 
			f[i*ny + j] = 0;
			g[i*ny + j] = 0;
		}
	}

    //-------------------------------------------------------------------------
  	// main loop for the time evolution
    //-------------------------------------------------------------------------
	for (tstep=1; tstep<=tmax; tstep++) {
    	//------------------------------------
    	// point source
    	//------------------------------------
		idx = (nx/2-1)*ny + (ny/2-1);
	    g[idx] = sin(0.02*tstep);


		advance_field(nx, ny, f, g);
		advance_field(nx, ny, g, f);
	}

    //-------------------------------------------------------------------------
  	// 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) {
	int i, j, idx;

	for (i=1; i<nx-1; i++) {
		for (j=1; j<ny-1; j++) {
			idx = i*ny + j;

			f[idx] = 0.49*(g[idx+ny] + g[idx-ny] + g[idx+1] + g[idx-1] - 4*g[idx]) 
			         + 2*g[idx] - f[idx];
		}
	}
}

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=1000;
	double *f, *g, *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));

	for (i=0; i<nx*ny; i++) {
		f[i] = 0;
		g[i] = 0;
	}

    //-------------------------------------------------------------------------
  	// main loop for the time evolution
    //-------------------------------------------------------------------------
	for (tstep=1; tstep<=tmax; tstep++) {
    	//------------------------------------
    	// point source
    	//------------------------------------
		if (nprocs%2 == 0) {
			if (myrank == nprocs/2)
	    		//g[1*ny+(ny/2-1)] = sin(0.02*tstep);
				for (j=0; j<ny; j++)
	    			g[1*ny+j] = sin(0.02*tstep);
		}
		else {
			if (myrank == nprocs/2)
	    		g[(nx/2-1)*ny+(ny/2-1)] = sin(0.02*tstep);
		}

		advance_field(nx, ny, f, g);
		exchange_boundary(nx, ny, f);

		advance_field(nx, ny, g, f);
		exchange_boundary(nx, ny, g);
	}

    //-------------------------------------------------------------------------
  	// gather fields and save as binary files
    //-------------------------------------------------------------------------
	ftot = (double*)malloc(tnx*tny*sizeof(double));
  	MPI_Gather(&f[1*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) {
	int i, j, idx;
	for (i=1; i<nx-1; i++) {
		for (j=1; j<ny-1; j++) {
			idx = i*ny + j;

			f[idx] = 0.49*(g[idx+ny] + g[idx-ny] + g[idx+1] + g[idx-1] - 4*g[idx]) 
			         + 2*g[idx] - f[idx];
		}
	}
}

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) {
	const int nx=8000, ny=8000;
	int tmax=1000;
	double *f, *g;
	int i, j, idx, tstep;
	double elapsed_time = -MPI_Wtime();

	int nprocs, myrank;
	int q, r, asize, xspos, xepos;
	//---------------------------------
	// 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));

	for (i=1; i<asize+1; i++) {
		for (j=0; j<ny; j++) { 
			f[i*ny + j] = 0;
			g[i*ny + j] = 0;
		}
	}
	

    //-------------------------------------------------------------------------
  	// main loop for the time evolution
    //-------------------------------------------------------------------------
	for (tstep=1; tstep<=tmax; tstep++) {
    	//------------------------------------
    	// point source
    	//------------------------------------
		idx = (nx/2-1)*ny + (ny/2-1);
		if(xspos <= idx/ny && idx/ny <= xepos){
			g[idx-xspos*ny+ny] = sin(0.02*tstep);
			if(tstep==1){
				printf("idx = %d, xspos = %d, xepos = %d, idx/ny = %d\n", idx, xspos, xepos, idx/ny);
				printf("pointing at g[%d] by rank %d\n", idx-xspos*ny+ny, myrank);
			}
		}

		exchange_boundary(ny, asize, g, nprocs, myrank);
		if(myrank == 0)
			advance_field(asize+1, ny, f+ny, g+ny);
		else if(myrank == nprocs-1)
			advance_field(asize+1, ny, f, g);
		else 
			advance_field(asize+2, ny, f, g);
		exchange_boundary(ny, asize, f, nprocs, myrank);
		if(myrank == 0)
			advance_field(asize+1, ny, g+ny, f+ny);
		else if(myrank == nprocs-1)
			advance_field(asize+1, ny, g, f);
		else 
			advance_field(asize+2, ny, g, f);
	}

	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("field2.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;
}