Actual source code: randomc.c

petsc-3.9.2 2018-05-20
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>
 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 = 0; 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:    Concepts: random numbers^seed

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

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

 81:    Not collective

 83:    Input Parameters:
 84: +  r  - The random number generator context
 85: -  seed - The random seed

 87:    Level: intermediate

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

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

 96:    Concepts: random numbers^seed

 98: .seealso: PetscRandomCreate(), PetscRandomGetSeed(), PetscRandomSeed()
 99: @*/
100: PetscErrorCode  PetscRandomSetSeed(PetscRandom r,unsigned long seed)
101: {

106:   r->seed = seed;
107:   PetscInfo1(NULL,"Setting seed to %d\n",(int)seed);
108:   return(0);
109: }

111: /* ------------------------------------------------------------------- */
112: /*
113:   PetscRandomSetTypeFromOptions_Private - Sets the type of random generator from user options. Defaults to type PETSCRAND48 or PETSCRAND.

115:   Collective on PetscRandom

117:   Input Parameter:
118: . rnd - The random number generator context

120:   Level: intermediate

122: .keywords: PetscRandom, set, options, database, type
123: .seealso: PetscRandomSetFromOptions(), PetscRandomSetType()
124: */
125: static PetscErrorCode PetscRandomSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,PetscRandom rnd)
126: {
127:   PetscBool      opt;
128:   const char     *defaultType;
129:   char           typeName[256];

133:   if (((PetscObject)rnd)->type_name) {
134:     defaultType = ((PetscObject)rnd)->type_name;
135:   } else {
136:     defaultType = PETSCRANDER48;
137:   }

139:   PetscRandomRegisterAll();
140:   PetscOptionsFList("-random_type","PetscRandom type","PetscRandomSetType",PetscRandomList,defaultType,typeName,256,&opt);
141:   if (opt) {
142:     PetscRandomSetType(rnd, typeName);
143:   } else {
144:     PetscRandomSetType(rnd, defaultType);
145:   }
146:   return(0);
147: }

149: /*@
150:   PetscRandomSetFromOptions - Configures the random number generator from the options database.

152:   Collective on PetscRandom

154:   Input Parameter:
155: . rnd - The random number generator context

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

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

165:   Level: beginner

167: .keywords: PetscRandom, set, options, database
168: .seealso: PetscRandomCreate(), PetscRandomSetType()
169: @*/
170: PetscErrorCode  PetscRandomSetFromOptions(PetscRandom rnd)
171: {
173:   PetscBool      set,noimaginary = PETSC_FALSE;
174:   PetscInt       seed;


179:   PetscObjectOptionsBegin((PetscObject)rnd);

181:   /* Handle PetscRandom type options */
182:   PetscRandomSetTypeFromOptions_Private(PetscOptionsObject,rnd);

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

210: #if defined(PETSC_HAVE_SAWS)
211:  #include <petscviewersaws.h>
212: #endif
213: /*@C
214:    PetscRandomView - Views a random number generator object.

216:    Collective on PetscRandom

218:    Input Parameters:
219: +  rnd - The random number generator context
220: -  viewer - an optional visualization context

222:    Notes:
223:    The available visualization contexts include
224: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
225: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
226:          output where only the first processor opens
227:          the file.  All other processors send their
228:          data to the first processor to print.

230:    You can change the format the vector is printed using the
231:    option PetscViewerPushFormat().

233:    Level: beginner

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

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

270:     PetscObjectGetName((PetscObject)rnd,&name);
271:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
272:     if (!((PetscObject)rnd)->amsmem && !rank) {
273:       char       dir[1024];

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

284: /*@
285:    PetscRandomCreate - Creates a context for generating random numbers,
286:    and initializes the random-number generator.

288:    Collective on MPI_Comm

290:    Input Parameters:
291: .  comm - MPI communicator

293:    Output Parameter:
294: .  r  - the random number generator context

296:    Level: intermediate

298:    Notes:
299:    The random type has to be set by PetscRandomSetType().

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

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

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

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

321:    Concepts: random numbers^creating

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

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

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

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

340:   MPI_Comm_rank(comm,&rank);

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

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

355:    Not collective

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

360:    Level: intermediate

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

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

369:    Concepts: random numbers^seed

371: .seealso: PetscRandomCreate(), PetscRandomGetSeed(), PetscRandomSetSeed()
372: @*/
373: PetscErrorCode  PetscRandomSeed(PetscRandom r)
374: {


381:   (*r->ops->seed)(r);
382:   PetscObjectStateIncrease((PetscObject)r);
383:   return(0);
384: }