1 2 static char help[] = "Tests `GarbageKeyAllReduceIntersect_Private()` in parallel\n\n"; 3 4 #include <petscsys.h> 5 #include <petsc/private/garbagecollector.h> 6 7 /* This program tests `GarbageKeyAllReduceIntersect_Private()`. 8 To test this routine in parallel, the sieve of Eratosthenes is 9 implemented. 10 */ 11 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; 21 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 } 39 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 } 46 47 /* Print out the contents of a set */ 48 PetscErrorCode PrintSet(MPI_Comm comm, PetscInt64 *set) 49 { 50 char text[64]; 51 PetscInt ii; 52 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 } 63 64 /* Check set equality */ 65 PetscErrorCode AssertSetsEqual(PetscInt64 *set, PetscInt64 *true_set) 66 { 67 PetscInt ii; 68 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 } 74 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; 82 83 PetscFunctionBeginUser; 84 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 85 PetscCallMPI(MPI_Comm_size(comm, &size)); 86 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))); 96 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)); 104 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)); 122 123 PetscCall(PetscFree(bootstrap_primes)); 124 125 /* Perform the intersection, testing parallel intersection routine */ 126 PetscCall(GarbageKeyAllReduceIntersect_Private(PETSC_COMM_WORLD, &local_set[1], (PetscInt *)&local_set[0])); 127 128 PetscCall(PetscPrintf(comm, "INTERSECTION:\n")); 129 PetscCall(PrintSet(comm, local_set)); 130 131 PetscCall(Prime(&truth, maximum)); 132 PetscCall(PetscPrintf(comm, "TRUTH:\n")); 133 PetscCall(PrintSet(comm, truth)); 134 135 /* Assert the intersection corresponds to primes calculated using 136 trial division 137 */ 138 PetscCall(AssertSetsEqual(local_set, truth)); 139 140 PetscCall(PetscFree(local_set)); 141 PetscCall(PetscFree(truth)); 142 PetscFunctionReturn(PETSC_SUCCESS); 143 } 144 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)); 150 151 PetscCall(test_sieve(PETSC_COMM_WORLD)); 152 153 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ALL PASSED\n")); 154 PetscCall(PetscFinalize()); 155 return 0; 156 } 157 158 /*TEST 159 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 173 174 TEST*/ 175