#include
<stdlib.h>
#include
<stdio.h>
#include
<math.h>
#include
<signal.h>
#include
"mpi.h"
double
f(double
a){
return
(4.0 / (1.0 + a*a));
}
double
mid_point_rule(int
start,int intervals, double
delta_x){
double
area = 0.0;
int
i;
for
(i = start; i <= start+intervals; i++) {
area += delta_x * f((double
) (delta_x*i+delta_x/2.));
}
return
area ;
}
double
rectangular_rule(int
start,int intervals, double
delta_x){
double
area = 0.0;
int
i;
for
(i = start; i <= start+intervals; i++) {
area += delta_x * f((double
) delta_x*i);
}
return
area ;
}
double
simpsons_rule(int
start,int intervals, double
delta_x){
double
area = f(delta_x*start)-f(delta_x*(start+intervals));
int
i;
for
(i = 1; i <= (intervals/2); i++) {
area += 4.0*f(delta_x*(start+ (2*i-1)))+2.0*f(delta_x*(start+ 2*i));
}
return
delta_x*area/3.0 ;
}
double
(*rule)(int
start,int intervals, double
delta_x);
int
main(int
argc,char *argv[]){
int
done = 0, n, myid, numprocs, mpi_error,namelen;
int
buffer[2],local_start,local_intervals;
double
PI_to_25dp = 3.141592653589793238462643;
double
mypi, pi, delta_x;
double
startwtime, endwtime,beginwtime;
char
processor_name[MPI_MAX_PROCESSOR_NAME];
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
MPI_Comm_rank(MPI_COMM_WORLD,&myid);
MPI_Get_processor_name(processor_name,&namelen);
printf("Process number %d on processor %s \n", myid, processor_name, numprocs);
fflush(stdout);
n = 0;
while
(!done) {
if
(myid == 0) {
beginwtime = MPI_Wtime();
printf("\nWhich rule do you want to use?\n 0= done \n 1=rectangular rule \n 2=mid point rule\n 3=simpsons rule (choose intervals/process(or)s to be even)\n ");
fflush(stdout);
scanf("%d",&buffer[0]);
if
((buffer[0]<1)||(3<buffer[0])) {
done=1;
buffer[0]=0;
}
if
(!done) {
printf("\nPlease input the number of intervals(0 quits).\n (Note: In general, the error is minimized if the number of intervals is a multiple of the number of processors)\n");
fflush(stdout); scanf("%d",&buffer[1]);
}
startwtime = MPI_Wtime();
printf("wall clock begin = %f\n",beginwtime);
fflush(stdout);
} MPI_Bcast(buffer, 2, MPI_INT, 0, MPI_COMM_WORLD);
if
(buffer[0] == 0)
done = 1;
else
{
n=buffer[1];
switch
(buffer[0]) {
case
0:
done=1;
break
;
case
1:
rule=&rectangular_rule;
break
;
case
2:
rule=&mid_point_rule;
break
;
case
3:
rule=&simpsons_rule;
break
;
}
delta_x = 1.0 / (double
) n;
local_start=(n/numprocs)*(myid);
local_intervals=n/numprocs;
printf("process=%d local start = %d local intervals = %d\n",myid,local_start,local_intervals);
fflush(stdout);
mypi=( *rule)(local_start,local_intervals,delta_x);
MPI_Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
if
(myid == 0) {
printf("pi is approximately %.16f, Relative Error is %16.8e\n",
pi, (double
)100 * (pi - PI_to_25dp)/PI_to_25dp);
fflush(stdout);
endwtime = MPI_Wtime();
printf("wall clock time = %f\n", endwtime-startwtime);
fflush(stdout);
}
}
}
mpi_error=MPI_Finalize();
if
(MPI_SUCCESS!=mpi_error)
return
mpi_error;
else
return
0;
}