Actual source code: random123.c
petsc-3.11.4 2019-09-28
1: #include <../src/sys/classes/random/randomimpl.h>
2: #include <Random123/threefry.h>
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;
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);
21: PetscErrorCode PetscRandomSeed_Random123(PetscRandom r)
22: {
23: threefry4x64_ukey_t ukey;
24: PetscRandom123 *r123 = (PetscRandom123 *) r->data;
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: return(0);
41: }
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;
50: ret = r123->result.v[mod] * scale + shift;
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: }
60: return ret;
61: }
63: PetscErrorCode PetscRandomGetValue_Random123(PetscRandom r,PetscScalar *val)
64: {
65: PetscRandom123 *r123 = (PetscRandom123 *) r->data;
66: PetscScalar rscal;
69: #if defined(PETSC_USE_COMPLEX)
70: {
71: PetscReal re = PetscRandom123Step(r123);
72: PetscReal im = PetscRandom123Step(r123);
74: if (r->iset) {
75: re = re * PetscRealPart(r->width) + PetscRealPart(r->low);
76: im = im * PetscImaginaryPart(r->width) + PetscImaginaryPart(r->low);
77: }
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: return(0);
87: }
89: PetscErrorCode PetscRandomGetValueReal_Random123(PetscRandom r,PetscReal *val)
90: {
91: PetscRandom123 *r123 = (PetscRandom123 *) r->data;
92: PetscReal rreal;
95: rreal = PetscRandom123Step(r123);
96: if (r->iset) rreal = rreal * PetscRealPart(r->width) + PetscRealPart(r->low);
97: *val = rreal;
98: return(0);
99: }
101: PetscErrorCode PetscRandomDestroy_Random123(PetscRandom r)
102: {
106: PetscFree(r->data);
107: return(0);
108: }
110: static struct _PetscRandomOps PetscRandomOps_Values = {
111: /* 0 */
112: PetscRandomSeed_Random123,
113: PetscRandomGetValue_Random123,
114: PetscRandomGetValueReal_Random123,
115: PetscRandomDestroy_Random123,
116: /* 5 */
117: 0
118: };
120: /*MC
121: PETSCRANDOM123- access to Random123 counter based pseudorandom number generators (currently threefry4x64)
123: Options Database Keys:
124: . -random_type <rand,rand48,sprng,random123>
126: Level: beginner
128: PETSc must have been ./configure with the option --download-random123 to use
129: this random number generator.
131: .seealso: RandomCreate(), RandomSetType(), PETSCRAND, PETSCRAND48, PETSCSPRNG
132: M*/
134: PETSC_EXTERN PetscErrorCode PetscRandomCreate_Random123(PetscRandom r)
135: {
136: PetscRandom123 *r123;
140: PetscNewLog(r,&r123);
141: r->data = r123;
142: PetscMemcpy(r->ops,&PetscRandomOps_Values,sizeof(PetscRandomOps_Values));
143: PetscObjectChangeTypeName((PetscObject)r,PETSCRANDOM123);
144: PetscRandomSeed(r);
145: return(0);
146: }