Actual source code: random.c

  1: #define PETSC_DLL
  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 src/sys/utils/random/randomimpl.h
 16: #if defined (PETSC_HAVE_STDLIB_H)
 17: #include <stdlib.h>
 18: #endif

 20: /* Logging support */
 21: PetscCookie  PETSC_RANDOM_COOKIE = 0;

 25: /*@
 26:    PetscRandomDestroy - Destroys a context that has been formed by 
 27:    PetscRandomCreate().

 29:    Collective on PetscRandom

 31:    Intput Parameter:
 32: .  r  - the random number generator context

 34:    Level: intermediate

 36: .seealso: PetscRandomGetValue(), PetscRandomCreate(), VecSetRandom()
 37: @*/
 38: PetscErrorCode  PetscRandomDestroy(PetscRandom r)
 39: {
 43:   if (--r->refct > 0) return(0);
 44:   PetscHeaderDestroy(r);
 45:   return(0);
 46: }

 50: /*@C
 51:    PetscRandomGetInterval - Gets the interval over which the random numbers
 52:    will be randomly distributed.  By default, this interval is [0,1).

 54:    Not collective

 56:    Input Parameters:
 57: .  r  - the random number generator context

 59:    Output Parameters:
 60: +  low - The lower bound of the interval
 61: -  high - The upper bound of the interval

 63:    Level: intermediate

 65:    Concepts: random numbers^range

 67: .seealso: PetscRandomCreate(), PetscRandomSetInterval()
 68: @*/
 69: PetscErrorCode  PetscRandomGetInterval(PetscRandom r,PetscScalar *low,PetscScalar *high)
 70: {
 73:   if (low) {
 75:     *low = r->low;
 76:   }
 77:   if (high) {
 79:     *high = r->low+r->width;
 80:   }
 81:   return(0);
 82: }

 86: /*@C
 87:    PetscRandomSetInterval - Sets the interval over which the random numbers
 88:    will be randomly distributed.  By default, this interval is [0,1).

 90:    Not collective

 92:    Input Parameters:
 93: +  r  - the random number generator context
 94: .  low - The lower bound of the interval
 95: -  high - The upper bound of the interval

 97:    Level: intermediate

 99:    Concepts: random numbers^range

101: .seealso: PetscRandomCreate(), PetscRandomGetInterval()
102: @*/
103: PetscErrorCode  PetscRandomSetInterval(PetscRandom r,PetscScalar low,PetscScalar high)
104: {
107: #if defined(PETSC_USE_COMPLEX)
108:   if (PetscRealPart(low) >= PetscRealPart(high))           SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"only low < high");
109:   if (PetscImaginaryPart(low) >= PetscImaginaryPart(high)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"only low < high");
110: #else
111:   if (low >= high) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"only low < high: Instead %G %G",low,high);
112: #endif
113:   r->low   = low;
114:   r->width = high-low;
115:   r->iset  = PETSC_TRUE;
116:   return(0);
117: }

121: /*@C
122:    PetscRandomGetSeed - Gets the random seed.

124:    Not collective

126:    Input Parameters:
127: .  r - The random number generator context

129:    Output Parameter:
130: .  seed - The random seed

132:    Level: intermediate

134:    Concepts: random numbers^seed

136: .seealso: PetscRandomCreate(), PetscRandomSetSeed(), PetscRandomSeed()
137: @*/
138: PetscErrorCode  PetscRandomGetSeed(PetscRandom r,unsigned long *seed)
139: {
142:   if (seed) {
144:     *seed = r->seed;
145:   }
146:   return(0);
147: }

151: /*@C
152:    PetscRandomSetSeed - Sets the random seed.

154:    Not collective

156:    Input Parameters:
157: +  r  - The random number generator context
158: -  seed - The random seed

160:    Level: intermediate

162:    Usage:
163:       PetscRandomSetSeed(r,a positive integer);
164:       PetscRandomSeed(r);  PetscRandomGetValue() will now start with the new seed.

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

169:    Concepts: random numbers^seed

171: .seealso: PetscRandomCreate(), PetscRandomGetSeed(), PetscRandomSeed()
172: @*/
173: PetscErrorCode  PetscRandomSetSeed(PetscRandom r,unsigned long seed)
174: {
177:   r->seed = seed;
178:   return(0);
179: }

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

