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: }