KSC 2013 1번 문제 및 답안

문제1

Simpson 적분공식을 활용한 1차원 정적분 프로그램을 고려한다. 주어진 시작 프로그램(C 또는 FORTRAN)을 병렬화(MPI/OpenMP)시켜서 주어진 시스템에서 적분값과 총 함수 호출 횟수를 출력하시오. 또한, 최대 병렬 효율 성능을 구하시오. 여기서 사용된 Simpson 정적분 공식은 구간 [a,b]에 대한 정적분을 의미하고, n+1개의 함수값을 호출하게 되어 있다. 여기서 n은 짝수 이다. 실제 정적분값은 다음의 수식으로 근사된다. (h/3)[f(x0)+4f(x1)+2f(x2)+....+4f(xn-1)+f(xn)]. 여기서, h=(b-a)/n 이며, xi=a+ih이고, f(x)=4/(1+x2)이다. 주의할 점은 함수 호출을 n+1번만 수행한다. 

Simpson 적분공식은 검색해보면 이와 같다. 다시말해, x0와 xn을 제외한 부분을 1차 FDM으로 계산하면 된다.

sequential.c

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<sys/time.h>

int ncall;
typedef enum { FALSE=0, TRUE=1} boolean;
float gettime();

double func(double x){
	ncall ++;
	return 4.L/(x*x+1.L);
}

double simpson0(double (*func1)(double),int n,double aa,double bb){
	double h,xx,rslt;
	boolean lodd;
	int j;

	if((n%2) !=0){
		fprintf(stderr,"input error, n must be even number %d\n",n);
		exit(9);
	}
	rslt = 0.L;
	h = (bb-aa)/(double)n;
	rslt = (func1(aa) + func1(bb));
	for(lodd = TRUE,j = 1; j<n;j++){
		xx = aa + h * (double)j;
		if(lodd) {
			rslt +=  4.L*func1(xx);
		}
		else {
			rslt += 2.L*func1(xx);
		}
		lodd = ((lodd)? (FALSE):(TRUE));
	}
	rslt = rslt*h/3.L;
	return rslt;
}

int main(int argc, char **argv){
	double rslt,aa,bb;
	int n;
	double time_start,time_end;

	time_start = gettime();
	n = 2000000000;
	aa = 0.L;
	bb = 1.L;
	ncall = 0;

	rslt = simpson0(func,n,aa,bb);
	time_end = gettime();
	printf("%d %18.16g n, rslt,%d\n",n,rslt,ncall);
	printf("Wallclock time %g\n",time_end-time_start);
	return 0;
}
struct timeval tv;
float gettime()
{
	static int startflag=1;
	static double tsecs0, tsecs1;	
	if( startflag ) {
		(void ) gettimeofday(&tv,0);
		tsecs0 = tv.tv_sec + tv.tv_usec*1.0e-6;
		startflag = 0;
	}
	(void) gettimeofday(&tv,0);
	tsecs1 = tv.tv_sec + tv.tv_usec*1.0e-6;
	return ((float) (tsecs1-tsecs0));
}

answer.c

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#include<mpi.h>


#define min(a,b) ((a)<(b)?(a):(b))


int ncall;

double simpson2(double (*func1)(double),int , double , double , int , int );

double func(double x){
	ncall ++;
	return 4.L/(x*x+1.L);
}

int main(int argc, char **argv){
	double rslt,aa,bb,rslt0;
	int n,ncall0;
	int ierr, kount, iroot;
	int myid,nproc;
	double time1,time2,time_start, time_end;

	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD,&myid);
	MPI_Comm_size(MPI_COMM_WORLD,&nproc);

	if(myid ==0 && nproc > 1) printf("%d processes are alive\n",nproc);
	if(myid ==0 && nproc == 1) printf("%d processes is alive\n",nproc);


	time_start = MPI_Wtime();
	n=20;
	n =6; 
	n  = 200000;
	n = 2000000000;

	aa = 0.L;
	bb = 1.L;
	ncall = 0;
	ncall = 0;

	rslt = simpson2(func,n,aa,bb,myid,nproc);

	iroot = 0; kount = 1;
	MPI_Reduce(&rslt,&rslt0, kount, MPI_DOUBLE, MPI_SUM,iroot, MPI_COMM_WORLD);
	if(myid==0) printf("%d %18.14g n, rslt0\n",n,rslt0);
	iroot = 0; kount = 1;
	MPI_Reduce(&ncall, &ncall0,kount,MPI_INT,MPI_SUM,iroot, MPI_COMM_WORLD);
	if(myid==0) printf("%d ncall0\n",ncall0);

	time_end = MPI_Wtime();
	if(myid==0) {
		double timer = time_end - time_start;
		double timerm = timer/60.L;
		double timerh = timer/3600.L;
		double timerd = timer/3600.L/24.L;
		printf("%14.5g s %14.5g m %14.5g h %14.5g d\n",timer, timerm,timerh,timerd);
	}
	MPI_Finalize();
	return 0;
}

