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: /*@
 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: /*@
198:   PetscRandomSetOptionsPrefix - Sets the prefix used for searching for all
199:   `PetscRandom` options in the database.

201:   Logically Collective

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

207:   Level: advanced

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

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

223: #if defined(PETSC_HAVE_SAWS)
224: #include <petscviewersaws.h>
225: #endif

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

230:   Collective

232:   Input Parameters:
233: + A    - the random number generator context
234: . obj  - Optional object
235: - name - command line option

237:   Level: intermediate

239: .seealso: `PetscRandom`, `PetscRandomView`, `PetscObjectViewFromOptions()`, `PetscRandomCreate()`
240: @*/
241: PetscErrorCode PetscRandomViewFromOptions(PetscRandom A, PetscObject obj, const char name[])
242: {
243:   PetscFunctionBegin;
245:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
246:   PetscFunctionReturn(PETSC_SUCCESS);
247: }

249: /*@
250:   PetscRandomView - Views a random number generator object.

252:   Collective

254:   Input Parameters:
255: + rnd    - The random number generator context
256: - viewer - an optional visualization context

258:   Level: beginner

260:   Note:
261:   The available visualization contexts include
262: +     `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
263: -     `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
264:   output where only the first processor opens
265:   the file.  All other processors send their
266:   data to the first processor to print.

268: .seealso: `PetscRandom`, `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`
269: @*/
270: PetscErrorCode PetscRandomView(PetscRandom rnd, PetscViewer viewer)
271: {
272:   PetscBool iascii;
273: #if defined(PETSC_HAVE_SAWS)
274:   PetscBool issaws;
275: #endif

277:   PetscFunctionBegin;
280:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rnd), &viewer));
282:   PetscCheckSameComm(rnd, 1, viewer, 2);
283:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
284: #if defined(PETSC_HAVE_SAWS)
285:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
286: #endif
287:   if (iascii) {
288:     PetscMPIInt rank;
289:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)rnd, viewer));
290:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)rnd), &rank));
291:     PetscCall(PetscViewerASCIIPushSynchronized(viewer));
292:     PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Random type %s, seed %lu\n", rank, ((PetscObject)rnd)->type_name, rnd->seed));
293:     PetscCall(PetscViewerFlush(viewer));
294:     PetscCall(PetscViewerASCIIPopSynchronized(viewer));
295: #if defined(PETSC_HAVE_SAWS)
296:   } else if (issaws) {
297:     PetscMPIInt rank;
298:     const char *name;

300:     PetscCall(PetscObjectGetName((PetscObject)rnd, &name));
301:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
302:     if (!((PetscObject)rnd)->amsmem && rank == 0) {
303:       char dir[1024];

305:       PetscCall(PetscObjectViewSAWs((PetscObject)rnd, viewer));
306:       PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/Low", name));
307:       PetscCallSAWs(SAWs_Register, (dir, &rnd->low, 1, SAWs_READ, SAWs_DOUBLE));
308:     }
309: #endif
310:   }
311:   PetscFunctionReturn(PETSC_SUCCESS);
312: }

314: /*@
315:   PetscRandomCreate - Creates an object for generating random numbers,
316:   and initializes the random-number generator.

318:   Collective

320:   Input Parameter:
321: . comm - MPI communicator

323:   Output Parameter:
324: . r - the random number generator object

326:   Level: intermediate

328:   Notes:
329:   The random type has to be set by `PetscRandomSetType()`.

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

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

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

342:   Example of Usage:
343: .vb
344:       PetscRandomCreate(PETSC_COMM_SELF,&r);
345:       PetscRandomSetType(r,PETSCRAND48);
346:       PetscRandomGetValue(r,&value1);
347:       PetscRandomGetValueReal(r,&value2);
348:       PetscRandomDestroy(&r);
349: .ve

351: .seealso: `PetscRandomSetType()`, `PetscRandomGetValue()`, `PetscRandomGetValueReal()`, `PetscRandomSetInterval()`,
352:           `PetscRandomDestroy()`, `VecSetRandom()`, `PetscRandomType`, `PetscRandom`
353: @*/
354: PetscErrorCode PetscRandomCreate(MPI_Comm comm, PetscRandom *r)
355: {
356:   PetscRandom rr;
357:   PetscMPIInt rank;

359:   PetscFunctionBegin;
360:   PetscAssertPointer(r, 2);
361:   PetscCall(PetscRandomInitializePackage());

363:   PetscCall(PetscHeaderCreate(rr, PETSC_RANDOM_CLASSID, "PetscRandom", "Random number generator", "Sys", comm, PetscRandomDestroy, PetscRandomView));
364:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
365:   rr->data  = NULL;
366:   rr->low   = 0.0;
367:   rr->width = 1.0;
368:   rr->iset  = PETSC_FALSE;
369:   rr->seed  = 0x12345678 + 76543 * rank;
370:   PetscCall(PetscRandomSetType(rr, PETSCRANDER48));
371:   *r = rr;
372:   PetscFunctionReturn(PETSC_SUCCESS);
373: }

375: /*@
376:   PetscRandomSeed - Seed the random number generator.

378:   Not collective

380:   Input Parameter:
381: . r - The random number generator context

383:   Level: intermediate

385:   Example Usage:
386: .vb
387:       PetscRandomSetSeed(r,a positive integer);
388:       PetscRandomSeed(r);
389:       PetscRandomGetValue() will now start with the new seed.

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

395: .seealso: `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSetSeed()`
396: @*/
397: PetscErrorCode PetscRandomSeed(PetscRandom r)
398: {
399:   PetscFunctionBegin;

403:   PetscUseTypeMethod(r, seed);
404:   PetscCall(PetscObjectStateIncrease((PetscObject)r));
405:   PetscFunctionReturn(PETSC_SUCCESS);
406: }