187:   Collective on PetscRandom

189:   Input Parameter:
190: . rnd - The random number generator context

192:   Level: intermediate

194: .keywords: PetscRandom, set, options, database, type
195: .seealso: PetscRandomSetFromOptions(), PetscRandomSetType()
196: */
197: static PetscErrorCode PetscRandomSetTypeFromOptions_Private(PetscRandom rnd)
198: {
199:   PetscTruth     opt;
200:   const char     *defaultType;
201:   char           typeName[256];

205:   if (rnd->type_name) {
206:     defaultType = rnd->type_name;
207:   } else {
208: #if defined(PETSC_HAVE_DRAND48)    
209:     defaultType = PETSCRAND48;
210: #elif defined(PETSC_HAVE_RAND)
211:     defaultType = PETSCRAND;
212: #endif
213:   }

215:   if (!PetscRandomRegisterAllCalled) {
216:     PetscRandomRegisterAll(PETSC_NULL);
217:   }
218:   PetscOptionsList("-random_type","PetscRandom type","PetscRandomSetType",PetscRandomList,defaultType,typeName,256,&opt);
219:   if (opt) {
220:     PetscRandomSetType(rnd, typeName);
221:   } else {
222:     PetscRandomSetType(rnd, defaultType);
223:   }
224:   return(0);
225: }

229: /*@
230:   PetscRandomSetFromOptions - Configures the random number generator from the options database.

232:   Collective on PetscRandom

234:   Input Parameter:
235: . rnd - The random number generator context

237:   Notes:  To see all options, run your program with the -help option, or consult the users manual.
238:           Must be called after PetscRandomCreate() but before the rnd is used.

240:   Level: beginner

242: .keywords: PetscRandom, set, options, database
243: .seealso: PetscRandomCreate(), PetscRandomSetType()
244: @*/
245: PetscErrorCode  PetscRandomSetFromOptions(PetscRandom rnd)
246: {


252:   PetscOptionsBegin(rnd->comm, rnd->prefix, "PetscRandom options", "PetscRandom");

254:     /* Handle PetscRandom type options */
255:     PetscRandomSetTypeFromOptions_Private(rnd);

257:     /* Handle specific random generator's options */
258:     if (rnd->ops->setfromoptions) {
259:       (*rnd->ops->setfromoptions)(rnd);
260:     }
261:   PetscOptionsEnd();
262:   PetscRandomViewFromOptions(rnd, rnd->name);
263:   return(0);
264: }

268: /*@C
269:    PetscRandomView - Views a random number generator object. 

271:    Collective on PetscRandom

273:    Input Parameters:
274: +  rnd - The random number generator context
275: -  viewer - an optional visualization context

277:    Notes:
278:    The available visualization contexts include
279: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
280: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
281:          output where only the first processor opens
282:          the file.  All other processors send their 
283:          data to the first processor to print. 

285:    You can change the format the vector is printed using the 
286:    option PetscViewerSetFormat().

288:    Level: beginner

290: .seealso:  PetscRealView(), PetscScalarView(), PetscIntView()
291: @*/
292: PetscErrorCode  PetscRandomView(PetscRandom rnd,PetscViewer viewer)
293: {
294:   PetscErrorCode    ierr;
295:   PetscTruth        iascii;

300:   if (!viewer) {
301:     PetscViewerASCIIGetStdout(rnd->comm,&viewer);
302:   }
305:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
306:   if (iascii) {
307:     PetscMPIInt rank;
308:     MPI_Comm_rank(rnd->comm,&rank);
309:     PetscViewerASCIISynchronizedPrintf(viewer,"[%D] Random type %s, seed %D\n",rank,rnd->type_name,rnd->seed);
310:     PetscViewerFlush(viewer);
311:   } else {
312:     const char *tname;
313:     PetscObjectGetName((PetscObject)viewer,&tname);
314:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for this object",tname);
315:   }
316:   return(0);
317: }

