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: static PetscErrorCode PetscRandomSeed_Random123(PetscRandom r)
     20: {
     21:   threefry4x64_ukey_t ukey;
     22:   PetscRandom123     *r123 = (PetscRandom123 *)r->data;

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

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

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

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

     58:   return ret;
     59: }

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

     66:   PetscFunctionBegin;
     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:   PetscFunctionReturn(PETSC_SUCCESS);
     85: }

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

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

     99: static PetscErrorCode PetscRandomGetValuesReal_Random123(PetscRandom r, PetscInt n, PetscReal vals[])
    100: {
    101:   PetscRandom123 *r123 = (PetscRandom123 *)r->data;
    102:   PetscInt        peel_start;
    103:   PetscInt        rem, lim;
    104:   PetscReal       scale = 1. / ((PetscReal)UINT64_MAX + 1.);
    105:   PetscReal       shift = .5 * scale;
    106:   PetscRandom123  r123_copy;

    108:   PetscFunctionBegin;
    109:   peel_start = (4 - (r123->count % 4)) % 4;
    110:   peel_start = PetscMin(n, peel_start);
    111:   for (PetscInt i = 0; i < peel_start; i++) PetscCall(PetscRandomGetValueReal(r, &vals[i]));
    112:   PetscAssert((r123->count % 4) == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad modular arithmetic");
    113:   n -= peel_start;
    114:   vals += peel_start;
    115:   rem = (n % 4);
    116:   lim = n - rem;
    117:   if (r->iset) {
    118:     scale *= PetscRealPart(r->width);
    119:     shift *= PetscRealPart(r->width);
    120:     shift += PetscRealPart(r->low);
    121:   }
    122:   r123_copy = *r123;
    123:   for (PetscInt i = 0; i < lim; i += 4, vals += 4) {
    124:     vals[0] = r123_copy.result.v[0] * scale + shift;
    125:     vals[1] = r123_copy.result.v[1] * scale + shift;
    126:     vals[2] = r123_copy.result.v[2] * scale + shift;
    127:     vals[3] = r123_copy.result.v[3] * scale + shift;
    128:     r123_copy.counter.v[0] += 4;
    129:     r123_copy.counter.v[1] += 4;
    130:     r123_copy.counter.v[2] += 4;
    131:     r123_copy.counter.v[3] += 4;
    132:     r123_copy.result = threefry4x64(r123->counter, r123->key);
    133:   }
    134:   r123_copy.count += lim;
    135:   *r123 = r123_copy;
    136:   for (PetscInt i = 0; i < rem; i++) PetscCall(PetscRandomGetValueReal(r, &vals[i]));
    137:   PetscFunctionReturn(PETSC_SUCCESS);
    138: }

    140: static PetscErrorCode PetscRandomGetValues_Random123(PetscRandom r, PetscInt n, PetscScalar vals[])
    141: {
    142:   PetscFunctionBegin;
    143: #if PetscDefined(USE_COMPLEX)
    144:   for (PetscInt i = 0; i < n; i++) PetscCall(PetscRandomGetValue_Random123(r, n, &vals[i]));
    145: #else
    146:   PetscCall(PetscRandomGetValuesReal_Random123(r, n, vals));
    147: #endif
    148:   PetscFunctionReturn(PETSC_SUCCESS);
    149: }

    151: static PetscErrorCode PetscRandomDestroy_Random123(PetscRandom r)
    152: {
    153:   PetscFunctionBegin;
    154:   PetscCall(PetscFree(r->data));
    155:   PetscFunctionReturn(PETSC_SUCCESS);
    156: }

    158: static struct _PetscRandomOps PetscRandomOps_Values = {
    159:   // clang-format off
    160:   PetscDesignatedInitializer(seed, PetscRandomSeed_Random123),
    161:   PetscDesignatedInitializer(getvalue, PetscRandomGetValue_Random123),
    162:   PetscDesignatedInitializer(getvaluereal, PetscRandomGetValueReal_Random123),
    163:   PetscDesignatedInitializer(getvalues, PetscRandomGetValues_Random123),
    164:   PetscDesignatedInitializer(getvaluesreal, PetscRandomGetValuesReal_Random123),
    165:   PetscDesignatedInitializer(destroy, PetscRandomDestroy_Random123),
    166:   // clang-format on
    167: };

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

    172:    Options Database Key:
    173: . -random_type (rander48|rand|rand48|sprng|random123|curand) - set the random number generator from the options database, see `PetscRandomType`.

    175:   Level: beginner

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

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

    183: PETSC_EXTERN PetscErrorCode PetscRandomCreate_Random123(PetscRandom r)
    184: {
    185:   PetscRandom123 *r123;

    187:   PetscFunctionBegin;
    188:   PetscCall(PetscNew(&r123));
    189:   r->data   = r123;
    190:   r->ops[0] = PetscRandomOps_Values;
    191:   PetscCall(PetscObjectChangeTypeName((PetscObject)r, PETSCRANDOM123));
    192:   PetscCall(PetscRandomSeed(r));
    193:   PetscFunctionReturn(PETSC_SUCCESS);
    194: }