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