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: {
  9:   threefry4x64_ctr_t  counter;
 10:   threefry4x64_key_t  key;
 11:   threefry4x64_ctr_t  result;
 12:   PetscInt            count;
 13: }
 14: PetscRandom123;

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

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

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

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

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

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

 59:   return ret;
 60: }

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

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

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

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

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

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

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

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

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

116:    Options Database Keys:
117: . -random_type <rand,rand48,sprng,random123>

119:   Level: beginner

121:    PETSc must have been ./configure with the option --download-random123 to use
122:    this random number generator.

124: .seealso: RandomCreate(), RandomSetType(), PETSCRAND, PETSCRAND48, PETSCSPRNG
125: M*/

127: PETSC_EXTERN PetscErrorCode PetscRandomCreate_Random123(PetscRandom r)
128: {
129:   PetscRandom123 *r123;

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