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 on PetscRandom

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

 30:    Level: intermediate

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

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

 49:    Not collective

 51:    Input Parameters:
 52: .  r - The random number generator context

 54:    Output Parameter:
 55: .  seed - The random seed

 57:    Level: intermediate

 59: .seealso: PetscRandomCreate(), PetscRandomSetSeed(), PetscRandomSeed()
 60: @*/
 61: PetscErrorCode  PetscRandomGetSeed(PetscRandom r,unsigned long *seed)
 62: {
 64:   if (seed) {
 66:     *seed = r->seed;
 67:   }
 68:   return 0;
 69: }

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

 74:    Not collective

 76:    Input Parameters:
 77: +  r  - The random number generator context
 78: -  seed - The random seed

 80:    Level: intermediate

 82:    Usage:
 83:       PetscRandomSetSeed(r,a positive integer);
 84:       PetscRandomSeed(r);  PetscRandomGetValue() will now start with the new seed.

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

 89: .seealso: PetscRandomCreate(), PetscRandomGetSeed(), PetscRandomSeed()
 90: @*/
 91: PetscErrorCode  PetscRandomSetSeed(PetscRandom r,unsigned long seed)
 92: {
 94:   r->seed = seed;
 95:   PetscInfo(NULL,"Setting seed to %d\n",(int)seed);
 96:   return 0;
 97: }

 99: /* ------------------------------------------------------------------- */
100: /*
101:   PetscRandomSetTypeFromOptions_Private - Sets the type of random generator from user options. Defaults to type PETSCRAND48 or PETSCRAND.

103:   Collective on PetscRandom

105:   Input Parameter:
106: . rnd - The random number generator context

108:   Level: intermediate

110: .seealso: PetscRandomSetFromOptions(), PetscRandomSetType()
111: */
112: static PetscErrorCode PetscRandomSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,PetscRandom rnd)
113: {
114:   PetscBool      opt;
115:   const char     *defaultType;
116:   char           typeName[256];

118:   if (((PetscObject)rnd)->type_name) {
119:     defaultType = ((PetscObject)rnd)->type_name;
120:   } else {
121:     defaultType = PETSCRANDER48;
122:   }

124:   PetscRandomRegisterAll();
125:   PetscOptionsFList("-random_type","PetscRandom type","PetscRandomSetType",PetscRandomList,defaultType,typeName,256,&opt);
126:   if (opt) {
127:     PetscRandomSetType(rnd, typeName);
128:   } else {
129:     PetscRandomSetType(rnd, defaultType);
130:   }
131:   return 0;
132: }

134: /*@
135:   PetscRandomSetFromOptions - Configures the random number generator from the options database.

137:   Collective on PetscRandom

139:   Input Parameter:
140: . rnd - The random number generator context

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

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

151:   Level: beginner

153: .seealso: PetscRandomCreate(), PetscRandomSetType()
154: @*/
155: PetscErrorCode  PetscRandomSetFromOptions(PetscRandom rnd)
156: {
158:   PetscBool      set,noimaginary = PETSC_FALSE;
159:   PetscInt       seed;


163:   PetscObjectOptionsBegin((PetscObject)rnd);

165:   /* Handle PetscRandom type options */
166:   PetscRandomSetTypeFromOptions_Private(PetscOptionsObject,rnd);

168:   /* Handle specific random generator's options */
169:   if (rnd->ops->setfromoptions) {
170:     (*rnd->ops->setfromoptions)(PetscOptionsObject,rnd);
171:   }
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 from Options

201:    Collective on PetscRandom

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 on PetscRandom

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

227:    Notes:
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:    You can change the format the vector is printed using the
236:    option PetscViewerPushFormat().

238:    Level: beginner

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

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

273:     PetscObjectGetName((PetscObject)rnd,&name);
274:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
275:     if (!((PetscObject)rnd)->amsmem && rank == 0) {
276:       char       dir[1024];

278:       PetscObjectViewSAWs((PetscObject)rnd,viewer);
279:       PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/Low",name);
280:       PetscStackCallSAWs(SAWs_Register,(dir,&rnd->low,1,SAWs_READ,SAWs_DOUBLE));
281:     }
282: #endif
283:   }
284:   return 0;
285: }

287: /*@
288:    PetscRandomCreate - Creates a context for generating random numbers,
289:    and initializes the random-number generator.

291:    Collective

293:    Input Parameters:
294: .  comm - MPI communicator

296:    Output Parameter:
297: .  r  - the random number generator context

299:    Level: intermediate

301:    Notes:
302:    The random type has to be set by PetscRandomSetType().

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

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

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

315:    Example of Usage:
316: .vb
317:       PetscRandomCreate(PETSC_COMM_SELF,&r);
318:       PetscRandomSetType(r,PETSCRAND48);
319:       PetscRandomGetValue(r,&value1);
320:       PetscRandomGetValueReal(r,&value2);
321:       PetscRandomDestroy(&r);
322: .ve

324: .seealso: PetscRandomSetType(), PetscRandomGetValue(), PetscRandomGetValueReal(), PetscRandomSetInterval(),
325:           PetscRandomDestroy(), VecSetRandom(), PetscRandomType
326: @*/

328: PetscErrorCode  PetscRandomCreate(MPI_Comm comm,PetscRandom *r)
329: {
330:   PetscRandom    rr;
331:   PetscMPIInt    rank;

334:   *r = NULL;
335:   PetscRandomInitializePackage();

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

339:   MPI_Comm_rank(comm,&rank);

341:   rr->data  = NULL;
342:   rr->low   = 0.0;
343:   rr->width = 1.0;
344:   rr->iset  = PETSC_FALSE;
345:   rr->seed  = 0x12345678 + 76543*rank;
346:   PetscRandomSetType(rr,PETSCRANDER48);
347:   *r = rr;
348:   return 0;
349: }

351: /*@
352:    PetscRandomSeed - Seed the generator.

354:    Not collective

356:    Input Parameters:
357: .  r - The random number generator context

359:    Level: intermediate

361:    Usage:
362:       PetscRandomSetSeed(r,a positive integer);
363:       PetscRandomSeed(r);  PetscRandomGetValue() will now start with the new seed.

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

368: .seealso: PetscRandomCreate(), PetscRandomGetSeed(), PetscRandomSetSeed()
369: @*/
370: PetscErrorCode  PetscRandomSeed(PetscRandom r)
371: {

375:   (*r->ops->seed)(r);
376:   PetscObjectStateIncrease((PetscObject)r);
377:   return 0;
378: }