1 #include "mpi.h"
  2 #include <math.h>
  3 #include <stdio.h>
  4 #define MIN(a,b)  ((a)<(b)?(a):(b))
  5 
  6 int main (int argc, char *argv[])
  7 {
  8    int    count;        /* Local prime count */
  9    double elapsed_time; /* Parallel execution time */
 10    int    first;        /* Index of first multiple */
 11    int    global_count; /* Global prime count */
 12    int    high_value;   /* Highest value on this proc */
 13    int    i;
 14    int    id;           /* Process ID number */
 15    int    index;        /* Index of current prime */
 16    int    low_value;    /* Lowest value on this proc */
 17    char  *marked;       /* Portion of 2,...,‘n‘ */
 18    int    n,m;            /* Sieving from 2, ..., ‘n‘ */
 19    int    p;            /* Number of processes */
 20    int    proc0_size;   /* Size of proc 0‘s subarray */
 21    int    prime;        /* Current prime */
 22    int    size;         /* Elements in ‘marked‘ */
 23 
 24    MPI_Init (&argc, &argv);
 25 
 26    /* Start the timer */
 27 
 28    MPI_Comm_rank (MPI_COMM_WORLD, &id);
 29    MPI_Comm_size (MPI_COMM_WORLD, &p);
 30    MPI_Barrier(MPI_COMM_WORLD);
 31    elapsed_time = -MPI_Wtime();
 32 
 33    if (argc != 2) {
 34       if (!id) printf ("Command line: %s <m>\n", argv[0]);
 35       MPI_Finalize();
 36       exit (1);
 37     }
 38 
 39    n = atoi(argv[1]);
 40    m=n;//
 41    n=(n%2==0)?(n/2-1):((n-1)/2);//将输入的整数n转换为存储奇数的数组大小，不包括奇数1
 42    //if (!id) printf ("Number of odd integers:%d    Maximum value of odd integers:%d\n",n+1,3+2*(n-1));
 43    if (n==0) {//输入2时，输出1 prime，结束
 44     if (!id) printf ("There are 1 prime less than or equal to %d\n",m);
 45       MPI_Finalize();
 46       exit (1);
 47     }
 48    /* Figure out this process‘s share of the array, as
 49       well as the integers represented by the first and
 50       last array elements */
 51 
 52    low_value = http://www.mamicode.com/3 + 2*(id*(n)/p);//进程的第一个数
 53    high_value = http://www.mamicode.com/3 + 2*((id+1)*(n)/p-1);//进程的最后一个数
 54    size = (high_value - low_value)/2 + 1;    //进程处理的数组大小
 55 
 56 
 57    /* Bail out if all the primes used for sieving are
 58       not all held by process 0 */
 59 
 60    proc0_size = (n-1)/p;
 61 
 62    if ((3 + 2*(proc0_size-1)) < (int) sqrt((double) (3+2*(n-1)))) {//
 63       if (!id) printf ("Too many processes\n");
 64       MPI_Finalize();
 65       exit (1);
 66     }
 67 
 68    /* Allocate this process‘s share of the array. */
 69 
 70    marked = (char *) malloc (size);
 71 
 72    if (marked == NULL) {
 73       printf ("Cannot allocate enough memory\n");
 74       MPI_Finalize();
 75       exit (1);
 76     }
 77 
 78    for (i = 0; i < size; i++) marked[i] = 0;
 79    if (!id) index = 0;
 80    prime = 3;//从素数3开始
 81    do {
 82     //确定奇数的第一个倍数的下标
 83       if (prime * prime > low_value)
 84          first = (prime * prime - low_value)/2;
 85       else {
 86          if (!(low_value % prime)) 
 87         first = 0;
 88          else 
 89         first=((prime-low_value%prime)%2==0)?((prime-low_value%prime)/2):((prime-low_value%prime+prime)/2);
 90         }
 91 
 92       for (i = first; i < size; i += prime)  marked[i] = 1;
 93       if (!id) {
 94          while (marked[++index]);
 95          prime = 2*index + 3;//下一个未被标记的素数
 96         }
 97       if (p > 1) MPI_Bcast (&prime,  1, MPI_INT, 0, MPI_COMM_WORLD);
 98    } while (prime * prime <= 3+2*(n-1));//
 99 
100    count = 0;
101    for (i = 0; i < size; i++)
102       if (!marked[i])  count++; 
103    if (p > 1) MPI_Reduce (&count, &global_count, 1, MPI_INT, MPI_SUM,
104       0, MPI_COMM_WORLD);
105 
106    /* Stop the timer */
107 
108    elapsed_time += MPI_Wtime();
109 
110 
111    /* Print the results */
112 
113    if (!id) {
114       printf ("There are %d primes less than or equal to %d\n",
115          global_count+1, m);//前面程序是从素数3开始标记，忽略了素数2，所以素数个数要加1
116       printf ("SIEVE (%d) %10.6f\n", p, elapsed_time);
117    }
118    MPI_Finalize ();
119    return 0;
120 }