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 <../src/sys/classes/random/randomimpl.h>
 16: #include <petscviewer.h>

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

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

 25:    Collective on PetscRandom

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

 30:    Level: intermediate

 32: .seealso: PetscRandomGetValue(), PetscRandomCreate(), VecSetRandom()
 33: @*/
 34: PetscErrorCode  PetscRandomDestroy(PetscRandom *r)
 35: {

 39:   if (!*r) return(0);
 41:   if (--((PetscObject)(*r))->refct > 0) {*r = NULL; return(0);}
 42:   if ((*r)->ops->destroy) {
 43:     (*(*r)->ops->destroy)(*r);
 44:   }
 45:   PetscHeaderDestroy(r);
 46:   return(0);
 47: }


 50: /*@C
 51:    PetscRandomGetSeed - Gets the random seed.

 53:    Not collective

 55:    Input Parameters:
 56: .  r - The random number generator context

 58:    Output Parameter:
 59: .  seed - The random seed

 61:    Level: intermediate

 63: .seealso: PetscRandomCreate(), PetscRandomSetSeed(), PetscRandomSeed()
 64: @*/
 65: PetscErrorCode  PetscRandomGetSeed(PetscRandom r,unsigned long *seed)
 66: {
 69:   if (seed) {
 71:     *seed = r->seed;
 72:   }
 73:   return(0);
 74: }

 76: /*@C
 77:    PetscRandomSetSeed - Sets the random seed. You MUST call PetscRandomSeed() after this call to have the new seed used.

 79:    Not collective

 81:    Input Parameters:
 82: +  r  - The random number generator context
 83: -  seed - The random seed

 85:    Level: intermediate

 87:    Usage:
 88:       PetscRandomSetSeed(r,a positive integer);
 89:       PetscRandomSeed(r);  PetscRandomGetValue() will now start with the new seed.

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

 94: .seealso: PetscRandomCreate(), PetscRandomGetSeed(), PetscRandomSeed()
 95: @*/
 96: PetscErrorCode  PetscRandomSetSeed(PetscRandom r,unsigned long seed)
 97: {

102:   r->seed = seed;
103:   PetscInfo1(NULL,"Setting seed to %d\n",(int)seed);
104:   return(0);
105: }

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

111:   Collective on PetscRandom

113:   Input Parameter:
114: . rnd - The random number generator context

116:   Level: intermediate

118: .seealso: PetscRandomSetFromOptions(), PetscRandomSetType()
119: */
120: static PetscErrorCode PetscRandomSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,PetscRandom rnd)
121: {
122:   PetscBool      opt;
123:   const char     *defaultType;
124:   char           typeName[256];

128:   if (((PetscObject)rnd)->type_name) {
129:     defaultType = ((PetscObject)rnd)->type_name;
130:   } else {
131:     defaultType = PETSCRANDER48;
132:   }

134:   PetscRandomRegisterAll();
135:   PetscOptionsFList("-random_type","PetscRandom type","PetscRandomSetType",PetscRandomList,defaultType,typeName,256,&opt);
136:   if (opt) {
137:     PetscRandomSetType(rnd, typeName);
138:   } else {
139:     PetscRandomSetType(rnd, defaultType);
140:   }
141:   return(0);
142: }

144: /*@
145:   PetscRandomSetFromOptions - Configures the random number generator from the options database.

147:   Collective on PetscRandom

149:   Input Parameter:
150: . rnd - The random number generator context

152:   Options Database:
153: + -random_seed <integer> - provide a seed to the random number generater
154: - -random_no_imaginary_part - makes the imaginary part of the random number zero, this is useful when you want the
155:                               same code to produce the same result when run with real numbers or complex numbers for regression testing purposes

157:   Notes:
158:     To see all options, run your program with the -help option.
159:           Must be called after PetscRandomCreate() but before the rnd is used.

161:   Level: beginner

163: .seealso: PetscRandomCreate(), PetscRandomSetType()
164: @*/
165: PetscErrorCode  PetscRandomSetFromOptions(PetscRandom rnd)
166: {
168:   PetscBool      set,noimaginary = PETSC_FALSE;
169:   PetscInt       seed;


174:   PetscObjectOptionsBegin((PetscObject)rnd);

176:   /* Handle PetscRandom type options */
177:   PetscRandomSetTypeFromOptions_Private(PetscOptionsObject,rnd);

179:   /* Handle specific random generator's options */
180:   if (rnd->ops->setfromoptions) {
181:     (*rnd->ops->setfromoptions)(PetscOptionsObject,rnd);
182:   }
183:   PetscOptionsInt("-random_seed","Seed to use to generate random numbers","PetscRandomSetSeed",0,&seed,&set);
184:   if (set) {
185:     PetscRandomSetSeed(rnd,(unsigned long int)seed);
186:     PetscRandomSeed(rnd);
187:   }
188:   PetscOptionsBool("-random_no_imaginary_part","The imaginary part of the random number will be zero","PetscRandomSetInterval",noimaginary,&noimaginary,&set);
189: #if defined(PETSC_HAVE_COMPLEX)
190:   if (set) {
191:     if (noimaginary) {
192:       PetscScalar low,high;
193:       PetscRandomGetInterval(rnd,&low,&high);
194:       low  = low - PetscImaginaryPart(low);
195:       high = high - PetscImaginaryPart(high);
196:       PetscRandomSetInterval(rnd,low,high);
197:     }
198:   }
199: #endif
200:   PetscOptionsEnd();
201:   PetscRandomViewFromOptions(rnd,NULL, "-random_view");
202:   return(0);
203: }

