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>

 16: /* Logging support */
 17: PetscClassId PETSC_RANDOM_CLASSID;

 19: /*@
 20:   PetscRandomDestroy - Destroys a `PetscRandom` object that was created by `PetscRandomCreate()`.

 22:   Collective

 24:   Input Parameter:
 25: . r - the random number generator object

 27:   Level: intermediate

 29: .seealso: `PetscRandom`, `PetscRandomGetValue()`, `PetscRandomCreate()`, `VecSetRandom()`
 30: @*/
 31: PetscErrorCode PetscRandomDestroy(PetscRandom *r)
 32: {
 33:   PetscFunctionBegin;
 34:   if (!*r) PetscFunctionReturn(PETSC_SUCCESS);
 36:   if (--((PetscObject)*r)->refct > 0) {
 37:     *r = NULL;
 38:     PetscFunctionReturn(PETSC_SUCCESS);
 39:   }
 40:   PetscTryTypeMethod(*r, destroy);
 41:   PetscCall(PetscHeaderDestroy(r));
 42:   PetscFunctionReturn(PETSC_SUCCESS);
 43: }

 45: /*@
 46:   PetscRandomGetSeed - Gets the random seed.

 48:   Not collective

 50:   Input Parameter:
 51: . r - The random number generator context

 53:   Output Parameter:
 54: . seed - The random seed

 56:   Level: intermediate

 58: .seealso: `PetscRandom`, `PetscRandomCreate()`, `PetscRandomSetSeed()`, `PetscRandomSeed()`
 59: @*/
 60: PetscErrorCode PetscRandomGetSeed(PetscRandom r, PetscInt64 *seed)
 61: {
 62:   PetscFunctionBegin;
 64:   if (seed) {
 65:     PetscAssertPointer(seed, 2);
 66:     *seed = (PetscInt64)r->seed;
 67:   }
 68:   PetscFunctionReturn(PETSC_SUCCESS);
 69: }

 71: /*@
 72:   PetscRandomSetSeed - Sets the random seed. You MUST call `PetscRandomSeed()` after this call to have the new seed used.

 74:   Not collective

 76:   Input Parameters:
 77: + r    - The random number generator context
 78: - seed - The random seed

 80:   Level: intermediate

 82:   Example Usage:
 83: .vb
 84:       PetscRandomSetSeed(r,a positive integer);
 85:       PetscRandomSeed(r);
 86:       PetscRandomGetValue() will now start with the new seed.

 88:       PetscRandomSeed(r) without a call to PetscRandomSetSeed() re-initializes
 89:       the seed. The random numbers generated will be the same as before.
 90: .ve

 92: .seealso: `PetscRandom`, `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSeed()`
 93: @*/
 94: PetscErrorCode PetscRandomSetSeed(PetscRandom r, PetscInt64 seed)
 95: {
 96:   PetscFunctionBegin;
 98:   r->seed = (unsigned long)seed;
 99:   PetscCall(PetscInfo(NULL, "Setting seed to %d\n", (int)seed));
100:   PetscFunctionReturn(PETSC_SUCCESS);
101: }

103: /*
104:   PetscRandomSetTypeFromOptions_Private - Sets the type of random generator from user options. Defaults to type PETSCRAND48 or PETSCRAND.

106:   Collective

108:   Input Parameter:
109: . rnd - The random number generator context

111:   Level: intermediate

113: .seealso: `PetscRandomSetFromOptions()`, `PetscRandomSetType()`
114: */
115: static PetscErrorCode PetscRandomSetTypeFromOptions_Private(PetscRandom rnd, PetscOptionItems PetscOptionsObject)
116: {
117:   PetscBool   opt;
118:   const char *defaultType;
119:   char        typeName[256];

121:   PetscFunctionBegin;
122:   if (((PetscObject)rnd)->type_name) {
123:     defaultType = ((PetscObject)rnd)->type_name;
124:   } else {
125:     defaultType = PETSCRANDER48;
126:   }

128:   PetscCall(PetscRandomRegisterAll());
129:   PetscCall(PetscOptionsFList("-random_type", "PetscRandom type", "PetscRandomSetType", PetscRandomList, defaultType, typeName, sizeof(typeName), &opt));
130:   if (opt) {
131:     PetscCall(PetscRandomSetType(rnd, typeName));
132:   } else {
133:     PetscCall(PetscRandomSetType(rnd, defaultType));
134:   }
135:   PetscFunctionReturn(PETSC_SUCCESS);
136: }

