1 #include <petsc/private/randomimpl.h> 2 #include <Random123/threefry.h> 3 4 /* The structure of the Random123 methods are similar enough that templates could be used to make the other CBRNGs in 5 * the package (aes, ars, philox) available, as well as different block sizes. But threefry4x64 is a good default, 6 * and I'd rather get a simple implementation up and working and come back if there's interest. */ 7 typedef struct _n_PetscRandom123 { 8 threefry4x64_ctr_t counter; 9 threefry4x64_key_t key; 10 threefry4x64_ctr_t result; 11 PetscInt count; 12 } PetscRandom123; 13 14 R123_ULONG_LONG PETSCR123_SEED_0 = R123_64BIT(0x615D333D2655FE14); 15 R123_ULONG_LONG PETSCR123_SEED_1 = R123_64BIT(0xAFF6369B3EE9FE96); 16 R123_ULONG_LONG PETSCR123_SEED_2 = R123_64BIT(0x5956EBC717B60E07); 17 R123_ULONG_LONG PETSCR123_SEED_3 = R123_64BIT(0xEE8612A0CBEABFF1); 18 19 static PetscErrorCode PetscRandomSeed_Random123(PetscRandom r) 20 { 21 threefry4x64_ukey_t ukey; 22 PetscRandom123 *r123 = (PetscRandom123 *)r->data; 23 24 PetscFunctionBegin; 25 ukey.v[0] = (R123_ULONG_LONG)r->seed; 26 ukey.v[1] = PETSCR123_SEED_1; 27 ukey.v[2] = PETSCR123_SEED_2; 28 ukey.v[3] = PETSCR123_SEED_3; 29 /* The point of seeding should be that every time the sequence is seeded you get the same output. In this CBRNG, 30 * that means we have to initialize the key and reset the counts */ 31 r123->key = threefry4x64keyinit(ukey); 32 r123->counter.v[0] = 0; 33 r123->counter.v[1] = 1; 34 r123->counter.v[2] = 2; 35 r123->counter.v[3] = 3; 36 r123->result = threefry4x64(r123->counter, r123->key); 37 r123->count = 0; 38 PetscFunctionReturn(PETSC_SUCCESS); 39 } 40 41 static PetscReal PetscRandom123Step(PetscRandom123 *r123) 42 { 43 PetscReal scale = 1. / ((PetscReal)UINT64_MAX + 1.); 44 PetscReal shift = .5 * scale; 45 PetscInt mod = (r123->count++) % 4; 46 PetscReal ret; 47 48 ret = r123->result.v[mod] * scale + shift; 49 50 if (mod == 3) { 51 r123->counter.v[0] += 4; 52 r123->counter.v[1] += 4; 53 r123->counter.v[2] += 4; 54 r123->counter.v[3] += 4; 55 r123->result = threefry4x64(r123->counter, r123->key); 56 } 57 58 return ret; 59 } 60 61 static PetscErrorCode PetscRandomGetValue_Random123(PetscRandom r, PetscScalar *val) 62 { 63 PetscRandom123 *r123 = (PetscRandom123 *)r->data; 64 PetscScalar rscal; 65 66 PetscFunctionBegin; 67 #if defined(PETSC_USE_COMPLEX) 68 { 69 PetscReal re = PetscRandom123Step(r123); 70 PetscReal im = PetscRandom123Step(r123); 71 72 if (r->iset) { 73 re = re * PetscRealPart(r->width) + PetscRealPart(r->low); 74 im = im * PetscImaginaryPart(r->width) + PetscImaginaryPart(r->low); 75 } 76 77 rscal = PetscCMPLX(re, im); 78 } 79 #else 80 rscal = PetscRandom123Step(r123); 81 if (r->iset) rscal = rscal * r->width + r->low; 82 #endif 83 *val = rscal; 84 PetscFunctionReturn(PETSC_SUCCESS); 85 } 86 87 static PetscErrorCode PetscRandomGetValueReal_Random123(PetscRandom r, PetscReal *val) 88 { 89 PetscRandom123 *r123 = (PetscRandom123 *)r->data; 90 PetscReal rreal; 91 92 PetscFunctionBegin; 93 rreal = PetscRandom123Step(r123); 94 if (r->iset) rreal = rreal * PetscRealPart(r->width) + PetscRealPart(r->low); 95 *val = rreal; 96 PetscFunctionReturn(PETSC_SUCCESS); 97 } 98 99 static PetscErrorCode PetscRandomGetValuesReal_Random123(PetscRandom r, PetscInt n, PetscReal vals[]) 100 { 101 PetscRandom123 *r123 = (PetscRandom123 *)r->data; 102 PetscInt peel_start; 103 PetscInt rem, lim; 104 PetscReal scale = 1. / ((PetscReal)UINT64_MAX + 1.); 105 PetscReal shift = .5 * scale; 106 PetscRandom123 r123_copy; 107 108 PetscFunctionBegin; 109 peel_start = (4 - (r123->count % 4)) % 4; 110 peel_start = PetscMin(n, peel_start); 111 for (PetscInt i = 0; i < peel_start; i++) PetscCall(PetscRandomGetValueReal(r, &vals[i])); 112 PetscAssert((r123->count % 4) == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad modular arithmetic"); 113 n -= peel_start; 114 vals += peel_start; 115 rem = (n % 4); 116 lim = n - rem; 117 if (r->iset) { 118 scale *= PetscRealPart(r->width); 119 shift *= PetscRealPart(r->width); 120 shift += PetscRealPart(r->low); 121 } 122 r123_copy = *r123; 123 for (PetscInt i = 0; i < lim; i += 4, vals += 4) { 124 vals[0] = r123_copy.result.v[0] * scale + shift; 125 vals[1] = r123_copy.result.v[1] * scale + shift; 126 vals[2] = r123_copy.result.v[2] * scale + shift; 127 vals[3] = r123_copy.result.v[3] * scale + shift; 128 r123_copy.counter.v[0] += 4; 129 r123_copy.counter.v[1] += 4; 130 r123_copy.counter.v[2] += 4; 131 r123_copy.counter.v[3] += 4; 132 r123_copy.result = threefry4x64(r123->counter, r123->key); 133 } 134 r123_copy.count += lim; 135 *r123 = r123_copy; 136 for (PetscInt i = 0; i < rem; i++) PetscCall(PetscRandomGetValueReal(r, &vals[i])); 137 PetscFunctionReturn(PETSC_SUCCESS); 138 } 139 140 static PetscErrorCode PetscRandomGetValues_Random123(PetscRandom r, PetscInt n, PetscScalar vals[]) 141 { 142 PetscFunctionBegin; 143 #if PetscDefined(USE_COMPLEX) 144 for (PetscInt i = 0; i < n; i++) PetscCall(PetscRandomGetValue_Random123(r, n, &vals[i])); 145 #else 146 PetscCall(PetscRandomGetValuesReal_Random123(r, n, vals)); 147 #endif 148 PetscFunctionReturn(PETSC_SUCCESS); 149 } 150 151 static PetscErrorCode PetscRandomDestroy_Random123(PetscRandom r) 152 { 153 PetscFunctionBegin; 154 PetscCall(PetscFree(r->data)); 155 PetscFunctionReturn(PETSC_SUCCESS); 156 } 157 158 static struct _PetscRandomOps PetscRandomOps_Values = { 159 // clang-format off 160 PetscDesignatedInitializer(seed, PetscRandomSeed_Random123), 161 PetscDesignatedInitializer(getvalue, PetscRandomGetValue_Random123), 162 PetscDesignatedInitializer(getvaluereal, PetscRandomGetValueReal_Random123), 163 PetscDesignatedInitializer(getvalues, PetscRandomGetValues_Random123), 164 PetscDesignatedInitializer(getvaluesreal, PetscRandomGetValuesReal_Random123), 165 PetscDesignatedInitializer(destroy, PetscRandomDestroy_Random123), 166 // clang-format on 167 }; 168 169 /*MC 170 PETSCRANDOM123 - access to Random123 counter based pseudorandom number generators (currently threefry4x64) 171 172 Options Database Key: 173 . -random_type <rand,rand48,sprng,random123> - select the random number generator at runtim 174 175 Level: beginner 176 177 Note: 178 PETSc must be ./configure with the option --download-random123 to use this random number generator. 179 180 .seealso: `RandomCreate()`, `RandomSetType()`, `PETSCRAND`, `PETSCRAND48`, `PETSCSPRNG`, `PetscRandomSetFromOptions()` 181 M*/ 182 183 PETSC_EXTERN PetscErrorCode PetscRandomCreate_Random123(PetscRandom r) 184 { 185 PetscRandom123 *r123; 186 187 PetscFunctionBegin; 188 PetscCall(PetscNew(&r123)); 189 r->data = r123; 190 r->ops[0] = PetscRandomOps_Values; 191 PetscCall(PetscObjectChangeTypeName((PetscObject)r, PETSCRANDOM123)); 192 PetscCall(PetscRandomSeed(r)); 193 PetscFunctionReturn(PETSC_SUCCESS); 194 } 195