查找素数Eratosthenes筛法的mpi程序
思路：
只保留奇数
(1)由输入的整数n确定存储奇数（不包括1）的数组大小：
n=(n%2==0)?(n/2-1):((n-1)/2);//n为存储奇数的数组大小，不包括基数1
(2)由数组大小n、进程号id和进程数p，确定每个进程负责的基数数组的第一个数、最后一个数和数组维度：
low_value = 3 + 2*(id*(n)/p);//进程的第一个数
high_value = 3 + 2*((id+1)*(n)/p-1);//进程的最后一个数
size = (high_value - low_value)/2 + 1; //进程处理的数组大小
(3)寻找奇数的第一个倍数的下标，经过反复思考，有如下规律：
if (prime * prime > low_value)
first = (prime * prime - low_value)/2;
else {
if (!(low_value % prime))
first = 0;
else
first=((prime-low_value%prime)%2==0)?((prime-low_value%prime)/2):((prime-low_value%prime+prime)/2);
}
code:
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 = 3 + 2*(id*(n)/p);//进程的第一个数 53 high_value = 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 }