138: /*@
139:   PetscRandomSetFromOptions - Configures the random number generator from the options database.

141:   Collective

143:   Input Parameter:
144: . rnd - The random number generator context

146:   Options Database Keys:
147: + -random_seed seed         - provide a seed to the random number generator
148: - -random_no_imaginary_part - makes the imaginary part of the random number zero, this is useful when you want the
149:                               same code to produce the same result when run with real numbers or complex numbers for regression testing purposes

151:   Level: beginner

153:   Note:
154:   Must be called after `PetscRandomCreate()` but before the rnd is used.

156: .seealso: `PetscRandom`, `PetscRandomCreate()`, `PetscRandomSetType()`
157: @*/
158: PetscErrorCode PetscRandomSetFromOptions(PetscRandom rnd)
159: {
160:   PetscBool set, noimaginary = PETSC_FALSE;
161:   PetscInt  seed;

163:   PetscFunctionBegin;

166:   PetscObjectOptionsBegin((PetscObject)rnd);

168:   /* Handle PetscRandom type options */
169:   PetscCall(PetscRandomSetTypeFromOptions_Private(rnd, PetscOptionsObject));

171:   /* Handle specific random generator's options */
172:   PetscTryTypeMethod(rnd, setfromoptions, PetscOptionsObject);
173:   PetscCall(PetscOptionsInt("-random_seed", "Seed to use to generate random numbers", "PetscRandomSetSeed", 0, &seed, &set));
174:   if (set) {
175:     PetscCall(PetscRandomSetSeed(rnd, (unsigned long int)seed));
176:     PetscCall(PetscRandomSeed(rnd));
177:   }
178:   PetscCall(PetscOptionsBool("-random_no_imaginary_part", "The imaginary part of the random number will be zero", "PetscRandomSetInterval", noimaginary, &noimaginary, &set));
179: #if PetscDefined(HAVE_COMPLEX)
180:   if (set) {
181:     if (noimaginary) {
182:       PetscScalar low, high;
183:       PetscCall(PetscRandomGetInterval(rnd, &low, &high));
184:       low  = low - PetscImaginaryPart(low);
185:       high = high - PetscImaginaryPart(high);
186:       PetscCall(PetscRandomSetInterval(rnd, low, high));
187:     }
188:   }
189: #endif
190:   PetscOptionsEnd();
191:   PetscCall(PetscRandomViewFromOptions(rnd, NULL, "-random_view"));
192:   PetscFunctionReturn(PETSC_SUCCESS);
193: }

195: /*@
196:   PetscRandomSetOptionsPrefix - Sets the prefix used for searching for all
197:   `PetscRandom` options in the database.

199:   Logically Collective

201:   Input Parameters:
202: + r      - the random number generator context
203: - prefix - the prefix to prepend to all option names

205:   Level: advanced

207:   Note:
208:   A hyphen (-) must NOT be given at the beginning of the prefix name.
209:   The first character of all runtime options is AUTOMATICALLY the hyphen.

211: .seealso: `PetscRandom`, `PetscRandomAppendOptionsPrefix()`, `PetscRandomSetFromOptions()`
212: @*/
213: PetscErrorCode PetscRandomSetOptionsPrefix(PetscRandom r, const char prefix[])
214: {
215:   PetscFunctionBegin;
217:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, prefix));
218:   PetscFunctionReturn(PETSC_SUCCESS);
219: }

