1 static char help[] = "Tests `GarbageKeyAllReduceIntersect_Private()` in parallel\n\n";
2
3 #include <petscsys.h>
4 #include <petsc/private/garbagecollector.h>
5
6 /* This program tests `GarbageKeyAllReduceIntersect_Private()`.
7 To test this routine in parallel, the sieve of Eratosthenes is
8 implemented.
9 */
10
11 /* Populate an array with Prime numbers <= n.
12 Primes are generated using trial division
13 */
Prime(PetscInt64 ** set,PetscInt64 n)14 PetscErrorCode Prime(PetscInt64 **set, PetscInt64 n)
15 {
16 size_t overestimate;
17 PetscBool is_prime;
18 PetscInt64 ii, jj, count = 0;
19 PetscInt64 *prime;
20
21 PetscFunctionBeginUser;
22 /* There will be fewer than ceil(1.26 * n/log(n)) primes <= n */
23 overestimate = (size_t)PetscCeilReal(((PetscReal)n) * 1.26 / PetscLogReal((PetscReal)n));
24 PetscCall(PetscMalloc1(overestimate, &prime));
25 for (ii = 2; ii < n + 1; ii++) {
26 is_prime = PETSC_TRUE;
27 for (jj = 2; jj <= PetscFloorReal(PetscSqrtReal(ii)); jj++) {
28 if (ii % jj == 0) {
29 is_prime = PETSC_FALSE;
30 break;
31 }
32 }
33 if (is_prime) {
34 prime[count] = ii;
35 count++;
36 }
37 }
38
39 PetscCall(PetscMalloc1((size_t)count + 1, set));
40 (*set)[0] = count;
41 for (ii = 1; ii < count + 1; ii++) (*set)[ii] = prime[ii - 1];
42 PetscCall(PetscFree(prime));
43 PetscFunctionReturn(PETSC_SUCCESS);
44 }
45
46 /* Print out the contents of a set */
PrintSet(MPI_Comm comm,PetscInt64 * set)47 PetscErrorCode PrintSet(MPI_Comm comm, PetscInt64 *set)
48 {
49 char text[64];
50 PetscInt ii;
51
52 PetscFunctionBeginUser;
53 PetscCall(PetscSynchronizedPrintf(comm, "["));
54 for (ii = 1; ii <= (PetscInt)set[0]; ii++) {
55 PetscCall(PetscFormatConvert(" %" PetscInt64_FMT ",", text));
56 PetscCall(PetscSynchronizedPrintf(comm, text, set[ii]));
57 }
58 PetscCall(PetscSynchronizedPrintf(comm, "]\n"));
59 PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
60 PetscFunctionReturn(PETSC_SUCCESS);
61 }
62
63 /* Check set equality */
AssertSetsEqual(PetscInt64 * set,PetscInt64 * true_set)64 PetscErrorCode AssertSetsEqual(PetscInt64 *set, PetscInt64 *true_set)
65 {
66 PetscInt ii;
67
68 PetscFunctionBeginUser;
69 PetscAssert(set[0] == true_set[0], PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Sets of different sizes");
70 for (ii = 1; ii < set[0] + 1; ii++) PetscAssert(set[ii] == true_set[ii], PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Sets are different");
71 PetscFunctionReturn(PETSC_SUCCESS);
72 }
73
74 /* Parallel implementation of the sieve of Eratosthenes */
test_sieve(MPI_Comm comm)75 PetscErrorCode test_sieve(MPI_Comm comm)
76 {
77 PetscInt64 ii, local_p, maximum, n;
78 PetscInt64 *local_set, *cursor, *bootstrap_primes, *truth;
79 PetscMPIInt size, rank;
80 PetscReal x;
81
82 PetscFunctionBeginUser;
83 PetscCallMPI(MPI_Comm_rank(comm, &rank));
84 PetscCallMPI(MPI_Comm_size(comm, &size));
85
86 /* There should be at least `size + 1` primes smaller than
87 (size + 1)*(log(size + 1) + log(log(size + 1))
88 once `size` >=6
89 This is sufficient for each rank to create its own sieve based off
90 a different prime and calculate the size of the sieve.
91 */
92 x = (PetscReal)(size > 6) ? size + 1 : 7;
93 x = x * (PetscLogReal(x) + PetscLogReal(PetscLogReal(x)));
94 PetscCall(Prime(&bootstrap_primes, PetscCeilReal(x)));
95
96 /* Calculate the maximum possible prime, select a prime number for
97 each rank and allocate memorty for the sieve
98 */
99 maximum = bootstrap_primes[size + 1] * bootstrap_primes[size + 1] - 1;
100 local_p = bootstrap_primes[rank + 1];
101 n = maximum - local_p - (maximum / local_p) + 1 + rank + 1;
102 PetscCall(PetscMalloc1(n + 1, &local_set));
103
104 /* Populate the sieve first with all primes <= `local_p`, followed by
105 all integers that are not a multiple of `local_p`
106 */
107 local_set[0] = n;
108 cursor = &local_set[1];
109 for (ii = 0; ii < rank + 1; ii++) {
110 *cursor = bootstrap_primes[ii + 1];
111 cursor++;
112 }
113 for (ii = local_p + 1; ii <= maximum; ii++) {
114 if (ii % local_p != 0) {
115 *cursor = ii;
116 cursor++;
117 }
118 }
119 PetscCall(PetscPrintf(comm, "SIEVES:\n"));
120 PetscCall(PrintSet(comm, local_set));
121
122 PetscCall(PetscFree(bootstrap_primes));
123
124 /* Perform the intersection, testing parallel intersection routine */
125 PetscCall(GarbageKeyAllReduceIntersect_Private(PETSC_COMM_WORLD, &local_set[1], (PetscInt *)&local_set[0]));
126
127 PetscCall(PetscPrintf(comm, "INTERSECTION:\n"));
128 PetscCall(PrintSet(comm, local_set));
129
130 PetscCall(Prime(&truth, maximum));
131 PetscCall(PetscPrintf(comm, "TRUTH:\n"));
132 PetscCall(PrintSet(comm, truth));
133
134 /* Assert the intersection corresponds to primes calculated using
135 trial division
136 */
137 PetscCall(AssertSetsEqual(local_set, truth));
138
139 PetscCall(PetscFree(local_set));
140 PetscCall(PetscFree(truth));
141 PetscFunctionReturn(PETSC_SUCCESS);
142 }
143
144 /* Main executes the individual tests in a predefined order */
main(int argc,char ** argv)145 int main(int argc, char **argv)
146 {
147 PetscFunctionBeginUser;
148 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
149
150 PetscCall(test_sieve(PETSC_COMM_WORLD));
151
152 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ALL PASSED\n"));
153 PetscCall(PetscFinalize());
154 return 0;
155 }
156
157 /*TEST
158
159 testset:
160 test:
161 nsize: 2
162 suffix: 2
163 output_file: output/ex63_2.out
164 test:
165 nsize: 3
166 suffix: 3
167 output_file: output/ex63_3.out
168 test:
169 nsize: 4
170 suffix: 4
171 output_file: output/ex63_4.out
172
173 TEST*/
174