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