221: /*@
222:   PetscRandomAppendOptionsPrefix - Appends to the prefix used for searching for all
223:   `PetscRandom` options in the database.

225:   Logically Collective

227:   Input Parameters:
228: + r      - the random number generator context
229: - prefix - the prefix string to prepend to all `PetscRandom` option requests

231:   Level: advanced

233:   Note:
234:   A hyphen (-) must NOT be given at the beginning of the prefix name.
235:   The first character of all runtime options is AUTOMATICALLY the hyphen.

237: .seealso: `PetscRandom`, `PetscRandomSetOptionsPrefix()`, `PetscRandomSetFromOptions()`
238: @*/
239: PetscErrorCode PetscRandomAppendOptionsPrefix(PetscRandom r, const char prefix[])
240: {
241:   PetscFunctionBegin;
243:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)r, prefix));
244:   PetscFunctionReturn(PETSC_SUCCESS);
245: }

247: /*@
248:   PetscRandomGetOptionsPrefix - Gets the prefix used for searching for all
249:   `PetscRandom` options in the database.

251:   Not Collective

253:   Input Parameter:
254: . r - the random number generator context

256:   Output Parameter:
257: . prefix - pointer to the prefix string used

259:   Level: advanced

261: .seealso: `PetscRandom`, `PetscRandomSetOptionsPrefix()`, `PetscRandomAppendOptionsPrefix()`
262: @*/
263: PetscErrorCode PetscRandomGetOptionsPrefix(PetscRandom r, const char *prefix[])
264: {
265:   PetscFunctionBegin;
267:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)r, prefix));
268:   PetscFunctionReturn(PETSC_SUCCESS);
269: }

271: #if PetscDefined(HAVE_SAWS)
272: #include <petscviewersaws.h>
273: #endif

275: /*@
276:   PetscRandomViewFromOptions - View a `PetscRandom` object based on the options database

278:   Collective

280:   Input Parameters:
281: + A    - the random number generator context
282: . obj  - optional object that provides the prefix for the option names, pass `NULL` to use the options prefix of `A`
283: - name - command line option

285:   Options Database Key:
286: . -name viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`

288:   Level: intermediate

290:   Note:
291:   This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
292:   rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.

294: .seealso: `PetscRandom`, `PetscRandomView()`, `PetscOptionsCreateViewer()`, `PetscObjectViewFromOptions()`, `PetscRandomCreate()`
295: @*/
296: PetscErrorCode PetscRandomViewFromOptions(PetscRandom A, PetscObject obj, const char name[])
297: {
298:   PetscFunctionBegin;
300:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
301:   PetscFunctionReturn(PETSC_SUCCESS);
302: }

304: /*@
305:   PetscRandomView - Views a random number generator object.

307:   Collective

309:   Input Parameters:
310: + rnd    - The random number generator context
311: - viewer - an optional visualization context

313:   Level: beginner

315:   Note:
316:   The available visualization contexts include
317: +     `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
318: -     `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
319:   output where only the first processor opens
320:   the file.  All other processors send their
321:   data to the first processor to print.

323: .seealso: `PetscRandom`, `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`
324: @*/
325: PetscErrorCode PetscRandomView(PetscRandom rnd, PetscViewer viewer)
326: {
327:   PetscBool isascii;
328: #if PetscDefined(HAVE_SAWS)
329:   PetscBool issaws;
330: #endif

332:   PetscFunctionBegin;
335:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rnd), &viewer));
337:   PetscCheckSameComm(rnd, 1, viewer, 2);
338:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
339: #if PetscDefined(HAVE_SAWS)
340:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
341: #endif
342:   if (isascii) {
343:     PetscMPIInt rank;
344:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)rnd, viewer));
345:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)rnd), &rank));
346:     PetscCall(PetscViewerASCIIPushSynchronized(viewer));
347:     PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Random type %s, seed %lu\n", rank, ((PetscObject)rnd)->type_name, rnd->seed));
348:     PetscCall(PetscViewerFlush(viewer));
349:     PetscCall(PetscViewerASCIIPopSynchronized(viewer));
350: #if PetscDefined(HAVE_SAWS)
351:   } else if (issaws) {
352:     PetscMPIInt rank;
353:     const char *name;

355:     PetscCall(PetscObjectGetName((PetscObject)rnd, &name));
356:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
357:     if (!((PetscObject)rnd)->amsmem && rank == 0) {
358:       char dir[1024];

360:       PetscCall(PetscObjectViewSAWs((PetscObject)rnd, viewer));
361:       PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/Low", name));
362:       PetscCallSAWs(SAWs_Register, dir, &rnd->low, 1, SAWs_READ, SAWs_DOUBLE);
363:     }
364: #endif
365:   }
366:   PetscFunctionReturn(PETSC_SUCCESS);
367: }

