diff --git a/intro_lab/README.md b/intro_lab/README.md index 8152b147d9bf925acac196c12fdedc5824dda10d..cb8f9d494c2298323c27a55b620d3047dc709792 100644 --- a/intro_lab/README.md +++ b/intro_lab/README.md @@ -118,20 +118,19 @@ where each rectangle has width Δx and height F(x<sub>i</sub>) at the middl A simple serial C code to calculate π is the following: ``` -static long num_steps = 100000000; -double step; -int main() -{ - int i; - double x, pi, sum = 0.0; - double start_time, run_time; - step = 1.0/(double) num_steps; - for (i=1; i<=num_steps; i++) { - x = (i-0.5)*step; + unsigned long nsteps = 1<<27; /* around 10^8 steps */ + double dx = 1.0 / nsteps; + + double pi = 0.0; + double start_time = omp_get_wtime(); + + unsigned long i; + for (i = 0; i < nsteps; i++) + { + double x = (i + 0.5) * dx; + pi += 1.0 / (1.0 + x * x); } - sum = sum + 4.0/(1.0+x*x); - pi = step * sum; -} + pi *= 4.0 * dx; ``` Instructions: Create a parallel version of the pi.c / pi.f90 program using a @@ -151,35 +150,7 @@ Hints: - Use a parallel construct: ``#pragma omp parallel``. - Divide loop iterations between threads (use the thread ID and the number of threads). - Create an accumulator for each thread to hold partial sums that you can later - combine to generate the global sum, i.e., - -``` -#define MAX_THREADS 4 -... -int main () -{ - double pi, full_sum = 0.0; - double sum[MAX_THREADS]; - for (j=1; j<=MAX_THREADS; j++) { - /* set the number of threads to j and start timing here */ - ... - #pragma omp parallel - { - ... - sum[id] = 0.0; - for (i ...) { - x = (i+0.5)*step; - sum[id] = sum[id] + 4.0/(1.0+x*x); - } - // calculate full_sum iterating over sum[id] - ... - pi = step * full_sum; - /* stop and report timing here */ - ... - } - } -} -``` + combine to generate the global sum. Questions: diff --git a/intro_lab/pi.c b/intro_lab/pi.c index 0d7bedb1db81185a91891f7947642d08e2aef6b3..fc01df493601d7e6c4c598048169d93ddb95ff75 100644 --- a/intro_lab/pi.c +++ b/intro_lab/pi.c @@ -1,43 +1,32 @@ -/* - -This program will numerically compute the integral of - - 4/(1+x*x) - -from 0 to 1. The value of this integral is pi -- which -is great since it gives us an easy way to check the answer. - -The is the original sequential program. It uses the timer -from the OpenMP runtime library - -History: Written by Tim Mattson, 11/99. - -*/ - #include <stdio.h> +#include <math.h> #include <omp.h> -static long num_steps = 100000000; -double step; - -int main() +/* + * This program computes pi as + * \pi = 4 arctan(1) + * = 4 \int _0 ^1 \frac{1} {1 + x^2} dx + */ +int main(int argc, char** argv) { - int i; - double x, pi, sum = 0.0; - double start_time, run_time; - - step = 1.0/(double) num_steps; + unsigned long nsteps = 1<<27; /* around 10^8 steps */ + double dx = 1.0 / nsteps; - start_time = omp_get_wtime(); + double pi = 0.0; + double start_time = omp_get_wtime(); - for (i=1;i<= num_steps; i++){ - x = (i-0.5)*step; - sum = sum + 4.0/(1.0+x*x); + unsigned long i; + for (i = 0; i < nsteps; i++) + { + double x = (i + 0.5) * dx; + pi += 1.0 / (1.0 + x * x); } + pi *= 4.0 * dx; - pi = step * sum; - run_time = omp_get_wtime() - start_time; - printf("\n pi with %ld steps is %lf in %lf seconds\n ",num_steps,pi,run_time); + double run_time = omp_get_wtime() - start_time; + double ref_pi = 4.0 * atan(1.0); + printf("\npi with %ld steps is %.10f in %.6f seconds (error=%e)\n", + nsteps, pi, run_time, fabs(ref_pi - pi)); return 0; }