xref: /petsc/src/sys/classes/random/interface/randomc.c (revision dae5874858be690a1464f3073a9d679bfd0ed673)
1 
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:
6 
7     Numbers can be shuffled for use as random tuples
8 
9     Multiple random number generators may be used
10 
11     We are still not sure what interface we want here.  There should be
12     one to reinitialize and set the seed.
13  */
14 
15 #include <../src/sys/classes/random/randomimpl.h>                              /*I "petscsys.h" I*/
16 #include <petscviewer.h>
17 
18 /* Logging support */
19 PetscClassId PETSC_RANDOM_CLASSID;
20 
21 #undef __FUNCT__
22 #define __FUNCT__ "PetscRandomDestroy"
23 /*@
24    PetscRandomDestroy - Destroys a context that has been formed by
25    PetscRandomCreate().
26 
27    Collective on PetscRandom
28 
29    Intput Parameter:
30 .  r  - the random number generator context
31 
32    Level: intermediate
33 
34 .seealso: PetscRandomGetValue(), PetscRandomCreate(), VecSetRandom()
35 @*/
36 PetscErrorCode  PetscRandomDestroy(PetscRandom *r)
37 {
38   PetscErrorCode ierr;
39 
40   PetscFunctionBegin;
41   if (!*r) PetscFunctionReturn(0);
42   PetscValidHeaderSpecific(*r,PETSC_RANDOM_CLASSID,1);
43   if (--((PetscObject)(*r))->refct > 0) {*r = 0; PetscFunctionReturn(0);}
44   ierr = PetscHeaderDestroy(r);CHKERRQ(ierr);
45   PetscFunctionReturn(0);
46 }
47 
48 
49 #undef __FUNCT__
50 #define __FUNCT__ "PetscRandomGetSeed"
51 /*@
52    PetscRandomGetSeed - Gets the random seed.
53 
54    Not collective
55 
56    Input Parameters:
57 .  r - The random number generator context
58 
59    Output Parameter:
60 .  seed - The random seed
61 
62    Level: intermediate
63 
64    Concepts: random numbers^seed
65 
66 .seealso: PetscRandomCreate(), PetscRandomSetSeed(), PetscRandomSeed()
67 @*/
68 PetscErrorCode  PetscRandomGetSeed(PetscRandom r,unsigned long *seed)
69 {
70   PetscFunctionBegin;
71   PetscValidHeaderSpecific(r,PETSC_RANDOM_CLASSID,1);
72   if (seed) {
73     PetscValidPointer(seed,2);
74     *seed = r->seed;
75   }
76   PetscFunctionReturn(0);
77 }
78 
79 #undef __FUNCT__
80 #define __FUNCT__ "PetscRandomSetSeed"
81 /*@
82    PetscRandomSetSeed - Sets the random seed. You MUST call PetscRandomSeed() after this call to have the new seed used.
83 
84    Not collective
85 
86    Input Parameters:
87 +  r  - The random number generator context
88 -  seed - The random seed
89 
90    Level: intermediate
91 
92    Usage:
93       PetscRandomSetSeed(r,a positive integer);
94       PetscRandomSeed(r);  PetscRandomGetValue() will now start with the new seed.
95 
96       PetscRandomSeed(r) without a call to PetscRandomSetSeed() re-initializes
97         the seed. The random numbers generated will be the same as before.
98 
99    Concepts: random numbers^seed
100 
101 .seealso: PetscRandomCreate(), PetscRandomGetSeed(), PetscRandomSeed()
102 @*/
103 PetscErrorCode  PetscRandomSetSeed(PetscRandom r,unsigned long seed)
104 {
105   PetscErrorCode ierr;
106 
107   PetscFunctionBegin;
108   PetscValidHeaderSpecific(r,PETSC_RANDOM_CLASSID,1);
109   r->seed = seed;
110   ierr    = PetscInfo1(NULL,"Setting seed to %d\n",(int)seed);CHKERRQ(ierr);
111   PetscFunctionReturn(0);
112 }
113 
114 /* ------------------------------------------------------------------- */
115 #undef __FUNCT__
116 #define __FUNCT__ "PetscRandomSetTypeFromOptions_Private"
117 /*
118   PetscRandomSetTypeFromOptions_Private - Sets the type of random generator from user options. Defaults to type PETSCRAND48 or PETSCRAND.
119 
120   Collective on PetscRandom
121 
122   Input Parameter:
123 . rnd - The random number generator context
124 
125   Level: intermediate
126 
127 .keywords: PetscRandom, set, options, database, type
128 .seealso: PetscRandomSetFromOptions(), PetscRandomSetType()
129 */
130 static PetscErrorCode PetscRandomSetTypeFromOptions_Private(PetscRandom rnd)
131 {
132   PetscBool      opt;
133   const char     *defaultType;
134   char           typeName[256];
135   PetscErrorCode ierr;
136 
137   PetscFunctionBegin;
138   if (((PetscObject)rnd)->type_name) {
139     defaultType = ((PetscObject)rnd)->type_name;
140   } else {
141 #if defined(PETSC_HAVE_DRAND48)
142     defaultType = PETSCRAND48;
143 #elif defined(PETSC_HAVE_RAND)
144     defaultType = PETSCRAND;
145 #endif
146   }
147 
148   if (!PetscRandomRegisterAllCalled) {ierr = PetscRandomRegisterAll();CHKERRQ(ierr);}
149   ierr = PetscOptionsList("-random_type","PetscRandom type","PetscRandomSetType",PetscRandomList,defaultType,typeName,256,&opt);CHKERRQ(ierr);
150   if (opt) {
151     ierr = PetscRandomSetType(rnd, typeName);CHKERRQ(ierr);
152   } else {
153     ierr = PetscRandomSetType(rnd, defaultType);CHKERRQ(ierr);
154   }
155   PetscFunctionReturn(0);
156 }
157 
158 #undef __FUNCT__
159 #define __FUNCT__ "PetscRandomSetFromOptions"
160 /*@
161   PetscRandomSetFromOptions - Configures the random number generator from the options database.
162 
163   Collective on PetscRandom
164 
165   Input Parameter:
166 . rnd - The random number generator context
167 
168   Options Database:
169 .  -random_seed <integer> - provide a seed to the random number generater
170 
171   Notes:  To see all options, run your program with the -help option.
172           Must be called after PetscRandomCreate() but before the rnd is used.
173 
174   Level: beginner
175 
176 .keywords: PetscRandom, set, options, database
177 .seealso: PetscRandomCreate(), PetscRandomSetType()
178 @*/
179 PetscErrorCode  PetscRandomSetFromOptions(PetscRandom rnd)
180 {
181   PetscErrorCode ierr;
182   PetscBool      set;
183   PetscInt       seed;
184 
185   PetscFunctionBegin;
186   PetscValidHeaderSpecific(rnd,PETSC_RANDOM_CLASSID,1);
187 
188   ierr = PetscObjectOptionsBegin((PetscObject)rnd);CHKERRQ(ierr);
189 
190   /* Handle PetscRandom type options */
191   ierr = PetscRandomSetTypeFromOptions_Private(rnd);CHKERRQ(ierr);
192 
193   /* Handle specific random generator's options */
194   if (rnd->ops->setfromoptions) {
195     ierr = (*rnd->ops->setfromoptions)(rnd);CHKERRQ(ierr);
196   }
197   ierr = PetscOptionsInt("-random_seed","Seed to use to generate random numbers","PetscRandomSetSeed",0,&seed,&set);CHKERRQ(ierr);
198   if (set) {
199     ierr = PetscRandomSetSeed(rnd,(unsigned long int)seed);CHKERRQ(ierr);
200     ierr = PetscRandomSeed(rnd);CHKERRQ(ierr);
201   }
202   ierr = PetscOptionsEnd();CHKERRQ(ierr);
203   ierr = PetscRandomViewFromOptions(rnd,NULL, "-random_view");CHKERRQ(ierr);
204   PetscFunctionReturn(0);
205 }
206 
207 #if defined(PETSC_HAVE_AMS)
208 #include <petscviewerams.h>
209 #endif
210 #undef __FUNCT__
211 #define __FUNCT__ "PetscRandomView"
212 /*@C
213    PetscRandomView - Views a random number generator object.
214 
215    Collective on PetscRandom
216 
217    Input Parameters:
218 +  rnd - The random number generator context
219 -  viewer - an optional visualization context
220 
221    Notes:
222    The available visualization contexts include
223 +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
224 -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
225          output where only the first processor opens
226          the file.  All other processors send their
227          data to the first processor to print.
228 
229    You can change the format the vector is printed using the
230    option PetscViewerSetFormat().
231 
232    Level: beginner
233 
234 .seealso:  PetscRealView(), PetscScalarView(), PetscIntView()
235 @*/
236 PetscErrorCode  PetscRandomView(PetscRandom rnd,PetscViewer viewer)
237 {
238   PetscErrorCode ierr;
239   PetscBool      iascii;
240 #if defined(PETSC_HAVE_AMS)
241   PetscBool      isams;
242 #endif
243 
244   PetscFunctionBegin;
245   PetscValidHeaderSpecific(rnd,PETSC_RANDOM_CLASSID,1);
246   PetscValidType(rnd,1);
247   if (!viewer) {
248     ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rnd),&viewer);CHKERRQ(ierr);
249   }
250   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
251   PetscCheckSameComm(rnd,1,viewer,2);
252   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
253 #if defined(PETSC_HAVE_AMS)
254   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERAMS,&isams);CHKERRQ(ierr);
255 #endif
256   if (iascii) {
257     PetscMPIInt rank;
258     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)rnd),&rank);CHKERRQ(ierr);
259     ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr);
260     ierr = PetscViewerASCIISynchronizedPrintf(viewer,"[%D] Random type %s, seed %D\n",rank,((PetscObject)rnd)->type_name,rnd->seed);CHKERRQ(ierr);
261     ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
262     ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr);
263 #if defined(PETSC_HAVE_AMS)
264   } else if (isams) {
265     if (((PetscObject)rnd)->amsmem == -1) {
266       ierr = PetscObjectViewAMS((PetscObject)rnd,viewer);CHKERRQ(ierr);
267       PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)rnd)->amsmem,"Low",&rnd->low,1,AMS_DOUBLE,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF));
268     }
269 #endif
270   }
271   PetscFunctionReturn(0);
272 }
273 
274 #undef  __FUNCT__
275 #define __FUNCT__ "PetscRandomViewFromOptions"
276 /*
277   PetscRandomViewFromOptions - This function visualizes the type and the seed of the generated random numbers based upon user options.
278 
279   Collective on PetscRandom
280 
281   Input Parameters:
282 + rnd   - The random number generator context
283 . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
284 - optionname - option to activate viewing
285 
286   Level: intermediate
287 
288 .keywords: PetscRandom, view, options, database
289 .seealso: PetscRandomSetFromOptions()
290 */
291 PetscErrorCode  PetscRandomViewFromOptions(PetscRandom rnd, const char prefix[], const char optionname[])
292 {
293   PetscBool         flg;
294   PetscViewer       viewer;
295   PetscErrorCode    ierr;
296   PetscViewerFormat format;
297 
298   PetscFunctionBegin;
299   if (prefix) {
300     ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)rnd),prefix,optionname,&viewer,&format,&flg);CHKERRQ(ierr);
301   } else {
302     ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)rnd),((PetscObject)rnd)->prefix,optionname,&viewer,&format,&flg);CHKERRQ(ierr);
303   }
304   if (flg) {
305     ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr);
306     ierr = PetscRandomView(rnd,viewer);CHKERRQ(ierr);
307     ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
308     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
309   }
310   PetscFunctionReturn(0);
311 }
312 
313 #undef __FUNCT__
314 #define __FUNCT__ "PetscRandomCreate"
315 /*@
316    PetscRandomCreate - Creates a context for generating random numbers,
317    and initializes the random-number generator.
318 
319    Collective on MPI_Comm
320 
321    Input Parameters:
322 +  comm - MPI communicator
323 
324    Output Parameter:
325 .  r  - the random number generator context
326 
327    Level: intermediate
328 
329    Notes:
330    The random type has to be set by PetscRandomSetType().
331 
332    This is only a primative "parallel" random number generator, it should NOT
333    be used for sophisticated parallel Monte Carlo methods since it will very likely
334    not have the correct statistics across processors. You can provide your own
335    parallel generator using PetscRandomRegister();
336 
337    If you create a PetscRandom() using PETSC_COMM_SELF on several processors then
338    the SAME random numbers will be generated on all those processors. Use PETSC_COMM_WORLD
339    or the appropriate parallel communicator to eliminate this issue.
340 
341    Use VecSetRandom() to set the elements of a vector to random numbers.
342 
343    Example of Usage:
344 .vb
345       PetscRandomCreate(PETSC_COMM_SELF,&r);
346       PetscRandomSetType(r,PETSCRAND48);
347       PetscRandomGetValue(r,&value1);
348       PetscRandomGetValueReal(r,&value2);
349       PetscRandomDestroy(&r);
350 .ve
351 
352    Concepts: random numbers^creating
353 
354 .seealso: PetscRandomSetType(), PetscRandomGetValue(), PetscRandomGetValueReal(), PetscRandomSetInterval(),
355           PetscRandomDestroy(), VecSetRandom(), PetscRandomType
356 @*/
357 
358 PetscErrorCode  PetscRandomCreate(MPI_Comm comm,PetscRandom *r)
359 {
360   PetscRandom    rr;
361   PetscErrorCode ierr;
362   PetscMPIInt    rank;
363 
364   PetscFunctionBegin;
365   PetscValidPointer(r,3);
366   *r = NULL;
367 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
368   ierr = PetscRandomInitializePackage();CHKERRQ(ierr);
369 #endif
370 
371   ierr = PetscHeaderCreate(rr,_p_PetscRandom,struct _PetscRandomOps,PETSC_RANDOM_CLASSID,"PetscRandom","Random number generator","Sys",comm,PetscRandomDestroy,0);CHKERRQ(ierr);
372 
373   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
374 
375   rr->data  = NULL;
376   rr->low   = 0.0;
377   rr->width = 1.0;
378   rr->iset  = PETSC_FALSE;
379   rr->seed  = 0x12345678 + 76543*rank;
380   *r = rr;
381   PetscFunctionReturn(0);
382 }
383 
384 #undef __FUNCT__
385 #define __FUNCT__ "PetscRandomSeed"
386 /*@
387    PetscRandomSeed - Seed the generator.
388 
389    Not collective
390 
391    Input Parameters:
392 .  r - The random number generator context
393 
394    Level: intermediate
395 
396    Usage:
397       PetscRandomSetSeed(r,a positive integer);
398       PetscRandomSeed(r);  PetscRandomGetValue() will now start with the new seed.
399 
400       PetscRandomSeed(r) without a call to PetscRandomSetSeed() re-initializes
401         the seed. The random numbers generated will be the same as before.
402 
403    Concepts: random numbers^seed
404 
405 .seealso: PetscRandomCreate(), PetscRandomGetSeed(), PetscRandomSetSeed()
406 @*/
407 PetscErrorCode  PetscRandomSeed(PetscRandom r)
408 {
409   PetscErrorCode ierr;
410 
411   PetscFunctionBegin;
412   PetscValidHeaderSpecific(r,PETSC_RANDOM_CLASSID,1);
413   PetscValidType(r,1);
414 
415   ierr = (*r->ops->seed)(r);CHKERRQ(ierr);
416   ierr = PetscObjectStateIncrease((PetscObject)r);CHKERRQ(ierr);
417   PetscFunctionReturn(0);
418 }
419 
420