369: /*@
370:   PetscRandomCreate - Creates an object for generating random numbers,
371:   and initializes the random-number generator.

373:   Collective

375:   Input Parameter:
376: . comm - MPI communicator

378:   Output Parameter:
379: . r - the random number generator object

381:   Level: intermediate

383:   Notes:
384:   The random type has to be set by `PetscRandomSetType()`.

386:   This is only a primitive "parallel" random number generator, it should NOT
387:   be used for sophisticated parallel Monte Carlo methods since it will very likely
388:   not have the correct statistics across processors. You can provide your own
389:   parallel generator using `PetscRandomRegister()`;

391:   If you create a `PetscRandom()` using `PETSC_COMM_SELF` on several processors then
392:   the SAME random numbers will be generated on all those processors. Use `PETSC_COMM_WORLD`
393:   or the appropriate parallel communicator to eliminate this issue.

395:   Use `VecSetRandom()` to set the elements of a vector to random numbers.

397:   Example of Usage:
398: .vb
399:       PetscRandomCreate(PETSC_COMM_SELF,&r);
400:       PetscRandomSetType(r,PETSCRAND48);
401:       PetscRandomGetValue(r,&value1);
402:       PetscRandomGetValueReal(r,&value2);
403:       PetscRandomDestroy(&r);
404: .ve

406: .seealso: `PetscRandomSetType()`, `PetscRandomGetValue()`, `PetscRandomGetValueReal()`, `PetscRandomSetInterval()`,
407:           `PetscRandomDestroy()`, `VecSetRandom()`, `PetscRandomType`, `PetscRandom`
408: @*/
409: PetscErrorCode PetscRandomCreate(MPI_Comm comm, PetscRandom *r)
410: {
411:   PetscRandom rr;
412:   PetscMPIInt rank;

414:   PetscFunctionBegin;
415:   PetscAssertPointer(r, 2);
416:   PetscCall(PetscRandomInitializePackage());

418:   PetscCall(PetscHeaderCreate(rr, PETSC_RANDOM_CLASSID, "PetscRandom", "Random number generator", "Sys", comm, PetscRandomDestroy, PetscRandomView));
419:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
420:   rr->data  = NULL;
421:   rr->low   = 0.0;
422:   rr->width = 1.0;
423:   rr->iset  = PETSC_FALSE;
424:   rr->seed  = 0x12345678 + 76543 * rank;
425:   PetscCall(PetscRandomSetType(rr, PETSCRANDER48));
426:   *r = rr;
427:   PetscFunctionReturn(PETSC_SUCCESS);
428: }

430: /*@
431:   PetscRandomSeed - Seed the random number generator.

433:   Not collective

435:   Input Parameter:
436: . r - The random number generator context

438:   Level: intermediate

440:   Example Usage:
441: .vb
442:       PetscRandomSetSeed(r,a positive integer);
443:       PetscRandomSeed(r);
444:       PetscRandomGetValue() will now start with the new seed.

446:       PetscRandomSeed(r) without a call to PetscRandomSetSeed() re-initializes
447:       the seed. The random numbers generated will be the same as before.
448: .ve

450: .seealso: `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSetSeed()`
451: @*/
452: PetscErrorCode PetscRandomSeed(PetscRandom r)
453: {
454:   PetscFunctionBegin;

458:   PetscUseTypeMethod(r, seed);
459:   PetscCall(PetscObjectStateIncrease((PetscObject)r));
460:   PetscFunctionReturn(PETSC_SUCCESS);
461: }