/* get_data2.c -- Parallel Trapezoidal Rule. Uses 3 calls to MPI_Bcast to * distribute input data. * * Input: * a, b: limits of integration. * n: number of trapezoids. * Output: Estimate of the integral from a to b of f(x) * using the trapezoidal rule and n trapezoids. * * Notes: * 1. f(x) is hardwired. * 2. the number of processes (p) should evenly divide * the number of trapezoids (n). * * See Chap. 5, pp. 69 & ff in PPMPI. */ #include /* We'll be using MPI routines, definitions, etc. */ #include "mpi.h" main(int argc, char** argv) { int my_rank; /* My process rank */ int p; /* The number of processes */ float a; /* Left endpoint */ float b; /* Right endpoint */ int n; /* Number of trapezoids */ float h; /* Trapezoid base length */ float local_a; /* Left endpoint my process */ float local_b; /* Right endpoint my process */ int local_n; /* Number of trapezoids for */ /* my calculation */ float integral; /* Integral over my interval */ float total; /* Total integral */ int source; /* Process sending integral */ int dest = 0; /* All messages go to 0 */ int tag = 0; MPI_Status status; void Get_data2(float* a_ptr, float* b_ptr, int* n_ptr, int my_rank); float Trap(float local_a, float local_b, int local_n, float h); /* Calculate local integral */ /* Let the system do what it needs to start up MPI */ MPI_Init(&argc, &argv); /* 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, &p); Get_data2(&a, &b, &n, my_rank); h = (b-a)/n; /* h is the same for all processes */ local_n = n/p; /* So is the number of trapezoids */ /* Length of each process' interval of * integration = local_n*h. So my interval * starts at: */ local_a = a + my_rank*local_n*h; local_b = local_a + local_n*h; integral = Trap(local_a, local_b, local_n, h); /* Add up the integrals calculated by each process */ if (my_rank == 0) { total = integral; for (source = 1; source < p; source++) { MPI_Recv(&integral, 1, MPI_FLOAT, source, tag, MPI_COMM_WORLD, &status); total = total + integral; } } else { MPI_Send(&integral, 1, MPI_FLOAT, dest, tag, MPI_COMM_WORLD); } /* Print the result */ if (my_rank == 0) { printf("With n = %d trapezoids, our estimate\n", n); printf("of the integral from %f to %f = %f\n", a, b, total); } /* Shut down MPI */ MPI_Finalize(); } /* main */ /********************************************************************/ /* Ceiling of log_2(x) is just the number of times * times x-1 can be divided by 2 until the quotient * is 0. Dividing by 2 is the same as right shift. */ int Ceiling_log2(int x /* in */) { /* Use unsigned so that right shift will fill * leftmost bit with 0 */ unsigned temp = (unsigned) x - 1; int result = 0; while (temp != 0) { temp = temp >> 1; result = result + 1 ; } return result; } /* Ceiling_log2 */ /********************************************************************/ int I_receive( int stage /* in */, int my_rank /* in */, int* source_ptr /* out */) { int power_2_stage; /* 2^stage = 1 << stage */ power_2_stage = 1 << stage; if ((power_2_stage <= my_rank) && (my_rank < 2*power_2_stage)){ *source_ptr = my_rank - power_2_stage; return 1; } else return 0; } /* I_receive */ /********************************************************************/ int I_send( int stage /* in */, int my_rank /* in */, int p /* in */, int* dest_ptr /* out */) { int power_2_stage; /* 2^stage = 1 << stage */ power_2_stage = 1 << stage; if (my_rank < power_2_stage){ *dest_ptr = my_rank + power_2_stage; if (*dest_ptr >= p) return 0; else return 1; } else return 0; } /* I_send */ /********************************************************************/ void Send( float a /* in */, float b /* in */, int n /* in */, int dest /* in */) { MPI_Send(&a, 1, MPI_FLOAT, dest, 0, MPI_COMM_WORLD); MPI_Send(&b, 1, MPI_FLOAT, dest, 1, MPI_COMM_WORLD); MPI_Send(&n, 1, MPI_INT, dest, 2, MPI_COMM_WORLD); } /* Send */ /********************************************************************/ void Receive( float* a_ptr /* out */, float* b_ptr /* out */, int* n_ptr /* out */, int source /* in */) { MPI_Status status; MPI_Recv(a_ptr, 1, MPI_FLOAT, source, 0, MPI_COMM_WORLD, &status); MPI_Recv(b_ptr, 1, MPI_FLOAT, source, 1, MPI_COMM_WORLD, &status); MPI_Recv(n_ptr, 1, MPI_INT, source, 2, MPI_COMM_WORLD, &status); } /* Receive */ /********************************************************************/ /* Function Get_data2 * Reads in the user input a, b, and n. * Input parameters: * 1. int my_rank: rank of current process. * 2. int p: number of processes. * Output parameters: * 1. float* a_ptr: pointer to left endpoint a. * 2. float* b_ptr: pointer to right endpoint b. * 3. int* n_ptr: pointer to number of trapezoids. * Algorithm: * 1. Process 0 prompts user for input and * reads in the values. * 2. Process 0 sends input values to other * processes using three calls to MPI_Bcast. */ void Get_data2( float* a_ptr /* out */, float* b_ptr /* out */, int* n_ptr /* out */, int my_rank /* in */) { if (my_rank == 0) { printf("Enter a, b, and n\n"); scanf("%f %f %d", a_ptr, b_ptr, n_ptr); } MPI_Bcast(a_ptr, 1, MPI_FLOAT, 0, MPI_COMM_WORLD); MPI_Bcast(b_ptr, 1, MPI_FLOAT, 0, MPI_COMM_WORLD); MPI_Bcast(n_ptr, 1, MPI_INT, 0, MPI_COMM_WORLD); } /* Get_data2 */ /********************************************************************/ float Trap( float local_a /* in */, float local_b /* in */, int local_n /* in */, float h /* in */) { float integral; /* Store result in integral */ float x; int i; float f(float x); /* function we're integrating */ integral = (f(local_a) + f(local_b))/2.0; x = local_a; for (i = 1; i <= local_n-1; i++) { x = x + h; integral = integral + f(x); } integral = integral*h; return integral; } /* Trap */ /********************************************************************/ float f(float x) { float return_val; /* Calculate f(x). */ /* Store calculation in return_val. */ return_val = x*x; return return_val; } /* f */