Actual source code: random123.c
1: #include <petsc/private/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: threefry4x64_ctr_t counter;
9: threefry4x64_key_t key;
10: threefry4x64_ctr_t result;
11: PetscInt count;
12: } PetscRandom123;
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);
19: PetscErrorCode PetscRandomSeed_Random123(PetscRandom r)
20: {
21: threefry4x64_ukey_t ukey;
22: PetscRandom123 *r123 = (PetscRandom123 *)r->data;
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: return 0;
38: }
40: static PetscReal PetscRandom123Step(PetscRandom123 *r123)
41: {
42: PetscReal scale = ((PetscReal)1.) / (UINT64_MAX + ((PetscReal)1.));
43: PetscReal shift = .5 * scale;
44: PetscInt mod = (r123->count++) % 4;
45: PetscReal ret;
47: ret = r123->result.v[mod] * scale + shift;
49: if (mod == 3) {
50: r123->counter.v[0] += 4;
51: r123->counter.v[1] += 4;
52: r123->counter.v[2] += 4;
53: r123->counter.v[3] += 4;
54: r123->result = threefry4x64(r123->counter, r123->key);
55: }
57: return ret;
58: }
60: PetscErrorCode PetscRandomGetValue_Random123(PetscRandom r, PetscScalar *val)
61: {
62: PetscRandom123 *r123 = (PetscRandom123 *)r->data;
63: PetscScalar rscal;
65: #if defined(PETSC_USE_COMPLEX)
66: {
67: PetscReal re = PetscRandom123Step(r123);
68: PetscReal im = PetscRandom123Step(r123);
70: if (r->iset) {
71: re = re * PetscRealPart(r->width) + PetscRealPart(r->low);
72: im = im * PetscImaginaryPart(r->width) + PetscImaginaryPart(r->low);
73: }
75: rscal = PetscCMPLX(re, im);
76: }
77: #else
78: rscal = PetscRandom123Step(r123);
79: if (r->iset) rscal = rscal * r->width + r->low;
80: #endif
81: *val = rscal;
82: return 0;
83: }
85: PetscErrorCode PetscRandomGetValueReal_Random123(PetscRandom r, PetscReal *val)
86: {
87: PetscRandom123 *r123 = (PetscRandom123 *)r->data;
88: PetscReal rreal;
90: rreal = PetscRandom123Step(r123);
91: if (r->iset) rreal = rreal * PetscRealPart(r->width) + PetscRealPart(r->low);
92: *val = rreal;
93: return 0;
94: }
96: PetscErrorCode PetscRandomDestroy_Random123(PetscRandom r)
97: {
98: PetscFree(r->data);
99: return 0;
100: }
102: static struct _PetscRandomOps PetscRandomOps_Values = {
103: PetscDesignatedInitializer(seed, PetscRandomSeed_Random123),
104: PetscDesignatedInitializer(getvalue, PetscRandomGetValue_Random123),
105: PetscDesignatedInitializer(getvaluereal, PetscRandomGetValueReal_Random123),
106: PetscDesignatedInitializer(getvalues, NULL),
107: PetscDesignatedInitializer(getvaluesreal, NULL),
108: PetscDesignatedInitializer(destroy, PetscRandomDestroy_Random123),
109: };
111: /*MC
112: PETSCRANDOM123 - access to Random123 counter based pseudorandom number generators (currently threefry4x64)
114: Options Database Keys:
115: . -random_type <rand,rand48,sprng,random123> - select the random number generator at runtim
117: Note:
118: PETSc must be ./configure with the option --download-random123 to use this random number generator.
120: Level: beginner
122: .seealso: `RandomCreate()`, `RandomSetType()`, `PETSCRAND`, `PETSCRAND48`, `PETSCSPRNG`, `PetscRandomSetFromOptions()`
123: M*/
125: PETSC_EXTERN PetscErrorCode PetscRandomCreate_Random123(PetscRandom r)
126: {
127: PetscRandom123 *r123;
129: PetscNew(&r123);
130: r->data = r123;
131: PetscMemcpy(r->ops, &PetscRandomOps_Values, sizeof(PetscRandomOps_Values));
132: PetscObjectChangeTypeName((PetscObject)r, PETSCRANDOM123);
133: PetscRandomSeed(r);
134: return 0;
135: }