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

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

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

 52:    Not collective

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

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

 60:    Level: intermediate

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

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

 78:    Not collective

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

 84:    Level: intermediate

 86:    Usage:
 87:       PetscRandomSetSeed(r,a positive integer);
 88:       PetscRandomSeed(r);  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.

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

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

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

110:   Collective on PetscRandom

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

115:   Level: intermediate

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

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

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

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

146:   Collective on PetscRandom

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

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

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

160:   Level: beginner

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


173:   PetscObjectOptionsBegin((PetscObject)rnd);

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

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

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

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

211:    Collective on PetscRandom

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

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

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

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

234:    Collective on PetscRandom

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

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

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

251:    Level: beginner

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

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

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

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

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

306:    Collective

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

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

314:    Level: intermediate

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

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

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

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

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

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

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

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

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

356:   MPI_Comm_rank(comm,&rank);

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

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

371:    Not collective

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

376:    Level: intermediate

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

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

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


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