Actual source code: randomc.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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>                              /*I "petscsys.h" I*/
 16: #include <petscviewer.h>

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

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

 27:    Collective on PetscRandom

 29:    Intput Parameter:
 30: .  r  - the random number generator context

 32:    Level: intermediate

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

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


 54: /*@
 55:    PetscRandomGetSeed - Gets the random seed.

 57:    Not collective

 59:    Input Parameters:
 60: .  r - The random number generator context

 62:    Output Parameter:
 63: .  seed - The random seed

 65:    Level: intermediate

 67:    Concepts: random numbers^seed

 69: .seealso: PetscRandomCreate(), PetscRandomSetSeed(), PetscRandomSeed()
 70: @*/
 71: PetscErrorCode  PetscRandomGetSeed(PetscRandom r,unsigned long *seed)
 72: {
 75:   if (seed) {
 77:     *seed = r->seed;
 78:   }
 79:   return(0);
 80: }

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

 87:    Not collective

 89:    Input Parameters:
 90: +  r  - The random number generator context
 91: -  seed - The random seed

 93:    Level: intermediate

 95:    Usage:
 96:       PetscRandomSetSeed(r,a positive integer);
 97:       PetscRandomSeed(r);  PetscRandomGetValue() will now start with the new seed.

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

102:    Concepts: random numbers^seed

104: .seealso: PetscRandomCreate(), PetscRandomGetSeed(), PetscRandomSeed()
105: @*/
106: PetscErrorCode  PetscRandomSetSeed(PetscRandom r,unsigned long seed)
107: {

112:   r->seed = seed;
113:   PetscInfo1(NULL,"Setting seed to %d\n",(int)seed);
114:   return(0);
115: }

117: /* ------------------------------------------------------------------- */
120: /*
121:   PetscRandomSetTypeFromOptions_Private - Sets the type of random generator from user options. Defaults to type PETSCRAND48 or PETSCRAND.

123:   Collective on PetscRandom

125:   Input Parameter:
126: . rnd - The random number generator context

128:   Level: intermediate

130: .keywords: PetscRandom, set, options, database, type
131: .seealso: PetscRandomSetFromOptions(), PetscRandomSetType()
132: */
133: static PetscErrorCode PetscRandomSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,PetscRandom rnd)
134: {
135:   PetscBool      opt;
136:   const char     *defaultType;
137:   char           typeName[256];

141:   if (((PetscObject)rnd)->type_name) {
142:     defaultType = ((PetscObject)rnd)->type_name;
143:   } else {
144:     defaultType = PETSCRANDER48;
145:   }

147:   PetscRandomRegisterAll();
148:   PetscOptionsFList("-random_type","PetscRandom type","PetscRandomSetType",PetscRandomList,defaultType,typeName,256,&opt);
149:   if (opt) {
150:     PetscRandomSetType(rnd, typeName);
151:   } else {
152:     PetscRandomSetType(rnd, defaultType);
153:   }
154:   return(0);
155: }

159: /*@
160:   PetscRandomSetFromOptions - Configures the random number generator from the options database.

162:   Collective on PetscRandom

164:   Input Parameter:
165: . rnd - The random number generator context

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

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

175:   Level: beginner

177: .keywords: PetscRandom, set, options, database
178: .seealso: PetscRandomCreate(), PetscRandomSetType()
179: @*/
180: PetscErrorCode  PetscRandomSetFromOptions(PetscRandom rnd)
181: {
183:   PetscBool      set,noimaginary = PETSC_FALSE;
184:   PetscInt       seed;


189:   PetscObjectOptionsBegin((PetscObject)rnd);

191:   /* Handle PetscRandom type options */
192:   PetscRandomSetTypeFromOptions_Private(PetscOptionsObject,rnd);

194:   /* Handle specific random generator's options */
195:   if (rnd->ops->setfromoptions) {
196:     (*rnd->ops->setfromoptions)(PetscOptionsObject,rnd);
197:   }
198:   PetscOptionsInt("-random_seed","Seed to use to generate random numbers","PetscRandomSetSeed",0,&seed,&set);
199:   if (set) {
200:     PetscRandomSetSeed(rnd,(unsigned long int)seed);
201:     PetscRandomSeed(rnd);
202:   }
203:   PetscOptionsBool("-random_no_imaginary_part","The imaginary part of the random number will be zero","PetscRandomSetInterval",noimaginary,&noimaginary,&set);
204: #if defined(PETSC_HAVE_COMPLEX)
205:   if (set) {
206:     if (noimaginary) {
207:       PetscScalar low,high;
208:       PetscRandomGetInterval(rnd,&low,&high);
209:       low  = low - PetscImaginaryPart(low);
210:       high = high - PetscImaginaryPart(high);
211:       PetscRandomSetInterval(rnd,low,high);
212:     }
213:   }
214: #endif
215:   PetscOptionsEnd();
216:   PetscRandomViewFromOptions(rnd,NULL, "-random_view");
217:   return(0);
218: }

220: #if defined(PETSC_HAVE_SAWS)
221: #include <petscviewersaws.h>
222: #endif
225: /*@C
226:    PetscRandomView - Views a random number generator object.

228:    Collective on PetscRandom

230:    Input Parameters:
231: +  rnd - The random number generator context
232: -  viewer - an optional visualization context

234:    Notes:
235:    The available visualization contexts include
236: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
237: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
238:          output where only the first processor opens
239:          the file.  All other processors send their
240:          data to the first processor to print.

242:    You can change the format the vector is printed using the
243:    option PetscViewerPushFormat().

245:    Level: beginner

247: .seealso:  PetscRealView(), PetscScalarView(), PetscIntView()
248: @*/
249: PetscErrorCode  PetscRandomView(PetscRandom rnd,PetscViewer viewer)
250: {
252:   PetscBool      iascii;
253: #if defined(PETSC_HAVE_SAWS)
254:   PetscBool      issaws;
255: #endif

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

282:     PetscObjectGetName((PetscObject)rnd,&name);
283:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
284:     if (!((PetscObject)rnd)->amsmem && !rank) {
285:       char       dir[1024];

287:       PetscObjectViewSAWs((PetscObject)rnd,viewer);
288:       PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/Low",name);
289:       PetscStackCallSAWs(SAWs_Register,(dir,&rnd->low,1,SAWs_READ,SAWs_DOUBLE));
290:     }
291: #endif
292:   }
293:   return(0);
294: }

298: /*@
299:    PetscRandomCreate - Creates a context for generating random numbers,
300:    and initializes the random-number generator.

302:    Collective on MPI_Comm

304:    Input Parameters:
305: +  comm - MPI communicator

307:    Output Parameter:
308: .  r  - the random number generator context

310:    Level: intermediate

312:    Notes:
313:    The random type has to be set by PetscRandomSetType().

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

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

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

326:    Example of Usage:
327: .vb
328:       PetscRandomCreate(PETSC_COMM_SELF,&r);
329:       PetscRandomSetType(r,PETSCRAND48);
330:       PetscRandomGetValue(r,&value1);
331:       PetscRandomGetValueReal(r,&value2);
332:       PetscRandomDestroy(&r);
333: .ve

335:    Concepts: random numbers^creating

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

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

349:   *r = NULL;
350:   PetscRandomInitializePackage();

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

354:   MPI_Comm_rank(comm,&rank);

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

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:    Concepts: random numbers^seed

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


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