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, PetscInt64 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, NULL, 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*/