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