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