Actual source code: randomc.c


  2: /*
  3:     This file contains routines for interfacing to random number generators.
  4:     This provides more than just an interface to some system random number
  5:     generator:

  7:     Numbers can be shuffled for use as random tuples

  9:     Multiple random number generators may be used

 11:     We are still not sure what interface we want here.  There should be
 12:     one to reinitialize and set the seed.
 13:  */

 15: #include <petsc/private/randomimpl.h>
 16: #include <petscviewer.h>

 18: /* Logging support */
 19: PetscClassId PETSC_RANDOM_CLASSID;

 21: /*@C
 22:    PetscRandomDestroy - Destroys a context that has been formed by
 23:    `PetscRandomCreate()`.

 25:    Collective

 27:    Input Parameter:
 28: .  r  - the random number generator context

 30:    Level: intermediate

 32: .seealso: `PetscRandom`, `PetscRandomGetValue()`, `PetscRandomCreate()`, `VecSetRandom()`
 33: @*/
 34: PetscErrorCode PetscRandomDestroy(PetscRandom *r)
 35: {
 36:   if (!*r) return 0;
 38:   if (--((PetscObject)(*r))->refct > 0) {
 39:     *r = NULL;
 40:     return 0;
 41:   }
 42:   if ((*r)->ops->destroy) (*(*r)->ops->destroy)(*r);
 43:   PetscHeaderDestroy(r);
 44:   return 0;
 45: }

 47: /*@C
 48:    PetscRandomGetSeed - Gets the random seed.

 50:    Not collective

 52:    Input Parameters:
 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: {
 65:   if (seed) {
 67:     *seed = r->seed;
 68:   }
 69:   return 0;
 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:    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: {
 98:   r->seed = seed;
 99:   PetscInfo(NULL, "Setting seed to %d\n", (int)seed);
100:   return 0;
101: }

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

107:   Collective

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

112:   Level: intermediate

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

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

128:   PetscRandomRegisterAll();
129:   PetscOptionsFList("-random_type", "PetscRandom type", "PetscRandomSetType", PetscRandomList, defaultType, typeName, 256, &opt);
130:   if (opt) {
131:     PetscRandomSetType(rnd, typeName);
132:   } else {
133:     PetscRandomSetType(rnd, defaultType);
134:   }
135:   return 0;
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 <integer> - 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:   Note:
152:   Must be called after `PetscRandomCreate()` but before the rnd is used.

154:   Level: beginner

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


165:   PetscObjectOptionsBegin((PetscObject)rnd);

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

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

194: #if defined(PETSC_HAVE_SAWS)
195: #include <petscviewersaws.h>
196: #endif

198: /*@C
199:    PetscRandomViewFromOptions - View a `PetscRandom` object based on the options database

201:    Collective

203:    Input Parameters:
204: +  A - the  random number generator context
205: .  obj - Optional object
206: -  name - command line option

208:    Level: intermediate
209: .seealso: `PetscRandom`, `PetscRandomView`, `PetscObjectViewFromOptions()`, `PetscRandomCreate()`
210: @*/
211: PetscErrorCode PetscRandomViewFromOptions(PetscRandom A, PetscObject obj, const char name[])
212: {
214:   PetscObjectViewFromOptions((PetscObject)A, obj, name);
215:   return 0;
216: }

