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 `PetscRandom` object that was created by `PetscRandomCreate()`.
23: Collective
25: Input Parameter:
26: . r - the random number generator object
28: Level: intermediate
30: .seealso: `PetscRandom`, `PetscRandomGetValue()`, `PetscRandomCreate()`, `VecSetRandom()`
31: @*/
32: PetscErrorCode PetscRandomDestroy(PetscRandom *r)
33: {
34: PetscFunctionBegin;
35: if (!*r) PetscFunctionReturn(PETSC_SUCCESS);
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: }
46: /*@C
47: PetscRandomGetSeed - Gets the random seed.
49: Not collective
51: Input Parameter:
52: . r - The random number generator context
54: Output Parameter:
55: . seed - The random seed
57: Level: intermediate
59: .seealso: `PetscRandom`, `PetscRandomCreate()`, `PetscRandomSetSeed()`, `PetscRandomSeed()`
60: @*/
61: PetscErrorCode PetscRandomGetSeed(PetscRandom r, unsigned long *seed)
62: {
63: PetscFunctionBegin;
65: if (seed) {
66: PetscAssertPointer(seed, 2);
67: *seed = r->seed;
68: }
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
72: /*@C
73: PetscRandomSetSeed - Sets the random seed. You MUST call `PetscRandomSeed()` after this call to have the new seed used.
75: Not collective
77: Input Parameters:
78: + r - The random number generator context
79: - seed - The random seed
81: Level: intermediate
83: Example Usage:
84: .vb
85: PetscRandomSetSeed(r,a positive integer);
86: PetscRandomSeed(r);
87: PetscRandomGetValue() will now start with the new seed.
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
93: .seealso: `PetscRandom`, `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSeed()`
94: @*/
95: PetscErrorCode PetscRandomSetSeed(PetscRandom r, unsigned long seed)
96: {
97: PetscFunctionBegin;
99: r->seed = seed;
100: PetscCall(PetscInfo(NULL, "Setting seed to %d\n", (int)seed));
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: /* ------------------------------------------------------------------- */
105: /*
106: PetscRandomSetTypeFromOptions_Private - Sets the type of random generator from user options. Defaults to type PETSCRAND48 or PETSCRAND.
108: Collective
110: Input Parameter:
111: . rnd - The random number generator context
113: Level: intermediate
115: .seealso: `PetscRandomSetFromOptions()`, `PetscRandomSetType()`
116: */
117: static PetscErrorCode PetscRandomSetTypeFromOptions_Private(PetscRandom rnd, PetscOptionItems *PetscOptionsObject)
118: {
119: PetscBool opt;
120: const char *defaultType;
121: char typeName[256];
123: PetscFunctionBegin;
124: if (((PetscObject)rnd)->type_name) {
125: defaultType = ((PetscObject)rnd)->type_name;
126: } else {
127: defaultType = PETSCRANDER48;
128: }
130: PetscCall(PetscRandomRegisterAll());
131: PetscCall(PetscOptionsFList("-random_type", "PetscRandom type", "PetscRandomSetType", PetscRandomList, defaultType, typeName, 256, &opt));
132: if (opt) {
133: PetscCall(PetscRandomSetType(rnd, typeName));
134: } else {
135: PetscCall(PetscRandomSetType(rnd, defaultType));
136: }
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: /*@
141: PetscRandomSetFromOptions - Configures the random number generator from the options database.
143: Collective
145: Input Parameter:
146: . rnd - The random number generator context
148: Options Database Keys:
149: + -random_seed <integer> - provide a seed to the random number generator
150: - -random_no_imaginary_part - makes the imaginary part of the random number zero, this is useful when you want the
151: same code to produce the same result when run with real numbers or complex numbers for regression testing purposes
153: Level: beginner
155: Note:
156: Must be called after `PetscRandomCreate()` but before the rnd is used.
158: .seealso: `PetscRandom`, `PetscRandomCreate()`, `PetscRandomSetType()`
159: @*/
160: PetscErrorCode PetscRandomSetFromOptions(PetscRandom rnd)
161: {
162: PetscBool set, noimaginary = PETSC_FALSE;
163: PetscInt seed;
165: PetscFunctionBegin;
168: PetscObjectOptionsBegin((PetscObject)rnd);
170: /* Handle PetscRandom type options */
171: PetscCall(PetscRandomSetTypeFromOptions_Private(rnd, PetscOptionsObject));
173: /* Handle specific random generator's options */
174: PetscTryTypeMethod(rnd, setfromoptions, PetscOptionsObject);
175: PetscCall(PetscOptionsInt("-random_seed", "Seed to use to generate random numbers", "PetscRandomSetSeed", 0, &seed, &set));
176: if (set) {
177: PetscCall(PetscRandomSetSeed(rnd, (unsigned long int)seed));
178: PetscCall(PetscRandomSeed(rnd));
179: }
180: PetscCall(PetscOptionsBool("-random_no_imaginary_part", "The imaginary part of the random number will be zero", "PetscRandomSetInterval", noimaginary, &noimaginary, &set));
181: #if defined(PETSC_HAVE_COMPLEX)
182: if (set) {
183: if (noimaginary) {
184: PetscScalar low, high;
185: PetscCall(PetscRandomGetInterval(rnd, &low, &high));
186: low = low - PetscImaginaryPart(low);
187: high = high - PetscImaginaryPart(high);
188: PetscCall(PetscRandomSetInterval(rnd, low, high));
189: }
190: }
191: #endif
192: PetscOptionsEnd();
193: PetscCall(PetscRandomViewFromOptions(rnd, NULL, "-random_view"));
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: #if defined(PETSC_HAVE_SAWS)
198: #include <petscviewersaws.h>
199: #endif
201: /*@C
202: PetscRandomViewFromOptions - View a `PetscRandom` object based on the options database
204: Collective
206: Input Parameters:
207: + A - the random number generator context
208: . obj - Optional object
209: - name - command line option
211: Level: intermediate
213: .seealso: `PetscRandom`, `PetscRandomView`, `PetscObjectViewFromOptions()`, `PetscRandomCreate()`
214: @*/
215: PetscErrorCode PetscRandomViewFromOptions(PetscRandom A, PetscObject obj, const char name[])
216: {
217: PetscFunctionBegin;
219: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*@C
224: PetscRandomView - Views a random number generator object.
226: Collective
228: Input Parameters:
229: + rnd - The random number generator context
230: - viewer - an optional visualization context
232: Level: beginner
234: Note:
235: The available visualization contexts include
236: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
237: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
238: output where only the first processor opens
239: the file. All other processors send their
240: data to the first processor to print.
242: .seealso: `PetscRandom`, `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`
243: @*/
244: PetscErrorCode PetscRandomView(PetscRandom rnd, PetscViewer viewer)
245: {
246: PetscBool iascii;
247: #if defined(PETSC_HAVE_SAWS)
248: PetscBool issaws;
249: #endif
251: PetscFunctionBegin;
254: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rnd), &viewer));
256: PetscCheckSameComm(rnd, 1, viewer, 2);
257: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
258: #if defined(PETSC_HAVE_SAWS)
259: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
260: #endif
261: if (iascii) {
262: PetscMPIInt rank;
263: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)rnd, viewer));
264: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)rnd), &rank));
265: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
266: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Random type %s, seed %lu\n", rank, ((PetscObject)rnd)->type_name, rnd->seed));
267: PetscCall(PetscViewerFlush(viewer));
268: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
269: #if defined(PETSC_HAVE_SAWS)
270: } else if (issaws) {
271: PetscMPIInt rank;
272: const char *name;
274: PetscCall(PetscObjectGetName((PetscObject)rnd, &name));
275: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
276: if (!((PetscObject)rnd)->amsmem && rank == 0) {
277: char dir[1024];
279: PetscCall(PetscObjectViewSAWs((PetscObject)rnd, viewer));
280: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/Low", name));
281: PetscCallSAWs(SAWs_Register, (dir, &rnd->low, 1, SAWs_READ, SAWs_DOUBLE));
282: }
283: #endif
284: }
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: /*@
289: PetscRandomCreate - Creates an object for generating random numbers,
290: and initializes the random-number generator.
292: Collective
294: Input Parameter:
295: . comm - MPI communicator
297: Output Parameter:
298: . r - the random number generator object
300: Level: intermediate
302: Notes:
303: The random type has to be set by `PetscRandomSetType()`.
305: This is only a primitive "parallel" random number generator, it should NOT
306: be used for sophisticated parallel Monte Carlo methods since it will very likely
307: not have the correct statistics across processors. You can provide your own
308: parallel generator using `PetscRandomRegister()`;
310: If you create a `PetscRandom()` using `PETSC_COMM_SELF` on several processors then
311: the SAME random numbers will be generated on all those processors. Use `PETSC_COMM_WORLD`
312: or the appropriate parallel communicator to eliminate this issue.
314: Use `VecSetRandom()` to set the elements of a vector to random numbers.
316: Example of Usage:
317: .vb
318: PetscRandomCreate(PETSC_COMM_SELF,&r);
319: PetscRandomSetType(r,PETSCRAND48);
320: PetscRandomGetValue(r,&value1);
321: PetscRandomGetValueReal(r,&value2);
322: PetscRandomDestroy(&r);
323: .ve
325: .seealso: `PetscRandomSetType()`, `PetscRandomGetValue()`, `PetscRandomGetValueReal()`, `PetscRandomSetInterval()`,
326: `PetscRandomDestroy()`, `VecSetRandom()`, `PetscRandomType`, `PetscRandom`
327: @*/
328: PetscErrorCode PetscRandomCreate(MPI_Comm comm, PetscRandom *r)
329: {
330: PetscRandom rr;
331: PetscMPIInt rank;
333: PetscFunctionBegin;
334: PetscAssertPointer(r, 2);
335: *r = NULL;
336: PetscCall(PetscRandomInitializePackage());
338: PetscCall(PetscHeaderCreate(rr, PETSC_RANDOM_CLASSID, "PetscRandom", "Random number generator", "Sys", comm, PetscRandomDestroy, PetscRandomView));
340: PetscCallMPI(MPI_Comm_rank(comm, &rank));
342: rr->data = NULL;
343: rr->low = 0.0;
344: rr->width = 1.0;
345: rr->iset = PETSC_FALSE;
346: rr->seed = 0x12345678 + 76543 * rank;
347: PetscCall(PetscRandomSetType(rr, PETSCRANDER48));
348: *r = rr;
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: /*@
353: PetscRandomSeed - Seed the random number generator.
355: Not collective
357: Input Parameter:
358: . r - The random number generator context
360: Level: intermediate
362: Example Usage:
363: .vb
364: PetscRandomSetSeed(r,a positive integer);
365: PetscRandomSeed(r);
366: PetscRandomGetValue() will now start with the new seed.
368: PetscRandomSeed(r) without a call to PetscRandomSetSeed() re-initializes
369: the seed. The random numbers generated will be the same as before.
370: .ve
372: .seealso: `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSetSeed()`
373: @*/
374: PetscErrorCode PetscRandomSeed(PetscRandom r)
375: {
376: PetscFunctionBegin;
380: PetscUseTypeMethod(r, seed);
381: PetscCall(PetscObjectStateIncrease((PetscObject)r));
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }