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