PETCs - How should one initialize PetscRandom?

The documentation of PETSc shows an example of usage:

  PetscRandomCreate(PETSC_COMM_SELF,&r);
PetscRandomSetType(r,PETSCRAND48);
PetscRandomGetValue(r,&value1);
PetscRandomGetValueReal(r,&value2);
PetscRandomDestroy(&r);


In the example random 1, we can find also the use of PetscRandomSeed.

However, it seems not to be enough for attaching a PetscRandom variable to an object.

The following code works well for any number of procs:

#undef __FUNCT__
#define __FUNCT__ "TSInjectCreate"
PetscErrorCode  TSInjectCreate(MPI_Comm comm,TSInject *ininject){
PetscErrorCode ierr;
size_t         t;
PetscInt       i;
PetscReal       r;
TSInject       inject;

PetscFunctionBegin;
PetscValidPointer(ininject,1);
*ininject = NULL;
ierr                       = PetscClassIdRegister("TSInject",&TSINJECT_CLASSID);CHKERRQ(ierr);
ierr = PetscHeaderCreate(inject,TSINJECT_CLASSID,"TSInject","Injection of SDC","TS",comm,TSInjectDestroy, TSInjectView);CHKERRQ(ierr);
ierr = PetscRandomCreate( comm,&inject->ran);CHKERRQ(ierr);
ierr = PetscRandomSetType(inject->ran, PETSCRAND48);CHKERRQ(ierr);
ierr = PetscRandomSeed(inject->ran);CHKERRQ(ierr);
ierr = PetscRandomGetValue(inject->ran, &r);CHKERRQ(ierr);
*ininject = inject;

PetscFunctionReturn(0);
}


But as soon as I use PetscRandomGetValue(inject->ran, &r); in another function, it works only when mpi is launched with one proc.

Why?

EDIT:

The error comes exactly (and weirdly) from the comparison:

PetscErrorCode  TSInjection(TSInject inject, Vec X)
{
PetscReal       r,err;
Vec             Y;
PetscErrorCode ierr;

PetscFunctionBegin;
ierr = PetscRandomGetValue(inject->ran, &r);CHKERRQ(ierr);
ierr = VecView(X, 0);CHKERRQ(ierr);    // works
if(r  > inject->proba){ PetscFunctionReturn(0);}
ierr = VecView(X, 0);CHKERRQ(ierr);      // stucks the execution
...
PetscFunctionReturn(0);
}


If I replace r by 0.1 for example, everything works fine. What is all the more surprising is that displaying r doesn't stuck the execution.

• Is there a reason your inject object should survive passing out of scope? – origimbo Mar 22 '16 at 23:17
• Yes the reason is that it receives parameters and these parameters will choose which function should be called. In this initialization function, non important lines were deleted. – guhur Mar 23 '16 at 13:43
• Moreover, the problem occurs on several machines. Linux/mac, same issue. – guhur Mar 23 '16 at 15:32

To solve this issue, a work-around is too employ a DIY random function. In my context, I have given up with PetscRandom to write:

PetscErrorCode TSInjectRandInt(TSInject inject, PetscInt *r)  // RAND_MAX assumed to be 32767
{
PetscFunctionBegin;
inject->next = inject->next * 1103515245 + 12345;
*r = PetscAbs((PetscInt)(inject->next/65536) % 32768);
PetscFunctionReturn(0);
}

PetscErrorCode TSInjectRandReal(TSInject inject, PetscReal *r)  {
PetscInt i;

PetscFunctionBegin;
TSInjectRandInt(inject, &i);
*r = ((PetscReal)i)/32767.0;
PetscFunctionReturn(0);
}

PetscErrorCode TSInjectRandSeed(TSInject inject, PetscInt seed)
{
PetscFunctionBegin;
inject->next = seed;
PetscFunctionReturn(0);
}


However, one should notice that RAND_MAX is fixed and that the generator is clearly very easy. It has the only advantage of working.