Actual source code: random123.c

  1: #include <petsc/private/randomimpl.h>
  2: #include <Random123/threefry.h>

  4: /* The structure of the Random123 methods are similar enough that templates could be used to make the other CBRNGs in
  5:  * the package (aes, ars, philox) available, as well as different block sizes.  But threefry4x64 is a good default,
  6:  * and I'd rather get a simple implementation up and working and come back if there's interest. */
  7: typedef struct _n_PetscRandom123 {
  8:   threefry4x64_ctr_t counter;
  9:   threefry4x64_key_t key;
 10:   threefry4x64_ctr_t result;
 11:   PetscInt           count;
 12: } PetscRandom123;

 14: R123_ULONG_LONG PETSCR123_SEED_0 = R123_64BIT(0x615D333D2655FE14);
 15: R123_ULONG_LONG PETSCR123_SEED_1 = R123_64BIT(0xAFF6369B3EE9FE96);
 16: R123_ULONG_LONG PETSCR123_SEED_2 = R123_64BIT(0x5956EBC717B60E07);
 17: R123_ULONG_LONG PETSCR123_SEED_3 = R123_64BIT(0xEE8612A0CBEABFF1);

 19: PetscErrorCode PetscRandomSeed_Random123(PetscRandom r)
 20: {
 21:   threefry4x64_ukey_t ukey;
 22:   PetscRandom123     *r123 = (PetscRandom123 *)r->data;

 24:   ukey.v[0] = (R123_ULONG_LONG)r->seed;
 25:   ukey.v[1] = PETSCR123_SEED_1;
 26:   ukey.v[2] = PETSCR123_SEED_2;
 27:   ukey.v[3] = PETSCR123_SEED_3;
 28:   /* The point of seeding should be that every time the sequence is seeded you get the same output.  In this CBRNG,
 29:    * that means we have to initialize the key and reset the counts */
 30:   r123->key          = threefry4x64keyinit(ukey);
 31:   r123->counter.v[0] = 0;
 32:   r123->counter.v[1] = 1;
 33:   r123->counter.v[2] = 2;
 34:   r123->counter.v[3] = 3;
 35:   r123->result       = threefry4x64(r123->counter, r123->key);
 36:   r123->count        = 0;
 37:   return 0;
 38: }

 40: static PetscReal PetscRandom123Step(PetscRandom123 *r123)
 41: {
 42:   PetscReal scale = ((PetscReal)1.) / (UINT64_MAX + ((PetscReal)1.));
 43:   PetscReal shift = .5 * scale;
 44:   PetscInt  mod   = (r123->count++) % 4;
 45:   PetscReal ret;

 47:   ret = r123->result.v[mod] * scale + shift;

 49:   if (mod == 3) {
 50:     r123->counter.v[0] += 4;
 51:     r123->counter.v[1] += 4;
 52:     r123->counter.v[2] += 4;
 53:     r123->counter.v[3] += 4;
 54:     r123->result = threefry4x64(r123->counter, r123->key);
 55:   }

 57:   return ret;
 58: }

 60: PetscErrorCode PetscRandomGetValue_Random123(PetscRandom r, PetscScalar *val)
 61: {
 62:   PetscRandom123 *r123 = (PetscRandom123 *)r->data;
 63:   PetscScalar     rscal;

 65: #if defined(PETSC_USE_COMPLEX)
 66:   {
 67:     PetscReal re = PetscRandom123Step(r123);
 68:     PetscReal im = PetscRandom123Step(r123);

 70:     if (r->iset) {
 71:       re = re * PetscRealPart(r->width) + PetscRealPart(r->low);
 72:       im = im * PetscImaginaryPart(r->width) + PetscImaginaryPart(r->low);
 73:     }

 75:     rscal = PetscCMPLX(re, im);
 76:   }
 77: #else
 78:   rscal = PetscRandom123Step(r123);
 79:   if (r->iset) rscal = rscal * r->width + r->low;
 80: #endif
 81:   *val = rscal;
 82:   return 0;
 83: }

 85: PetscErrorCode PetscRandomGetValueReal_Random123(PetscRandom r, PetscReal *val)
 86: {
 87:   PetscRandom123 *r123 = (PetscRandom123 *)r->data;
 88:   PetscReal       rreal;

 90:   rreal = PetscRandom123Step(r123);
 91:   if (r->iset) rreal = rreal * PetscRealPart(r->width) + PetscRealPart(r->low);
 92:   *val = rreal;
 93:   return 0;
 94: }

 96: PetscErrorCode PetscRandomDestroy_Random123(PetscRandom r)
 97: {
 98:   PetscFree(r->data);
 99:   return 0;
100: }

102: static struct _PetscRandomOps PetscRandomOps_Values = {
103:   PetscDesignatedInitializer(seed, PetscRandomSeed_Random123),
104:   PetscDesignatedInitializer(getvalue, PetscRandomGetValue_Random123),
105:   PetscDesignatedInitializer(getvaluereal, PetscRandomGetValueReal_Random123),
106:   PetscDesignatedInitializer(getvalues, NULL),
107:   PetscDesignatedInitializer(getvaluesreal, NULL),
108:   PetscDesignatedInitializer(destroy, PetscRandomDestroy_Random123),
109: };

111: /*MC
112:    PETSCRANDOM123 - access to Random123 counter based pseudorandom number generators (currently threefry4x64)

114:    Options Database Keys:
115: . -random_type <rand,rand48,sprng,random123> - select the random number generator at runtim

117:    Note:
118:    PETSc must be ./configure with the option --download-random123 to use this random number generator.

120:   Level: beginner

122: .seealso: `RandomCreate()`, `RandomSetType()`, `PETSCRAND`, `PETSCRAND48`, `PETSCSPRNG`, `PetscRandomSetFromOptions()`
123: M*/

125: PETSC_EXTERN PetscErrorCode PetscRandomCreate_Random123(PetscRandom r)
126: {
127:   PetscRandom123 *r123;

129:   PetscNew(&r123);
130:   r->data = r123;
131:   PetscMemcpy(r->ops, &PetscRandomOps_Values, sizeof(PetscRandomOps_Values));
132:   PetscObjectChangeTypeName((PetscObject)r, PETSCRANDOM123);
133:   PetscRandomSeed(r);
134:   return 0;
135: }