// area.c #include <stdio.h> #include <stdlib.h> #include <math.h> #include <omp.h> # define NPOINTS 2000 # define MAXITER 2000 struct complex{ double real; double imag; }; int main(){ int i, j, iter, numoutside = 0, myid, num_threads, start, end, chunk; double area, error, ztemp; struct complex z, c; /* * * * Outer loops run over npoints, initialise z=c * * Inner loop has the iteration z=z*z+c, and threshold test */ #pragma omp parallel default(none) \ private(j,iter,c,ztemp,z,myid,start,end,chunk) \ shared(num_threads) \ reduction(+:numoutside) { myid = omp_get_thread_num(); num_threads = omp_get_num_threads(); chunk = (NPOINTS + num_threads - 1) / num_threads; start = myid * chunk; end = start + chunk < NPOINTS ? start + chunk : NPOINTS; for (i=start; i<end; i++) { for (j=0; j<NPOINTS; j++) { c.real = -2.0+2.5*(double)(i)/(double)(NPOINTS)+1.0e-7; c.imag = 1.125*(double)(j)/(double)(NPOINTS)+1.0e-7; z=c; for (iter=0; iter<MAXITER; iter++){ ztemp=(z.real*z.real)-(z.imag*z.imag)+c.real; z.imag=z.real*z.imag*2+c.imag; z.real=ztemp; if ((z.real*z.real+z.imag*z.imag)>4.0e0) { numoutside++; break; } } } } } /* parallel */ /* * Calculate area and error and output the results */ area=2.0*2.5*1.125*(double)(NPOINTS*NPOINTS-numoutside)/(double)(NPOINTS*NPOINTS); error=area/(double)NPOINTS; printf("Area of Mandlebrot set = %12.8f +/- %12.8f\n",area,error); }