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