diff --git a/intro_lab/pi.c b/intro_lab/pi.c index fc01df493601d7e6c4c598048169d93ddb95ff75..4307c6afa22777edeba4e37ee6b85b0469fdd26b 100644 --- a/intro_lab/pi.c +++ b/intro_lab/pi.c @@ -9,7 +9,7 @@ */ int main(int argc, char** argv) { - unsigned long nsteps = 1<<27; /* around 10^8 steps */ + unsigned long nsteps = 100000000; double dx = 1.0 / nsteps; double pi = 0.0; @@ -25,7 +25,7 @@ int main(int argc, char** argv) 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", + printf("pi with %ld steps is %.10f in %.6f seconds (error=%e)\n", nsteps, pi, run_time, fabs(ref_pi - pi)); return 0; diff --git a/intro_lab/pi.f90 b/intro_lab/pi.f90 index 4f4e1def58ea447aec5b3b2e6f2adc9ee8bed931..c04529230f241c88e9ff69176e61ef4371f31d0c 100644 --- a/intro_lab/pi.f90 +++ b/intro_lab/pi.f90 @@ -1,41 +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. - -program calc_pi +! This program computes pi as +! \pi = 4 arctan(1) +! = 4 \int _0 ^1 \frac{1} {1 + x^2} dx + +program compute_pi use omp_lib implicit none -integer(kind=8) :: num_steps = 100000000 -real(kind=8) step - -integer i -real(kind=8) x, pi, raw_sum +integer(kind=8) :: nsteps, i +real(kind=8) :: dx, x, pi, ref_pi real(kind=8) start_time, run_time -step = 1.0D0 / num_steps +nsteps = 100000000 +dx = 1.0D0 / nsteps +pi = 0.0D0 start_time = OMP_GET_WTIME() - -raw_sum = 0.0D0 -do i = 1, num_steps - x = (i-0.5D0)*step - raw_sum = raw_sum + 4.0D0/(1.0D0+x*x) +do i = 1, nsteps + x = (i - 0.5D0) * dx + pi = pi + 1.0D0 / (1.0D0 + x * x) enddo -pi = step * raw_sum +pi = pi * 4.0D0 * dx + run_time = OMP_GET_WTIME() - start_time -print '(" pi with ", i0, " steps is ", f12.6, " in ", f12.6, " seconds")', num_steps, pi, run_time +ref_pi = 4.0D0 * atan(1.0D0) +print '("pi with ", i0, " steps is ", f16.10, " in ", f12.6, " seconds (error=", e12.6, ")")', nsteps, pi, run_time, abs(ref_pi - pi) end program