-
Notifications
You must be signed in to change notification settings - Fork 0
/
mpi_trap1.c
146 lines (120 loc) · 4.18 KB
/
mpi_trap1.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
#include <stdio.h>
/* We'll be using MPI routines, definitions, etc. */
#include <mpi.h>
/* Get the input values */
void Get_input(int my_rank, int comm_sz, double* a_p, double* b_p,
int* n_p);
/* Calculate local integral */
double Trap(double left_endpt, double right_endpt, int trap_count,
double base_len);
/* Function we're integrating */
double f(double x);
int main(void)
{
int my_rank, comm_sz, n, local_n;
double a, b, h, local_a, local_b;
double local_int, total_int;
int source;
/* Let the system do what it needs to start up MPI */
MPI_Init(NULL, NULL);
/* Get my process rank */
MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
/* Find out how many processes are being used */
MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
Get_input(my_rank, comm_sz, &a, &b, &n);
h = (b - a) / n; /* h is the same for all processes */
local_n = n / comm_sz; /* So is the number of trapezoids */
local_a = a + my_rank * local_n * h;
local_b = local_a + local_n * h;
local_int = Trap(local_a, local_b, local_n, h);
/* Add up the integrals calculated by each process */
if (my_rank != 0)
{
MPI_Send(&local_int, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
} else {
total_int = local_int;
for ( source = 1; source < comm_sz; source++)
{
MPI_Recv(&local_int, 1, MPI_DOUBLE, source, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
total_int += local_int;
}
}
/* Print the result */
if (my_rank == 0) {
printf("With n = %d trapezoids, our estimate\n", n);
printf("of the integral from %f to %f = %.15e\n",
a, b, total_int);
}
/* Shut down MPI */
MPI_Finalize();
return 0;
} /* main */
/*------------------------------------------------------------------
* Function: Get_input
* Purpose: Get the user input: the left and right endpoints
* and the number of trapezoids
* Input args: my_rank: process rank in MPI_COMM_WORLD
* comm_sz: number of processes in MPI_COMM_WORLD
* Output args: a_p: pointer to left endpoint
* b_p: pointer to right endpoint
* n_p: pointer to number of trapezoids
*/
void Get_input(
int my_rank, /* in */
int comm_sz, /* in */
double* a_p, /* out */
double* b_p, /* out */
int* n_p /* out */) {
int dest;
if (my_rank == 0) {
printf("Enter a, b, and n \n");
scanf("%lf %lf %d", a_p, b_p, n_p);
for ( dest = 1; dest < comm_sz; dest++) {
MPI_Send(a_p, 1, MPI_DOUBLE, dest, 0, MPI_COMM_WORLD);
MPI_Send(b_p, 1, MPI_DOUBLE, dest, 0, MPI_COMM_WORLD);
MPI_Send(n_p, 1, MPI_INT, dest, 0, MPI_COMM_WORLD);
}
} else { /* my_rank !=0 */
MPI_Recv(a_p, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
MPI_Recv(b_p, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
MPI_Recv(n_p, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
}
} /* Get_input */
/*------------------------------------------------------------------
* Function: Trap
* Purpose: Serial function for estimating a definite integral
* using the trapezoidal rule
* Input args: left_endpt
* right_endpt
* trap_count
* base_len
* Return val: Trapezoidal rule estimate of integral from
* left_endpt to right_endpt using trap_count
* trapezoids
*/
double Trap(
double left_endpt /* in */,
double right_endpt /* in */,
int trap_count /* in */,
double base_len /* in */)
{
double estimate, x;
int i;
estimate = (f(left_endpt) + f(right_endpt)) / 2.0;
for (i = 1; i <= trap_count - 1; i++)
{
x = left_endpt + i * base_len;
estimate += f(x);
}
estimate = estimate * base_len;
return estimate;
} /* Trap */
/*------------------------------------------------------------------
* Function: f
* Purpose: Compute value of function to be integrated
* Input args: x
*/
double f(double x)
{
return x*x;
} /* f */