218: /*@C
219:    PetscRandomView - Views a random number generator object.

221:    Collective

223:    Input Parameters:
224: +  rnd - The random number generator context
225: -  viewer - an optional visualization context

227:    Note:
228:    The available visualization contexts include
229: +     `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
230: -     `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
231:          output where only the first processor opens
232:          the file.  All other processors send their
233:          data to the first processor to print.

235:    Level: beginner

237: .seealso: `PetscRandom`, `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`
238: @*/
239: PetscErrorCode PetscRandomView(PetscRandom rnd, PetscViewer viewer)
240: {
241:   PetscBool iascii;
242: #if defined(PETSC_HAVE_SAWS)
243:   PetscBool issaws;
244: #endif

248:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rnd), &viewer);
251:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
252: #if defined(PETSC_HAVE_SAWS)
253:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws);
254: #endif
255:   if (iascii) {
256:     PetscMPIInt rank;
257:     PetscObjectPrintClassNamePrefixType((PetscObject)rnd, viewer);
258:     MPI_Comm_rank(PetscObjectComm((PetscObject)rnd), &rank);
259:     PetscViewerASCIIPushSynchronized(viewer);
260:     PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Random type %s, seed %lu\n", rank, ((PetscObject)rnd)->type_name, rnd->seed);
261:     PetscViewerFlush(viewer);
262:     PetscViewerASCIIPopSynchronized(viewer);
263: #if defined(PETSC_HAVE_SAWS)
264:   } else if (issaws) {
265:     PetscMPIInt rank;
266:     const char *name;

268:     PetscObjectGetName((PetscObject)rnd, &name);
269:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
270:     if (!((PetscObject)rnd)->amsmem && rank == 0) {
271:       char dir[1024];

273:       PetscObjectViewSAWs((PetscObject)rnd, viewer);
274:       PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/Low", name);
275:       SAWs_Register, (dir, &rnd->low, 1, SAWs_READ, SAWs_DOUBLE);
276:     }
277: #endif
278:   }
279:   return 0;
280: }

282: /*@
283:    PetscRandomCreate - Creates a context for generating random numbers,
284:    and initializes the random-number generator.

286:    Collective

288:    Input Parameter:
289: .  comm - MPI communicator

291:    Output Parameter:
292: .  r  - the random number generator context

294:    Level: intermediate

296:    Notes:
297:    The random type has to be set by `PetscRandomSetType()`.

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

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

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

310:    Example of Usage:
311: .vb
312:       PetscRandomCreate(PETSC_COMM_SELF,&r);
313:       PetscRandomSetType(r,PETSCRAND48);
314:       PetscRandomGetValue(r,&value1);
315:       PetscRandomGetValueReal(r,&value2);
316:       PetscRandomDestroy(&r);
317: .ve

319: .seealso: `PetscRandomSetType()`, `PetscRandomGetValue()`, `PetscRandomGetValueReal()`, `PetscRandomSetInterval()`,
320:           `PetscRandomDestroy()`, `VecSetRandom()`, `PetscRandomType`, `PetscRandom`
321: @*/
322: PetscErrorCode PetscRandomCreate(MPI_Comm comm, PetscRandom *r)
323: {
324:   PetscRandom rr;
325:   PetscMPIInt rank;

328:   *r = NULL;
329:   PetscRandomInitializePackage();

331:   PetscHeaderCreate(rr, PETSC_RANDOM_CLASSID, "PetscRandom", "Random number generator", "Sys", comm, PetscRandomDestroy, PetscRandomView);

333:   MPI_Comm_rank(comm, &rank);

335:   rr->data  = NULL;
336:   rr->low   = 0.0;
337:   rr->width = 1.0;
338:   rr->iset  = PETSC_FALSE;
339:   rr->seed  = 0x12345678 + 76543 * rank;
340:   PetscRandomSetType(rr, PETSCRANDER48);
341:   *r = rr;
342:   return 0;
343: }

345: /*@
346:    PetscRandomSeed - Seed the random number generator.

348:    Not collective

350:    Input Parameters:
351: .  r - The random number generator context

353:    Level: intermediate

355:    Usage:
356: .vb
357:       PetscRandomSetSeed(r,a positive integer);
358:       PetscRandomSeed(r);
359:       PetscRandomGetValue() will now start with the new seed.

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

365: .seealso: `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSetSeed()`
366: @*/
367: PetscErrorCode PetscRandomSeed(PetscRandom r)
368: {

372:   PetscUseTypeMethod(r, seed);
373:   PetscObjectStateIncrease((PetscObject)r);
374:   return 0;
375: }