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 PetscErrorCode PetscRandomSeed_Random123(PetscRandom r) { 20 threefry4x64_ukey_t ukey; 21 PetscRandom123 *r123 = (PetscRandom123 *)r->data; 22 23 PetscFunctionBegin; 24 ukey.v[0] = (R123_ULONG_LONG)r->seed; 25 ukey.v[1] = PETSCR123_SEED_1; 26 ukey.v[2] = PETSCR123_SEED_2; 27 ukey.v[3] = PETSCR123_SEED_3; 28 /* The point of seeding should be that every time the sequence is seeded you get the same output. In this CBRNG, 29 * that means we have to initialize the key and reset the counts */ 30 r123->key = threefry4x64keyinit(ukey); 31 r123->counter.v[0] = 0; 32 r123->counter.v[1] = 1; 33 r123->counter.v[2] = 2; 34 r123->counter.v[3] = 3; 35 r123->result = threefry4x64(r123->counter, r123->key); 36 r123->count = 0; 37 PetscFunctionReturn(0); 38 } 39 40 static PetscReal PetscRandom123Step(PetscRandom123 *r123) { 41 PetscReal scale = ((PetscReal)1.) / (UINT64_MAX + ((PetscReal)1.)); 42 PetscReal shift = .5 * scale; 43 PetscInt mod = (r123->count++) % 4; 44 PetscReal ret; 45 46 ret = r123->result.v[mod] * scale + shift; 47 48 if (mod == 3) { 49 r123->counter.v[0] += 4; 50 r123->counter.v[1] += 4; 51 r123->counter.v[2] += 4; 52 r123->counter.v[3] += 4; 53 r123->result = threefry4x64(r123->counter, r123->key); 54 } 55 56 return ret; 57 } 58 59 PetscErrorCode PetscRandomGetValue_Random123(PetscRandom r, PetscScalar *val) { 60 PetscRandom123 *r123 = (PetscRandom123 *)r->data; 61 PetscScalar rscal; 62 63 PetscFunctionBegin; 64 #if defined(PETSC_USE_COMPLEX) 65 { 66 PetscReal re = PetscRandom123Step(r123); 67 PetscReal im = PetscRandom123Step(r123); 68 69 if (r->iset) { 70 re = re * PetscRealPart(r->width) + PetscRealPart(r->low); 71 im = im * PetscImaginaryPart(r->width) + PetscImaginaryPart(r->low); 72 } 73 74 rscal = PetscCMPLX(re, im); 75 } 76 #else 77 rscal = PetscRandom123Step(r123); 78 if (r->iset) rscal = rscal * r->width + r->low; 79 #endif 80 *val = rscal; 81 PetscFunctionReturn(0); 82 } 83 84 PetscErrorCode PetscRandomGetValueReal_Random123(PetscRandom r, PetscReal *val) { 85 PetscRandom123 *r123 = (PetscRandom123 *)r->data; 86 PetscReal rreal; 87 88 PetscFunctionBegin; 89 rreal = PetscRandom123Step(r123); 90 if (r->iset) rreal = rreal * PetscRealPart(r->width) + PetscRealPart(r->low); 91 *val = rreal; 92 PetscFunctionReturn(0); 93 } 94 95 PetscErrorCode PetscRandomDestroy_Random123(PetscRandom r) { 96 PetscFunctionBegin; 97 PetscCall(PetscFree(r->data)); 98 PetscFunctionReturn(0); 99 } 100 101 static struct _PetscRandomOps PetscRandomOps_Values = { 102 PetscDesignatedInitializer(seed, PetscRandomSeed_Random123), 103 PetscDesignatedInitializer(getvalue, PetscRandomGetValue_Random123), 104 PetscDesignatedInitializer(getvaluereal, PetscRandomGetValueReal_Random123), 105 PetscDesignatedInitializer(getvalues, NULL), 106 PetscDesignatedInitializer(getvaluesreal, NULL), 107 PetscDesignatedInitializer(destroy, PetscRandomDestroy_Random123), 108 }; 109 110 /*MC 111 PETSCRANDOM123 - access to Random123 counter based pseudorandom number generators (currently threefry4x64) 112 113 Options Database Keys: 114 . -random_type <rand,rand48,sprng,random123> - select the random number generator at runtim 115 116 Note: 117 PETSc must be ./configure with the option --download-random123 to use this random number generator. 118 119 Level: beginner 120 121 .seealso: `RandomCreate()`, `RandomSetType()`, `PETSCRAND`, `PETSCRAND48`, `PETSCSPRNG`, `PetscRandomSetFromOptions()` 122 M*/ 123 124 PETSC_EXTERN PetscErrorCode PetscRandomCreate_Random123(PetscRandom r) { 125 PetscRandom123 *r123; 126 127 PetscFunctionBegin; 128 PetscCall(PetscNewLog(r, &r123)); 129 r->data = r123; 130 PetscCall(PetscMemcpy(r->ops, &PetscRandomOps_Values, sizeof(PetscRandomOps_Values))); 131 PetscCall(PetscObjectChangeTypeName((PetscObject)r, PETSCRANDOM123)); 132 PetscCall(PetscRandomSeed(r)); 133 PetscFunctionReturn(0); 134 } 135