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: static PetscErrorCode PetscRandomSeed_Random123(PetscRandom r)
20: {
21: threefry4x64_ukey_t ukey;
22: PetscRandom123 *r123 = (PetscRandom123 *)r->data;
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: }
41: static PetscReal PetscRandom123Step(PetscRandom123 *r123)
42: {
43: PetscReal scale = ((PetscReal)1.) / (((PetscReal)UINT64_MAX) + ((PetscReal)1.));
44: PetscReal shift = .5 * scale;
45: PetscInt mod = (r123->count++) % 4;
46: PetscReal ret;
48: ret = r123->result.v[mod] * scale + shift;
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: }
58: return ret;
59: }
61: static PetscErrorCode PetscRandomGetValue_Random123(PetscRandom r, PetscScalar *val)
62: {
63: PetscRandom123 *r123 = (PetscRandom123 *)r->data;
64: PetscScalar rscal;
66: PetscFunctionBegin;
67: #if defined(PETSC_USE_COMPLEX)
68: {
69: PetscReal re = PetscRandom123Step(r123);
70: PetscReal im = PetscRandom123Step(r123);
72: if (r->iset) {
73: re = re * PetscRealPart(r->width) + PetscRealPart(r->low);
74: im = im * PetscImaginaryPart(r->width) + PetscImaginaryPart(r->low);
75: }
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: }
87: static PetscErrorCode PetscRandomGetValueReal_Random123(PetscRandom r, PetscReal *val)
88: {
89: PetscRandom123 *r123 = (PetscRandom123 *)r->data;
90: PetscReal rreal;
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: }
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 = ((PetscReal)1.) / (UINT64_MAX + ((PetscReal)1.));
105: PetscReal shift = .5 * scale;
106: PetscRandom123 r123_copy;
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: }
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: }
151: static PetscErrorCode PetscRandomDestroy_Random123(PetscRandom r)
152: {
153: PetscFunctionBegin;
154: PetscCall(PetscFree(r->data));
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
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: };
169: /*MC
170: PETSCRANDOM123 - access to Random123 counter based pseudorandom number generators (currently threefry4x64)
172: Options Database Key:
173: . -random_type <rand,rand48,sprng,random123> - select the random number generator at runtim
175: Level: beginner
177: Note:
178: PETSc must be ./configure with the option --download-random123 to use this random number generator.
180: .seealso: `RandomCreate()`, `RandomSetType()`, `PETSCRAND`, `PETSCRAND48`, `PETSCSPRNG`, `PetscRandomSetFromOptions()`
181: M*/
183: PETSC_EXTERN PetscErrorCode PetscRandomCreate_Random123(PetscRandom r)
184: {
185: PetscRandom123 *r123;
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: }