1 /*
2 This file contains routines for interfacing to random number generators.
3 This provides more than just an interface to some system random number
4 generator:
5
6 Numbers can be shuffled for use as random tuples
7
8 Multiple random number generators may be used
9
10 We are still not sure what interface we want here. There should be
11 one to reinitialize and set the seed.
12 */
13
14 #include <petsc/private/randomimpl.h> /*I "petscsys.h" I*/
15 #include <petscviewer.h>
16
17 /* Logging support */
18 PetscClassId PETSC_RANDOM_CLASSID;
19
20 /*@
21 PetscRandomDestroy - Destroys a `PetscRandom` object that was created by `PetscRandomCreate()`.
22
23 Collective
24
25 Input Parameter:
26 . r - the random number generator object
27
28 Level: intermediate
29
30 .seealso: `PetscRandom`, `PetscRandomGetValue()`, `PetscRandomCreate()`, `VecSetRandom()`
31 @*/
PetscRandomDestroy(PetscRandom * r)32 PetscErrorCode PetscRandomDestroy(PetscRandom *r)
33 {
34 PetscFunctionBegin;
35 if (!*r) PetscFunctionReturn(PETSC_SUCCESS);
36 PetscValidHeaderSpecific(*r, PETSC_RANDOM_CLASSID, 1);
37 if (--((PetscObject)*r)->refct > 0) {
38 *r = NULL;
39 PetscFunctionReturn(PETSC_SUCCESS);
40 }
41 PetscTryTypeMethod(*r, destroy);
42 PetscCall(PetscHeaderDestroy(r));
43 PetscFunctionReturn(PETSC_SUCCESS);
44 }
45
46 /*@C
47 PetscRandomGetSeed - Gets the random seed.
48
49 Not collective
50
51 Input Parameter:
52 . r - The random number generator context
53
54 Output Parameter:
55 . seed - The random seed
56
57 Level: intermediate
58
59 .seealso: `PetscRandom`, `PetscRandomCreate()`, `PetscRandomSetSeed()`, `PetscRandomSeed()`
60 @*/
PetscRandomGetSeed(PetscRandom r,PetscInt64 * seed)61 PetscErrorCode PetscRandomGetSeed(PetscRandom r, PetscInt64 *seed)
62 {
63 PetscFunctionBegin;
64 PetscValidHeaderSpecific(r, PETSC_RANDOM_CLASSID, 1);
65 if (seed) {
66 PetscAssertPointer(seed, 2);
67 *seed = (PetscInt64)r->seed;
68 }
69 PetscFunctionReturn(PETSC_SUCCESS);
70 }
71
72 /*@C
73 PetscRandomSetSeed - Sets the random seed. You MUST call `PetscRandomSeed()` after this call to have the new seed used.
74
75 Not collective
76
77 Input Parameters:
78 + r - The random number generator context
79 - seed - The random seed
80
81 Level: intermediate
82
83 Example Usage:
84 .vb
85 PetscRandomSetSeed(r,a positive integer);
86 PetscRandomSeed(r);
87 PetscRandomGetValue() will now start with the new seed.
88
89 PetscRandomSeed(r) without a call to PetscRandomSetSeed() re-initializes
90 the seed. The random numbers generated will be the same as before.
91 .ve
92
93 .seealso: `PetscRandom`, `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSeed()`
94 @*/
PetscRandomSetSeed(PetscRandom r,PetscInt64 seed)95 PetscErrorCode PetscRandomSetSeed(PetscRandom r, PetscInt64 seed)
96 {
97 PetscFunctionBegin;
98 PetscValidHeaderSpecific(r, PETSC_RANDOM_CLASSID, 1);
99 r->seed = (unsigned long)seed;
100 PetscCall(PetscInfo(NULL, "Setting seed to %d\n", (int)seed));
101 PetscFunctionReturn(PETSC_SUCCESS);
102 }
103
104 /*
105 PetscRandomSetTypeFromOptions_Private - Sets the type of random generator from user options. Defaults to type PETSCRAND48 or PETSCRAND.
106
107 Collective
108
109 Input Parameter:
110 . rnd - The random number generator context
111
112 Level: intermediate
113
114 .seealso: `PetscRandomSetFromOptions()`, `PetscRandomSetType()`
115 */
PetscRandomSetTypeFromOptions_Private(PetscRandom rnd,PetscOptionItems PetscOptionsObject)116 static PetscErrorCode PetscRandomSetTypeFromOptions_Private(PetscRandom rnd, PetscOptionItems PetscOptionsObject)
117 {
118 PetscBool opt;
119 const char *defaultType;
120 char typeName[256];
121
122 PetscFunctionBegin;
123 if (((PetscObject)rnd)->type_name) {
124 defaultType = ((PetscObject)rnd)->type_name;
125 } else {
126 defaultType = PETSCRANDER48;
127 }
128
129 PetscCall(PetscRandomRegisterAll());
130 PetscCall(PetscOptionsFList("-random_type", "PetscRandom type", "PetscRandomSetType", PetscRandomList, defaultType, typeName, 256, &opt));
131 if (opt) {
132 PetscCall(PetscRandomSetType(rnd, typeName));
133 } else {
134 PetscCall(PetscRandomSetType(rnd, defaultType));
135 }
136 PetscFunctionReturn(PETSC_SUCCESS);
137 }
138
139 /*@
140 PetscRandomSetFromOptions - Configures the random number generator from the options database.
141
142 Collective
143
144 Input Parameter:
145 . rnd - The random number generator context
146
147 Options Database Keys:
148 + -random_seed <integer> - provide a seed to the random number generator
149 - -random_no_imaginary_part - makes the imaginary part of the random number zero, this is useful when you want the
150 same code to produce the same result when run with real numbers or complex numbers for regression testing purposes
151
152 Level: beginner
153
154 Note:
155 Must be called after `PetscRandomCreate()` but before the rnd is used.
156
157 .seealso: `PetscRandom`, `PetscRandomCreate()`, `PetscRandomSetType()`
158 @*/
PetscRandomSetFromOptions(PetscRandom rnd)159 PetscErrorCode PetscRandomSetFromOptions(PetscRandom rnd)
160 {
161 PetscBool set, noimaginary = PETSC_FALSE;
162 PetscInt seed;
163
164 PetscFunctionBegin;
165 PetscValidHeaderSpecific(rnd, PETSC_RANDOM_CLASSID, 1);
166
167 PetscObjectOptionsBegin((PetscObject)rnd);
168
169 /* Handle PetscRandom type options */
170 PetscCall(PetscRandomSetTypeFromOptions_Private(rnd, PetscOptionsObject));
171
172 /* Handle specific random generator's options */
173 PetscTryTypeMethod(rnd, setfromoptions, PetscOptionsObject);
174 PetscCall(PetscOptionsInt("-random_seed", "Seed to use to generate random numbers", "PetscRandomSetSeed", 0, &seed, &set));
175 if (set) {
176 PetscCall(PetscRandomSetSeed(rnd, (unsigned long int)seed));
177 PetscCall(PetscRandomSeed(rnd));
178 }
179 PetscCall(PetscOptionsBool("-random_no_imaginary_part", "The imaginary part of the random number will be zero", "PetscRandomSetInterval", noimaginary, &noimaginary, &set));
180 #if defined(PETSC_HAVE_COMPLEX)
181 if (set) {
182 if (noimaginary) {
183 PetscScalar low, high;
184 PetscCall(PetscRandomGetInterval(rnd, &low, &high));
185 low = low - PetscImaginaryPart(low);
186 high = high - PetscImaginaryPart(high);
187 PetscCall(PetscRandomSetInterval(rnd, low, high));
188 }
189 }
190 #endif
191 PetscOptionsEnd();
192 PetscCall(PetscRandomViewFromOptions(rnd, NULL, "-random_view"));
193 PetscFunctionReturn(PETSC_SUCCESS);
194 }
195
196 /*@
197 PetscRandomSetOptionsPrefix - Sets the prefix used for searching for all
198 `PetscRandom` options in the database.
199
200 Logically Collective
201
202 Input Parameters:
203 + r - the random number generator context
204 - prefix - the prefix to prepend to all option names
205
206 Level: advanced
207
208 Note:
209 A hyphen (-) must NOT be given at the beginning of the prefix name.
210 The first character of all runtime options is AUTOMATICALLY the hyphen.
211
212 .seealso: `PetscRandom`, `PetscRandomSetFromOptions()`
213 @*/
PetscRandomSetOptionsPrefix(PetscRandom r,const char prefix[])214 PetscErrorCode PetscRandomSetOptionsPrefix(PetscRandom r, const char prefix[])
215 {
216 PetscFunctionBegin;
217 PetscValidHeaderSpecific(r, PETSC_RANDOM_CLASSID, 1);
218 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, prefix));
219 PetscFunctionReturn(PETSC_SUCCESS);
220 }
221
222 #if defined(PETSC_HAVE_SAWS)
223 #include <petscviewersaws.h>
224 #endif
225
226 /*@
227 PetscRandomViewFromOptions - View a `PetscRandom` object based on the options database
228
229 Collective
230
231 Input Parameters:
232 + A - the random number generator context
233 . obj - Optional object
234 - name - command line option
235
236 Level: intermediate
237
238 .seealso: `PetscRandom`, `PetscRandomView`, `PetscObjectViewFromOptions()`, `PetscRandomCreate()`
239 @*/
PetscRandomViewFromOptions(PetscRandom A,PetscObject obj,const char name[])240 PetscErrorCode PetscRandomViewFromOptions(PetscRandom A, PetscObject obj, const char name[])
241 {
242 PetscFunctionBegin;
243 PetscValidHeaderSpecific(A, PETSC_RANDOM_CLASSID, 1);
244 PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
245 PetscFunctionReturn(PETSC_SUCCESS);
246 }
247
248 /*@
249 PetscRandomView - Views a random number generator object.
250
251 Collective
252
253 Input Parameters:
254 + rnd - The random number generator context
255 - viewer - an optional visualization context
256
257 Level: beginner
258
259 Note:
260 The available visualization contexts include
261 + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
262 - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
263 output where only the first processor opens
264 the file. All other processors send their
265 data to the first processor to print.
266
267 .seealso: `PetscRandom`, `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`
268 @*/
PetscRandomView(PetscRandom rnd,PetscViewer viewer)269 PetscErrorCode PetscRandomView(PetscRandom rnd, PetscViewer viewer)
270 {
271 PetscBool isascii;
272 #if defined(PETSC_HAVE_SAWS)
273 PetscBool issaws;
274 #endif
275
276 PetscFunctionBegin;
277 PetscValidHeaderSpecific(rnd, PETSC_RANDOM_CLASSID, 1);
278 PetscValidType(rnd, 1);
279 if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rnd), &viewer));
280 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
281 PetscCheckSameComm(rnd, 1, viewer, 2);
282 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
283 #if defined(PETSC_HAVE_SAWS)
284 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
285 #endif
286 if (isascii) {
287 PetscMPIInt rank;
288 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)rnd, viewer));
289 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)rnd), &rank));
290 PetscCall(PetscViewerASCIIPushSynchronized(viewer));
291 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Random type %s, seed %lu\n", rank, ((PetscObject)rnd)->type_name, rnd->seed));
292 PetscCall(PetscViewerFlush(viewer));
293 PetscCall(PetscViewerASCIIPopSynchronized(viewer));
294 #if defined(PETSC_HAVE_SAWS)
295 } else if (issaws) {
296 PetscMPIInt rank;
297 const char *name;
298
299 PetscCall(PetscObjectGetName((PetscObject)rnd, &name));
300 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
301 if (!((PetscObject)rnd)->amsmem && rank == 0) {
302 char dir[1024];
303
304 PetscCall(PetscObjectViewSAWs((PetscObject)rnd, viewer));
305 PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/Low", name));
306 PetscCallSAWs(SAWs_Register, (dir, &rnd->low, 1, SAWs_READ, SAWs_DOUBLE));
307 }
308 #endif
309 }
310 PetscFunctionReturn(PETSC_SUCCESS);
311 }
312
313 /*@
314 PetscRandomCreate - Creates an object for generating random numbers,
315 and initializes the random-number generator.
316
317 Collective
318
319 Input Parameter:
320 . comm - MPI communicator
321
322 Output Parameter:
323 . r - the random number generator object
324
325 Level: intermediate
326
327 Notes:
328 The random type has to be set by `PetscRandomSetType()`.
329
330 This is only a primitive "parallel" random number generator, it should NOT
331 be used for sophisticated parallel Monte Carlo methods since it will very likely
332 not have the correct statistics across processors. You can provide your own
333 parallel generator using `PetscRandomRegister()`;
334
335 If you create a `PetscRandom()` using `PETSC_COMM_SELF` on several processors then
336 the SAME random numbers will be generated on all those processors. Use `PETSC_COMM_WORLD`
337 or the appropriate parallel communicator to eliminate this issue.
338
339 Use `VecSetRandom()` to set the elements of a vector to random numbers.
340
341 Example of Usage:
342 .vb
343 PetscRandomCreate(PETSC_COMM_SELF,&r);
344 PetscRandomSetType(r,PETSCRAND48);
345 PetscRandomGetValue(r,&value1);
346 PetscRandomGetValueReal(r,&value2);
347 PetscRandomDestroy(&r);
348 .ve
349
350 .seealso: `PetscRandomSetType()`, `PetscRandomGetValue()`, `PetscRandomGetValueReal()`, `PetscRandomSetInterval()`,
351 `PetscRandomDestroy()`, `VecSetRandom()`, `PetscRandomType`, `PetscRandom`
352 @*/
PetscRandomCreate(MPI_Comm comm,PetscRandom * r)353 PetscErrorCode PetscRandomCreate(MPI_Comm comm, PetscRandom *r)
354 {
355 PetscRandom rr;
356 PetscMPIInt rank;
357
358 PetscFunctionBegin;
359 PetscAssertPointer(r, 2);
360 PetscCall(PetscRandomInitializePackage());
361
362 PetscCall(PetscHeaderCreate(rr, PETSC_RANDOM_CLASSID, "PetscRandom", "Random number generator", "Sys", comm, PetscRandomDestroy, PetscRandomView));
363 PetscCallMPI(MPI_Comm_rank(comm, &rank));
364 rr->data = NULL;
365 rr->low = 0.0;
366 rr->width = 1.0;
367 rr->iset = PETSC_FALSE;
368 rr->seed = 0x12345678 + 76543 * rank;
369 PetscCall(PetscRandomSetType(rr, PETSCRANDER48));
370 *r = rr;
371 PetscFunctionReturn(PETSC_SUCCESS);
372 }
373
374 /*@
375 PetscRandomSeed - Seed the random number generator.
376
377 Not collective
378
379 Input Parameter:
380 . r - The random number generator context
381
382 Level: intermediate
383
384 Example Usage:
385 .vb
386 PetscRandomSetSeed(r,a positive integer);
387 PetscRandomSeed(r);
388 PetscRandomGetValue() will now start with the new seed.
389
390 PetscRandomSeed(r) without a call to PetscRandomSetSeed() re-initializes
391 the seed. The random numbers generated will be the same as before.
392 .ve
393
394 .seealso: `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSetSeed()`
395 @*/
PetscRandomSeed(PetscRandom r)396 PetscErrorCode PetscRandomSeed(PetscRandom r)
397 {
398 PetscFunctionBegin;
399 PetscValidHeaderSpecific(r, PETSC_RANDOM_CLASSID, 1);
400 PetscValidType(r, 1);
401
402 PetscUseTypeMethod(r, seed);
403 PetscCall(PetscObjectStateIncrease((PetscObject)r));
404 PetscFunctionReturn(PETSC_SUCCESS);
405 }
406