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;

 22:   PetscFunctionBeginUser;
 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:   PetscCall(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:   PetscCall(PetscMalloc1((size_t)count + 1, set));
 41:   (*set)[0] = count;
 42:   for (ii = 1; ii < count + 1; ii++) { (*set)[ii] = prime[ii - 1]; }
 43:   PetscCall(PetscFree(prime));
 44:   PetscFunctionReturn(PETSC_SUCCESS);
 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;

 53:   PetscFunctionBeginUser;
 54:   PetscCall(PetscSynchronizedPrintf(comm, "["));
 55:   for (ii = 1; ii <= (PetscInt)set[0]; ii++) {
 56:     PetscCall(PetscFormatConvert(" %" PetscInt64_FMT ",", text));
 57:     PetscCall(PetscSynchronizedPrintf(comm, text, set[ii]));
 58:   }
 59:   PetscCall(PetscSynchronizedPrintf(comm, "]\n"));
 60:   PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
 61:   PetscFunctionReturn(PETSC_SUCCESS);
 62: }

 64: /* Check set equality */
 65: PetscErrorCode AssertSetsEqual(PetscInt64 *set, PetscInt64 *true_set)
 66: {
 67:   PetscInt ii;

 69:   PetscFunctionBeginUser;
 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:   PetscFunctionReturn(PETSC_SUCCESS);
 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;

 83:   PetscFunctionBeginUser;
 84:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
 85:   PetscCallMPI(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:   PetscCall(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:   PetscCall(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:   PetscCall(PetscPrintf(comm, "SIEVES:\n"));
121:   PetscCall(PrintSet(comm, local_set));

123:   PetscCall(PetscFree(bootstrap_primes));

125:   /* Perform the intersection, testing parallel intersection routine */
126:   PetscCall(GarbageKeyAllReduceIntersect_Private(PETSC_COMM_WORLD, &local_set[1], (PetscInt *)&local_set[0]));

128:   PetscCall(PetscPrintf(comm, "INTERSECTION:\n"));
129:   PetscCall(PrintSet(comm, local_set));

131:   PetscCall(Prime(&truth, maximum));
132:   PetscCall(PetscPrintf(comm, "TRUTH:\n"));
133:   PetscCall(PrintSet(comm, truth));

135:   /* Assert the intersection corresponds to primes calculated using
136:      trial division
137:   */
138:   PetscCall(AssertSetsEqual(local_set, truth));

140:   PetscCall(PetscFree(local_set));
141:   PetscCall(PetscFree(truth));
142:   PetscFunctionReturn(PETSC_SUCCESS);
143: }

145: /* Main executes the individual tests in a predefined order */
146: int main(int argc, char **argv)
147: {
148:   PetscFunctionBeginUser;
149:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

151:   PetscCall(test_sieve(PETSC_COMM_WORLD));

153:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ALL PASSED\n"));
154:   PetscCall(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*/