You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
119 lines
3.5 KiB
C
119 lines
3.5 KiB
C
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#include <math.h>
|
|
#include <stdbool.h>
|
|
#include <mpi.h>
|
|
|
|
const bool debug = false;
|
|
const double a = 0.0;
|
|
const double b = 1.0;
|
|
|
|
/* Function declarations */
|
|
void readinput(int argc, char* argv[], int iam, double* n_p);
|
|
double trapezoidal_rule(double left_endpt, double right_endpt, int trap_count,
|
|
double base_len);
|
|
double f(double x);
|
|
|
|
int main(int argc, char** argv) {
|
|
int iam, np, local_n;
|
|
double n, h, local_a, local_b;
|
|
double local_int, total_int;
|
|
double start, finish, loc_elapsed, elapsed;
|
|
|
|
/* Star the MPI environment and discover who we are */
|
|
MPI_Init(NULL,NULL);
|
|
MPI_Comm_rank(MPI_COMM_WORLD, &iam);
|
|
MPI_Comm_size(MPI_COMM_WORLD, &np);
|
|
|
|
readinput(argc, argv, iam, &n); // Read the user input from the command line
|
|
|
|
/*Note: h and local_n are the same for all processes (load balance!)*/
|
|
h = (b-a)/n; /* basis of each trapezoid */
|
|
local_n = n/np; /* number of trapezoids per process */
|
|
|
|
/* Length of each process' interval of integration = local_n*h. */
|
|
local_a = a + iam*local_n*h;
|
|
local_b = local_a + local_n*h;
|
|
if (debug) printf("Process %d a = %f b = %f\n",iam,local_a,local_b);
|
|
/* Start a timer */
|
|
MPI_Barrier(MPI_COMM_WORLD);
|
|
start = MPI_Wtime();
|
|
/* We compute the local integral on each process using the local endpoints*/
|
|
local_int = trapezoidal_rule(local_a, local_b, local_n, h);
|
|
if (debug) printf("Process %d local integral = %f\n",iam,local_int);
|
|
finish = MPI_Wtime();
|
|
loc_elapsed = finish-start;
|
|
/* Do a max reduce to get the time */
|
|
MPI_Reduce(&loc_elapsed, &elapsed, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
|
|
|
|
/* Do a sum reduce of the integrals computed on each process */
|
|
MPI_Reduce(&local_int, &total_int, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
|
|
|
|
|
|
if (iam == 0) {
|
|
printf("Employing n = %.0f trapezoids\n", n);
|
|
printf("Trapezoid per process n = %d trapezoids\n", local_n);
|
|
printf("Number of processes: %d\n",np);
|
|
printf("of the integral of f(x) from %.0f to %.0f = %1.16f\n", a, b, total_int);
|
|
printf("Total elapsed time = %f milliseconds \n", elapsed * 1000);
|
|
}
|
|
|
|
/* Shut down MPI */
|
|
MPI_Finalize();
|
|
|
|
return 0;
|
|
}
|
|
|
|
void readinput(int argc, char* argv[], int iam, double* n_p){
|
|
/*READINPUT gets the user input: the number of trapezoids
|
|
* Input args: iam: process rank in MPI_COMM_WORLD
|
|
* np: number of processes in MPI_COMM_WORLD
|
|
* Output args: n_p: pointer to number of trapezoids
|
|
*/
|
|
if (iam == 0) {
|
|
if (argc!= 2){
|
|
fprintf(stderr, "usage: srun %s <number of trapezoids> \n", argv[0]);
|
|
fflush(stderr);
|
|
*n_p = -1;
|
|
} else {
|
|
*n_p = atoi(argv[1]);
|
|
}
|
|
}
|
|
/* Broadcasts value of n to each process */
|
|
MPI_Bcast(n_p, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
|
|
|
|
// negative n ends the program
|
|
if (*n_p <= 0) {
|
|
MPI_Finalize();
|
|
exit(-1); // we return an error!
|
|
}
|
|
}
|
|
|
|
double trapezoidal_rule(double left_endpt, double right_endpt, int trap_count, double base_len) {
|
|
/* TRAPEZOIDAL_RULE the serial function applying the trapezoidal rule
|
|
* Input: left_endpt
|
|
* right_endpt
|
|
* trap_count
|
|
* base_len
|
|
* Output: Approximated value of the integral
|
|
*/
|
|
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;
|
|
}
|
|
|
|
|
|
double f(double x) {
|
|
/*F function to be integrated */
|
|
return (4.0*sqrt(1.0 - x*x));
|
|
}
|