void equal_load(int n1, int n2, int nproc, int myid, int *istart, int *ifinish){
	int iw1,iw2;
	iw1 = (n2-n1+1)/nproc;
	iw2 = (n2-n1+1)%nproc;
	*istart = myid*iw1+n1+min(myid,iw2);
	*ifinish= *istart +iw1 -1;
	if(iw2>myid) *ifinish = *ifinish + 1;
	if(n2<*istart) *ifinish = *istart - 1;
}

double simpson2(double (*func1)(double),int n, double aa, double bb, int myid, int nproc){
	double rslt;
	double h,xx;
	int j,n1,n2;
	int istart,ifinish;

	rslt = 0.; 
	if((n%2) != 0) {
		fprintf(stderr,"input error, n must be even number %d\n",n);
		exit(9);
	}

	n1 = 1; n2 = n-1;
	equal_load(n1,n2,nproc,myid,&istart,&ifinish);


	h = (bb-aa)/(double)n;
	if(myid==0) rslt = (func1(aa)+func1(bb));
	for(j=istart;j<=ifinish;j++){
		xx = aa + h*(double)j;
		/*
		if((j%2) == 1) {
		*/
		if((j&1)) {
			rslt += 4.L*func1(xx);
		}
		else {
			rslt += 2.L*func1(xx);
		}
	}
	rslt = rslt * h/3.L;
	return rslt;
}

myparallel.c

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<sys/time.h>
#include<mpi.h>

int ncall;
typedef enum { FALSE=0, TRUE=1} boolean;
float gettime();

double func(double x){
	ncall ++;
	return 4.L/(x*x+1.L);
}

double simpson0(double (*func1)(double),int n,double aa,double bb, int ispos, int iepos){
	double h,xx,rslt;
	int j;

	rslt = 0.L;
	h = (bb-aa)/(double)n;
	for(j = ispos; j <= iepos; j++){
		xx = aa + h * (double)j;
		if(j%2 == 1) {
			rslt +=  4.L*func1(xx);
		}
		else {
			rslt += 2.L*func1(xx);
		}
	}
	//rslt = rslt*h/3.L;
	return rslt;
}

int main(int argc, char **argv){
	double rslt,aa,bb,h,totrslt;
	int n, totcall;
	double time_start,time_end;
	time_start = gettime();
	n = 2000000000;
	aa = 0.L;
	bb = 1.L;
	ncall = 0;

	int nprocs, myrank;
	int q, r, asize, ispos, iepos;
	MPI_Init(&argc, &argv);
	MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

	if(myrank == 0 && (n%2) !=0){
		fprintf(stderr,"input error, n must be even number %d\n",n);
		exit(9);
	}

	q = (n+1)/nprocs;
	r = (n+1)%nprocs;
	asize = myrank<r?(q+1):q;
	ispos = myrank<r?(q+1)*myrank:(q*myrank+r);
	iepos = ispos + asize - 1;
	if(myrank == 0) ++ispos;
	if(myrank == nprocs-1) --iepos;

	rslt = 0.L;
	h = (bb-aa)/(double)n;
	if(myrank == 0) rslt += (func(aa) + func(bb));
	rslt += simpson0(func,n,aa,bb,ispos,iepos);
	MPI_Reduce(&rslt, &totrslt, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
	MPI_Reduce(&ncall, &totcall, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);

	if(myrank == 0)
		totrslt = totrslt*h/3.L;
	
	time_end = gettime();
	if(myrank == 0){
		printf("%d %18.16g n, rslt,%d\n",n,totrslt,totcall);
		printf("Wallclock time %g\n",time_end-time_start);
	}
	MPI_Finalize();
	return 0;
}
struct timeval tv;
float gettime()
{
	static int startflag=1;
	static double tsecs0, tsecs1;	
	if( startflag ) {
		(void ) gettimeofday(&tv,0);
		tsecs0 = tv.tv_sec + tv.tv_usec*1.0e-6;
		startflag = 0;
	}
	(void) gettimeofday(&tv,0);
	tsecs1 = tv.tv_sec + tv.tv_usec*1.0e-6;
	return ((float) (tsecs1-tsecs0));
}