思路:
只保留奇數
(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 }