205: #if defined(PETSC_HAVE_SAWS)
206: #include <petscviewersaws.h>
207: #endif

209: /*@C
210:    PetscRandomViewFromOptions - View from Options

212:    Collective on PetscRandom

214:    Input Parameters:
215: +  A - the  random number generator context
216: .  obj - Optional object
217: -  name - command line option

219:    Level: intermediate
220: .seealso:  PetscRandom, PetscRandomView, PetscObjectViewFromOptions(), PetscRandomCreate()
221: @*/
222: PetscErrorCode  PetscRandomViewFromOptions(PetscRandom A,PetscObject obj,const char name[])
223: {

228:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
229:   return(0);
230: }

232: /*@C
233:    PetscRandomView - Views a random number generator object.

235:    Collective on PetscRandom

237:    Input Parameters:
238: +  rnd - The random number generator context
239: -  viewer - an optional visualization context

241:    Notes:
242:    The available visualization contexts include
243: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
244: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
245:          output where only the first processor opens
246:          the file.  All other processors send their
247:          data to the first processor to print.

249:    You can change the format the vector is printed using the
250:    option PetscViewerPushFormat().

252:    Level: beginner

254: .seealso:  PetscRealView(), PetscScalarView(), PetscIntView()
255: @*/
256: PetscErrorCode  PetscRandomView(PetscRandom rnd,PetscViewer viewer)
257: {
259:   PetscBool      iascii;
260: #if defined(PETSC_HAVE_SAWS)
261:   PetscBool      issaws;
262: #endif

267:   if (!viewer) {
268:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rnd),&viewer);
269:   }
272:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
273: #if defined(PETSC_HAVE_SAWS)
274:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
275: #endif
276:   if (iascii) {
277:     PetscMPIInt rank;
278:     PetscObjectPrintClassNamePrefixType((PetscObject)rnd,viewer);
279:     MPI_Comm_rank(PetscObjectComm((PetscObject)rnd),&rank);
280:     PetscViewerASCIIPushSynchronized(viewer);
281:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Random type %s, seed %lu\n",rank,((PetscObject)rnd)->type_name,rnd->seed);
282:     PetscViewerFlush(viewer);
283:     PetscViewerASCIIPopSynchronized(viewer);
284: #if defined(PETSC_HAVE_SAWS)
285:   } else if (issaws) {
286:     PetscMPIInt rank;
287:     const char  *name;

289:     PetscObjectGetName((PetscObject)rnd,&name);
290:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
291:     if (!((PetscObject)rnd)->amsmem && !rank) {
292:       char       dir[1024];

294:       PetscObjectViewSAWs((PetscObject)rnd,viewer);
295:       PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/Low",name);
296:       PetscStackCallSAWs(SAWs_Register,(dir,&rnd->low,1,SAWs_READ,SAWs_DOUBLE));
297:     }
298: #endif
299:   }
300:   return(0);
301: }

303: /*@
304:    PetscRandomCreate - Creates a context for generating random numbers,
305:    and initializes the random-number generator.

307:    Collective

309:    Input Parameters:
310: .  comm - MPI communicator

312:    Output Parameter:
313: .  r  - the random number generator context

315:    Level: intermediate

317:    Notes:
318:    The random type has to be set by PetscRandomSetType().

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

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

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

331:    Example of Usage:
332: .vb
333:       PetscRandomCreate(PETSC_COMM_SELF,&r);
334:       PetscRandomSetType(r,PETSCRAND48);
335:       PetscRandomGetValue(r,&value1);
336:       PetscRandomGetValueReal(r,&value2);
337:       PetscRandomDestroy(&r);
338: .ve

340: .seealso: PetscRandomSetType(), PetscRandomGetValue(), PetscRandomGetValueReal(), PetscRandomSetInterval(),
341:           PetscRandomDestroy(), VecSetRandom(), PetscRandomType
342: @*/

344: PetscErrorCode  PetscRandomCreate(MPI_Comm comm,PetscRandom *r)
345: {
346:   PetscRandom    rr;
348:   PetscMPIInt    rank;

352:   *r = NULL;
353:   PetscRandomInitializePackage();

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

357:   MPI_Comm_rank(comm,&rank);

359:   rr->data  = NULL;
360:   rr->low   = 0.0;
361:   rr->width = 1.0;
362:   rr->iset  = PETSC_FALSE;
363:   rr->seed  = 0x12345678 + 76543*rank;
364:   PetscRandomSetType(rr,PETSCRANDER48);
365:   *r = rr;
366:   return(0);
367: }

369: /*@
370:    PetscRandomSeed - Seed the generator.

372:    Not collective

374:    Input Parameters:
375: .  r - The random number generator context

377:    Level: intermediate

379:    Usage:
380:       PetscRandomSetSeed(r,a positive integer);
381:       PetscRandomSeed(r);  PetscRandomGetValue() will now start with the new seed.

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

386: .seealso: PetscRandomCreate(), PetscRandomGetSeed(), PetscRandomSetSeed()
387: @*/
388: PetscErrorCode  PetscRandomSeed(PetscRandom r)
389: {


396:   (*r->ops->seed)(r);
397:   PetscObjectStateIncrease((PetscObject)r);
398:   return(0);
399: }