319: #undef  __FUNCT__
321: /*@
322:   PetscRandomViewFromOptions - This function visualizes the type and the seed of the generated random numbers based upon user options.

324:   Collective on PetscRandom

326:   Input Parameters:
327: . rnd   - The random number generator context
328: . title - The title

330:   Level: intermediate

332: .keywords: PetscRandom, view, options, database
333: .seealso: PetscRandomSetFromOptions()
334: @*/
335: PetscErrorCode  PetscRandomViewFromOptions(PetscRandom rnd, char *title)
336: {
337:   PetscTruth     opt;
338:   PetscViewer    viewer;
339:   char           typeName[1024];
340:   char           fileName[PETSC_MAX_PATH_LEN];
341:   size_t         len;
342: 

346:   PetscOptionsHasName(rnd->prefix, "-random_view", &opt);
347:   if (opt) {
348:     PetscOptionsGetString(rnd->prefix, "-random_view", typeName, 1024, &opt);
349:     PetscStrlen(typeName, &len);
350:     if (len > 0) {
351:       PetscViewerCreate(rnd->comm, &viewer);
352:       PetscViewerSetType(viewer, typeName);
353:       PetscOptionsGetString(rnd->prefix, "-random_view_file", fileName, 1024, &opt);
354:       if (opt) {
355:         PetscViewerFileSetName(viewer, fileName);
356:       } else {
357:         PetscViewerFileSetName(viewer, rnd->name);
358:       }
359:       PetscRandomView(rnd, viewer);
360:       PetscViewerFlush(viewer);
361:       PetscViewerDestroy(viewer);
362:     } else {
363:       PetscViewer viewer;
364:       PetscViewerASCIIGetStdout(rnd->comm,&viewer);
365:       PetscRandomView(rnd, viewer);
366:     }
367:   }
368:   return(0);
369: }


374: /*@
375:    PetscRandomCreate - Creates a context for generating random numbers,
376:    and initializes the random-number generator.

378:    Collective on MPI_Comm

380:    Input Parameters:
381: +  comm - MPI communicator

383:    Output Parameter:
384: .  r  - the random number generator context

386:    Level: intermediate

388:    Notes:
389:    The random type has to be set by PetscRandomSetType().

391:    This is only a primative "parallel" random number generator, it should NOT
392:    be used for sophisticated parallel Monte Carlo methods since it will very likely
393:    not have the correct statistics across processors. You can provide your own
394:    parallel generator using PetscRandomRegister();

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

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

402:    Example of Usage:
403: .vb
404:       PetscRandomCreate(PETSC_COMM_SELF,&r);
405:       PetscRandomSetType(r,PETSCRAND48);
406:       PetscRandomGetValue(r,&value1);
407:       PetscRandomGetValueReal(r,&value2);
408:       PetscRandomGetValueImaginary(r,&value3);
409:       PetscRandomDestroy(r);
410: .ve

412:    Concepts: random numbers^creating

414: .seealso: PetscRandomSetType(), PetscRandomGetValue(), PetscRandomGetValueReal(), PetscRandomGetValueImaginary(), PetscRandomSetInterval(), PetscRandomDestroy(), VecSetRandom()
415: @*/

417: PetscErrorCode  PetscRandomCreate(MPI_Comm comm,PetscRandom *r)
418: {
419:   PetscRandom    rr;
421:   PetscMPIInt    rank;

425:   *r = PETSC_NULL;
426: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
427:   PetscRandomInitializePackage(PETSC_NULL);
428: #endif

430:   PetscHeaderCreate(rr,_p_PetscRandom,struct _PetscRandomOps,PETSC_RANDOM_COOKIE,-1,"PetscRandom",comm,PetscRandomDestroy,0);
431:   PetscLogObjectMemory(rr, sizeof(struct _p_PetscRandom));
432:   PetscMemzero(rr->ops, sizeof(struct _PetscRandomOps));
433:   rr->bops->publish = PETSC_NULL;
434:   rr->type_name     = PETSC_NULL;

436:   MPI_Comm_rank(comm,&rank);
437:   rr->data  = PETSC_NULL;
438:   rr->low   = 0.0;
439:   rr->width = 1.0;
440:   rr->iset  = PETSC_FALSE;
441:   rr->seed  = 0x12345678 + 76543*rank;
442:   *r = rr;
443:   return(0);
444: }

