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