1/** An OpenMP example.
2 *  Based on the example listed on the following web page:
3 *  http://developers.sun.com/sunstudio/downloads/ssx/tha/tha_using.html
4 */
5
6
7#include <assert.h>
8#include <math.h>
9#include <omp.h>
10#include <stdio.h>
11#include <stdlib.h>
12#include <unistd.h>  // getopt()
13#include "../../drd/drd.h"
14
15
16static int is_prime(int* const pflag, int v)
17{
18  int i;
19  int bound = floor(sqrt ((double)v)) + 1;
20
21  for (i = 2; i < bound; i++)
22  {
23    /* No need to check against known composites */
24    if (!pflag[i])
25      continue;
26    if (v % i == 0)
27    {
28      pflag[v] = 0;
29      return 0;
30    }
31  }
32  return (v > 1);
33}
34
35int main(int argc, char **argv)
36{
37  int i;
38  int total = 0;
39  int trace_total = 0;
40  int silent = 0;
41  int n;
42  int num_threads = 2;
43  int optchar;
44  int* primes;
45  int* pflag;
46
47  while ((optchar = getopt(argc, argv, "qt:v")) != EOF)
48  {
49    switch (optchar)
50    {
51    case 'q':
52      silent = 1;
53      break;
54    case 't':
55      num_threads = atoi(optarg);
56      break;
57    case 'v':
58      trace_total = 1;
59      break;
60    default:
61      fprintf(stderr, "Error: unknown option '%c'.\n", optchar);
62      return 1;
63    }
64  }
65
66  if (optind + 1 != argc)
67  {
68    fprintf(stderr, "Error: wrong number of arguments.\n");
69    return 1;
70  }
71  n = atoi(argv[optind]);
72
73  // Not the most user-friendly way to do error checking, but better than
74  // nothing.
75  assert(n > 2);
76  assert(num_threads >= 1);
77
78  primes = malloc(n * sizeof(primes[0]));
79  pflag  = malloc(n * sizeof(pflag[0]));
80
81  omp_set_num_threads(num_threads);
82  omp_set_dynamic(0);
83
84  for (i = 0; i < n; i++) {
85    pflag[i] = 1;
86  }
87
88  if (trace_total)
89    DRD_TRACE_VAR(total);
90
91#pragma omp parallel for
92  for (i = 2; i < n; i++)
93  {
94    if (is_prime(pflag, i))
95    {
96      primes[total] = i;
97      total++;
98    }
99  }
100  if (! silent)
101  {
102    printf("Number of prime numbers between 2 and %d: %d\n",
103           n, total);
104    for (i = 0; i < total; i++)
105    {
106      printf("%d\n", primes[i]);
107    }
108  }
109
110  free(pflag);
111  free(primes);
112
113  return 0;
114}
115