448: /*@C
449:    PetscRandomSeed - Seed the generator.

451:    Not collective

453:    Input Parameters:
454: .  r - The random number generator context

456:    Level: intermediate

458:    Usage:
459:       PetscRandomSetSeed(r,a positive integer);
460:       PetscRandomSeed(r);  PetscRandomGetValue() will now start with the new seed.

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

465:    Concepts: random numbers^seed

467: .seealso: PetscRandomCreate(), PetscRandomGetSeed(), PetscRandomSetSeed()
468: @*/
469: PetscErrorCode  PetscRandomSeed(PetscRandom r)
470: {


477:   (*r->ops->seed)(r);
478:   PetscObjectStateIncrease((PetscObject)r);
479:   return(0);
480: }

484: /*@
485:    PetscRandomGetValue - Generates a random number.  Call this after first calling
486:    PetscRandomCreate().

488:    Not Collective

490:    Intput Parameter:
491: .  r  - the random number generator context

493:    Output Parameter:
494: .  val - the value

496:    Level: intermediate

498:    Notes:
499:    Use VecSetRandom() to set the elements of a vector to random numbers.

501:    Example of Usage:
502: .vb
503:       PetscRandomCreate(PETSC_COMM_WORLD,&r);
504:       PetscRandomGetValue(r,&value1);
505:       PetscRandomGetValue(r,&value2);
506:       PetscRandomGetValue(r,&value3);
507:       PetscRandomDestroy(r);
508: .ve

510:    Concepts: random numbers^getting

512: .seealso: PetscRandomCreate(), PetscRandomDestroy(), VecSetRandom()
513: @*/
514: PetscErrorCode  PetscRandomGetValue(PetscRandom r,PetscScalar *val)
515: {


523:   (*r->ops->getvalue)(r,val);
524:   PetscObjectStateIncrease((PetscObject)r);
525:   return(0);
526: }

530: /*@
531:    PetscRandomGetValueReal - Generates a random number.  Call this after first calling
532:    PetscRandomCreate().

534:    Not Collective

536:    Intput Parameter:
537: .  r  - the random number generator context

539:    Output Parameter:
540: .  val - the value

542:    Level: intermediate

544:    Notes:
545:    Use VecSetRandom() to set the elements of a vector to random numbers.

547:    Example of Usage:
548: .vb
549:       PetscRandomCreate(PETSC_COMM_WORLD,&r);
550:       PetscRandomGetValueReal(r,&value1);
551:       PetscRandomGetValueReal(r,&value2);
552:       PetscRandomGetValueReal(r,&value3);
553:       PetscRandomDestroy(r);
554: .ve

556:    Concepts: random numbers^getting

558: .seealso: PetscRandomCreate(), PetscRandomDestroy(), VecSetRandom()
559: @*/
560: PetscErrorCode  PetscRandomGetValueReal(PetscRandom r,PetscReal *val)
561: {


569:   (*r->ops->getvaluereal)(r,val);
570:   PetscObjectStateIncrease((PetscObject)r);
571:   return(0);
572: }

576: /*@
577:    PetscRandomGetValueImaginary - Generates a random number.  Call this after first calling
578:    PetscRandomCreate().

580:    Not Collective

582:    Intput Parameter:
583: .  r  - the random number generator context

585:    Output Parameter:
586: .  val - the value

588:    Options Database Keys:
589: +    -random_type petscrand48
590: .    -random_type petscrand
591: -    -random_type sprng, uses SPRNG package

593:    Level: intermediate

595:    Notes:
596:    Use VecSetRandom() to set the elements of a vector to random numbers.

598:    Example of Usage:
599: .vb
600:       PetscRandomCreate(PETSC_COMM_WORLD,&r);
601:       PetscRandomGetValueImaginary(r,&value1);
602:       PetscRandomGetValueImaginary(r,&value2);
603:       PetscRandomGetValueImaginary(r,&value3);
604:       PetscRandomDestroy(r);
605: .ve

607:    Concepts: random numbers^getting

609: .seealso: PetscRandomCreate(), PetscRandomDestroy(), VecSetRandom()
610: @*/
611: PetscErrorCode  PetscRandomGetValueImaginary(PetscRandom r,PetscScalar *val)
612: {


620:   (*r->ops->getvalueimaginary)(r,val);
621:   PetscObjectStateIncrease((PetscObject)r);
622:   return(0);
623: }