Actual source code: ex63.c
1: static char help[] = "Tests `GarbageKeyAllReduceIntersect_Private()` in parallel\n\n";
3: #include <petscsys.h>
4: #include <petsc/private/garbagecollector.h>
6: /* This program tests `GarbageKeyAllReduceIntersect_Private()`.
7: To test this routine in parallel, the sieve of Eratosthenes is
8: implemented.
9: */
11: /* Populate an array with Prime numbers <= n.
12: Primes are generated using trial division
13: */
14: PetscErrorCode Prime(PetscInt64 **set, PetscInt n)
15: {
16: size_t overestimate;
17: PetscBool is_prime;
18: PetscInt64 ii, jj, count = 0;
19: PetscInt64 *prime;
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: }
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: }
46: /* Print out the contents of a set */
47: PetscErrorCode PrintSet(MPI_Comm comm, PetscInt64 *set)
48: {
49: char text[64];
50: PetscInt ii;
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: }
63: /* Check set equality */
64: PetscErrorCode AssertSetsEqual(PetscInt64 *set, PetscInt64 *true_set)
65: {
66: PetscInt ii;
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: }
74: /* Parallel implementation of the sieve of Eratosthenes */
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;
82: PetscFunctionBeginUser;
83: PetscCallMPI(MPI_Comm_rank(comm, &rank));
84: PetscCallMPI(MPI_Comm_size(comm, &size));
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)));
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));
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));
122: PetscCall(PetscFree(bootstrap_primes));
124: /* Perform the intersection, testing parallel intersection routine */
125: PetscCall(GarbageKeyAllReduceIntersect_Private(PETSC_COMM_WORLD, &local_set[1], (PetscInt *)&local_set[0]));
127: PetscCall(PetscPrintf(comm, "INTERSECTION:\n"));
128: PetscCall(PrintSet(comm, local_set));
130: PetscCall(Prime(&truth, maximum));
131: PetscCall(PetscPrintf(comm, "TRUTH:\n"));
132: PetscCall(PrintSet(comm, truth));
134: /* Assert the intersection corresponds to primes calculated using
135: trial division
136: */
137: PetscCall(AssertSetsEqual(local_set, truth));
139: PetscCall(PetscFree(local_set));
140: PetscCall(PetscFree(truth));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: /* Main executes the individual tests in a predefined order */
145: int main(int argc, char **argv)
146: {
147: PetscFunctionBeginUser;
148: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
150: PetscCall(test_sieve(PETSC_COMM_WORLD));
152: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ALL PASSED\n"));
153: PetscCall(PetscFinalize());
154: return 0;
155: }
157: /*TEST
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
173: TEST*/