Actual source code: random.c

      1: #define PETSC_DLL
      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.h
     16:  #include petscsys.h
     17: #if defined (PETSC_HAVE_STDLIB_H)
     18: #include <stdlib.h>
     19: #else
     20: /* maybe the protypes are missing */
     21: #if defined(PETSC_HAVE_DRAND48)
     26: #else
     28: #endif
     29: #endif

     31: PetscCookie PETSC_DLLEXPORT PETSC_RANDOM_COOKIE = 0;

     33: /* Private data */
     34: struct _p_PetscRandom {
     35:   PETSCHEADER(int);
     36:   unsigned    long seed;
     37:   PetscScalar low,width;       /* lower bound and width of the interval over
     38:                                   which the random numbers are distributed */
     39:   PetscTruth  iset;             /* if true, indicates that the user has set the interval */
     40:   /* array for shuffling ??? */
     41: };

     45: /*@C
     46:    PetscRandomDestroy - Destroys a context that has been formed by
     47:    PetscRandomCreate().

     49:    Collective on PetscRandom

     51:    Intput Parameter:
     52: .  r  - the random number generator context

     54:    Level: intermediate

     56: .seealso: PetscRandomGetValue(), PetscRandomCreate(), VecSetRandom()
     57: @*/
     58: PetscErrorCode PETSC_DLLEXPORT PetscRandomDestroy(PetscRandom r)
     59: {
     63:   if (--r->refct > 0) return(0);
     64:   PetscHeaderDestroy(r);
     65:   return(0);
     66: }

     70: /*@C
     71:    PetscRandomSetInterval - Sets the interval over which the random numbers
     72:    will be randomly distributed.  By default, this interval is [0,1).

     74:    Collective on PetscRandom

     76:    Input Parameters:
     77: .  r  - the random number generator context

     79:    Example of Usage:
     80: .vb
     81:       PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&r);
     82:       PetscRandomSetInterval(RANDOM_DEFAULT,&r);
     83:       PetscRandomGetValue(r,&value1);
     84:       PetscRandomGetValue(r,&value2);
     85:       PetscRandomDestroy(r);
     86: .ve

     88:    Level: intermediate

     90:    Concepts: random numbers^range

     92: .seealso: PetscRandomCreate()
     93: @*/
     94: PetscErrorCode PETSC_DLLEXPORT PetscRandomSetInterval(PetscRandom r,PetscScalar low,PetscScalar high)
     95: {
     98: #if defined(PETSC_USE_COMPLEX)
     99:   if (PetscRealPart(low) >= PetscRealPart(high))           SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"only low < high");
    100:   if (PetscImaginaryPart(low) >= PetscImaginaryPart(high)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"only low < high");
    101: #else
    102:   if (low >= high) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"only low < high: Instead %g %g",low,high);
    103: #endif
    104:   r->low   = low;
    105:   r->width = high-low;
    106:   r->iset  = PETSC_TRUE;
    107:   return(0);
    108: }

    110: /*
    111:    For now we have set up using the DRAND48() generater. We need to deal
    112:    with other variants of random number generators. We should also add
    113:    a routine to enable restarts [seed48()]
    114: */
    115: #if defined(PETSC_HAVE_DRAND48)

    119: /*@C
    120:    PetscRandomCreate - Creates a context for generating random numbers,
    121:    and initializes the random-number generator.

    123:    Collective on MPI_Comm

    125:    Input Parameters:
    126: +  comm - MPI communicator
    127: -  type - the type of random numbers to be generated, usually RANDOM_DEFAULT

    129:    Output Parameter:
    130: .  r  - the random number generator context

    132:    Level: intermediate

    134:    Notes:
    135:    By default, we generate random numbers via srand48()/drand48() that
    136:    are uniformly distributed over [0,1).  The user can shift and stretch
    137:    this interval by calling PetscRandomSetInterval().
    138:
    139:    Currently three types of random numbers are supported. These types
    140:    are equivalent when working with real numbers.
    141: .     RANDOM_DEFAULT - both real and imaginary components are random
    142: .     RANDOM_DEFAULT_REAL - real component is random; imaginary component is 0
    143: .     RANDOM_DEFAULT_IMAGINARY - imaginary component is random; real component is 0

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

    147:    Example of Usage:
    148: .vb
    149:       PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&r);
    150:       PetscRandomGetValue(r,&value1);
    151:       PetscRandomGetValue(r,&value2);
    152:       PetscRandomGetValue(r,&value3);
    153:       PetscRandomDestroy(r);
    154: .ve

    156:    Concepts: random numbers^creating

    158: .seealso: PetscRandomGetValue(), PetscRandomSetInterval(), PetscRandomDestroy(), VecSetRandom()
    159: @*/
    160: PetscErrorCode PETSC_DLLEXPORT PetscRandomCreate(MPI_Comm comm,PetscRandomType type,PetscRandom *r)
    161: {
    162:   PetscRandom    rr;
    164:   PetscMPIInt    rank;

    167:   *r = 0;
    168:   if (type != RANDOM_DEFAULT && type != RANDOM_DEFAULT_REAL && type != RANDOM_DEFAULT_IMAGINARY){
    169:     SETERRQ(PETSC_ERR_SUP,"Not for this random number type");
    170:   }
    171:   PetscHeaderCreate(rr,_p_PetscRandom,int,PETSC_RANDOM_COOKIE,type,"random",comm,PetscRandomDestroy,0);
    172:   rr->low   = 0.0;
    173:   rr->width = 1.0;
    174:   rr->iset  = PETSC_FALSE;
    175:   rr->seed  = 0;
    176:   MPI_Comm_rank(comm,&rank);
    177:   srand48(0x12345678+rank);
    178:   *r = rr;
    179:   return(0);
    180: }

    184: /*@C
    185:    PetscRandomGetValue - Generates a random number.  Call this after first calling
    186:    PetscRandomCreate().

    188:    Not Collective

    190:    Intput Parameter:
    191: .  r  - the random number generator context

    193:    Output Parameter:
    194: .  val - the value

    196:    Level: intermediate

    198:    Notes:
    199:    Use VecSetRandom() to set the elements of a vector to random numbers.

    201:    Example of Usage:
    202: .vb
    203:       PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&r);
    204:       PetscRandomGetValue(r,&value1);
    205:       PetscRandomGetValue(r,&value2);
    206:       PetscRandomGetValue(r,&value3);
    207:       PetscRandomDestroy(r);
    208: .ve

    210:    Concepts: random numbers^getting

    212: .seealso: PetscRandomCreate(), PetscRandomDestroy(), VecSetRandom()
    213: @*/
    214: PetscErrorCode PETSC_DLLEXPORT PetscRandomGetValue(PetscRandom r,PetscScalar *val)
    215: {
    219: #if defined(PETSC_USE_COMPLEX)
    220:   if (r->type == RANDOM_DEFAULT) {
    221:     if (r->iset) {
    222:          *val = PetscRealPart(r->width)*drand48() + PetscRealPart(r->low) +
    223:                 (PetscImaginaryPart(r->width)*drand48() + PetscImaginaryPart(r->low)) * PETSC_i;
    224:     }
    225:     else *val = drand48() + drand48()*PETSC_i;
    226:   } else if (r->type == RANDOM_DEFAULT_REAL) {
    227:     if (r->iset) *val = PetscRealPart(r->width)*drand48() + PetscRealPart(r->low);
    228:     else                       *val = drand48();
    229:   } else if (r->type == RANDOM_DEFAULT_IMAGINARY) {
    230:     if (r->iset) *val = (PetscImaginaryPart(r->width)*drand48()+PetscImaginaryPart(r->low))*PETSC_i;
    231:     else         *val = drand48()*PETSC_i;
    232:   } else {
    233:     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid random number type");
    234:   }
    235: #else
    236:   if (r->iset) *val = r->width * drand48() + r->low;
    237:   else         *val = drand48();
    238: #endif
    239:   return(0);
    240: }

    242: #elif defined(PETSC_HAVE_RAND)

    246: PetscErrorCode PETSC_DLLEXPORT PetscRandomCreate(MPI_Comm comm,PetscRandomType type,PetscRandom *r)
    247: {
    248:   PetscRandom    rr;
    250:   PetscMPIInt    rank;

    253:   PetscLogInfo((0,"PetscRandomCreate: using rand(). not as efficinet as dran48\n"));
    254:   *r = 0;
    255:   if (type != RANDOM_DEFAULT && type != RANDOM_DEFAULT_REAL && type != RANDOM_DEFAULT_IMAGINARY) {
    256:     SETERRQ(PETSC_ERR_SUP,"Not for this random number type");
    257:   }
    258:   PetscHeaderCreate(rr,_p_PetscRandom,int,PETSC_RANDOM_COOKIE,type,"random",comm,PetscRandomDestroy,0);
    259:   rr->low   = 0.0;
    260:   rr->width = 1.0;
    261:   rr->iset  = PETSC_FALSE;
    262:   rr->seed  = 0;
    263:   MPI_Comm_rank(comm,&rank);
    264:   srand(0x12345678+rank);
    265:   *r = rr;
    266:   return(0);
    267: }

    269: #define RAND_WRAP() (rand()/(double)((unsigned int)RAND_MAX+1))
    272: PetscErrorCode PETSC_DLLEXPORT PetscRandomGetValue(PetscRandom r,PetscScalar *val)
    273: {
    277: #if defined(PETSC_USE_COMPLEX)
    278:   if (r->type == RANDOM_DEFAULT) {
    279:     if (r->iset)
    280:          *val = PetscRealPart(r->width)*RAND_WRAP() + PetscRealPart(r->low) +
    281:                 (PetscImaginaryPart(r->width)*RAND_WRAP() + PetscImaginaryPart(r->low)) * PETSC_i;
    282:     else *val = RAND_WRAP() + RAND_WRAP()*PETSC_i;
    283:   } else if (r->type == RANDOM_DEFAULT_REAL) {
    284:     if (r->iset) *val = PetscRealPart(r->width)*RAND_WRAP() + PetscRealPart(r->low);
    285:     else         *val = RAND_WRAP();
    286:   } else if (r->type == RANDOM_DEFAULT_IMAGINARY) {
    287:     if (r->iset) *val = (PetscImaginaryPart(r->width)*RAND_WRAP()+PetscImaginaryPart(r->low))*PETSC_i;
    288:     else         *val = RAND_WRAP()*PETSC_i;
    289:   } else SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid random number type");
    290: #else
    291:   if (r->iset) *val = r->width * RAND_WRAP() + r->low;
    292:   else         *val = RAND_WRAP();
    293: #endif
    294:   return(0);
    295: }

    297: #else
    298: /* Should put a simple, portable random number generator here? */

    302: PetscErrorCode PETSC_DLLEXPORT PetscRandomCreate(MPI_Comm comm,PetscRandomType type,PetscRandom *r)
    303: {
    304:   PetscRandom    rr;
    305:   char           arch[10];

    309:   *r = 0;
    310:   if (type != RANDOM_DEFAULT) SETERRQ(PETSC_ERR_SUP,"Not for this random number type");
    311:   PetscHeaderCreate(rr,_p_PetscRandom,int,PETSC_RANDOM_COOKIE,type,"random",comm,PetscRandomDestroy,0);
    312:   *r = rr;
    313:   PetscGetArchType(arch,10);
    314:   PetscPrintf(comm,"PetscRandomCreate: Warning: Random number generator not set for machine %s; using fake random numbers.\n",arch);
    315:   return(0);
    316: }

    320: PetscErrorCode PETSC_DLLEXPORT PetscRandomGetValue(PetscRandom r,PetscScalar *val)
    321: {
    325: #if defined(PETSC_USE_COMPLEX)
    326:   *val = (0.5,0.5);
    327: #else
    328:   *val = 0.5;
    329: #endif
    330:   return(0);
    331: }
    332: #endif