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 <petsc/private/randomimpl.h> /*I "petscsys.h" I*/ 16 #include <petscviewer.h> 17 18 /* Logging support */ 19 PetscClassId PETSC_RANDOM_CLASSID; 20 21 /*@C 22 PetscRandomDestroy - Destroys a context that has been formed by 23 PetscRandomCreate(). 24 25 Collective on PetscRandom 26 27 Input Parameter: 28 . r - the random number generator context 29 30 Level: intermediate 31 32 .seealso: `PetscRandomGetValue()`, `PetscRandomCreate()`, `VecSetRandom()` 33 @*/ 34 PetscErrorCode PetscRandomDestroy(PetscRandom *r) 35 { 36 PetscFunctionBegin; 37 if (!*r) PetscFunctionReturn(0); 38 PetscValidHeaderSpecific(*r,PETSC_RANDOM_CLASSID,1); 39 if (--((PetscObject)(*r))->refct > 0) {*r = NULL; PetscFunctionReturn(0);} 40 if ((*r)->ops->destroy) { 41 PetscCall((*(*r)->ops->destroy)(*r)); 42 } 43 PetscCall(PetscHeaderDestroy(r)); 44 PetscFunctionReturn(0); 45 } 46 47 /*@C 48 PetscRandomGetSeed - Gets the random seed. 49 50 Not collective 51 52 Input Parameters: 53 . r - The random number generator context 54 55 Output Parameter: 56 . seed - The random seed 57 58 Level: intermediate 59 60 .seealso: `PetscRandomCreate()`, `PetscRandomSetSeed()`, `PetscRandomSeed()` 61 @*/ 62 PetscErrorCode PetscRandomGetSeed(PetscRandom r,unsigned long *seed) 63 { 64 PetscFunctionBegin; 65 PetscValidHeaderSpecific(r,PETSC_RANDOM_CLASSID,1); 66 if (seed) { 67 PetscValidPointer(seed,2); 68 *seed = r->seed; 69 } 70 PetscFunctionReturn(0); 71 } 72 73 /*@C 74 PetscRandomSetSeed - Sets the random seed. You MUST call PetscRandomSeed() after this call to have the new seed used. 75 76 Not collective 77 78 Input Parameters: 79 + r - The random number generator context 80 - seed - The random seed 81 82 Level: intermediate 83 84 Usage: 85 PetscRandomSetSeed(r,a positive integer); 86 PetscRandomSeed(r); PetscRandomGetValue() will now start with the new seed. 87 88 PetscRandomSeed(r) without a call to PetscRandomSetSeed() re-initializes 89 the seed. The random numbers generated will be the same as before. 90 91 .seealso: `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSeed()` 92 @*/ 93 PetscErrorCode PetscRandomSetSeed(PetscRandom r,unsigned long seed) 94 { 95 PetscFunctionBegin; 96 PetscValidHeaderSpecific(r,PETSC_RANDOM_CLASSID,1); 97 r->seed = seed; 98 PetscCall(PetscInfo(NULL,"Setting seed to %d\n",(int)seed)); 99 PetscFunctionReturn(0); 100 } 101 102 /* ------------------------------------------------------------------- */ 103 /* 104 PetscRandomSetTypeFromOptions_Private - Sets the type of random generator from user options. Defaults to type PETSCRAND48 or PETSCRAND. 105 106 Collective on PetscRandom 107 108 Input Parameter: 109 . rnd - The random number generator context 110 111 Level: intermediate 112 113 .seealso: `PetscRandomSetFromOptions()`, `PetscRandomSetType()` 114 */ 115 static PetscErrorCode PetscRandomSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,PetscRandom rnd) 116 { 117 PetscBool opt; 118 const char *defaultType; 119 char typeName[256]; 120 121 PetscFunctionBegin; 122 if (((PetscObject)rnd)->type_name) { 123 defaultType = ((PetscObject)rnd)->type_name; 124 } else { 125 defaultType = PETSCRANDER48; 126 } 127 128 PetscCall(PetscRandomRegisterAll()); 129 PetscCall(PetscOptionsFList("-random_type","PetscRandom type","PetscRandomSetType",PetscRandomList,defaultType,typeName,256,&opt)); 130 if (opt) { 131 PetscCall(PetscRandomSetType(rnd, typeName)); 132 } else { 133 PetscCall(PetscRandomSetType(rnd, defaultType)); 134 } 135 PetscFunctionReturn(0); 136 } 137 138 /*@ 139 PetscRandomSetFromOptions - Configures the random number generator from the options database. 140 141 Collective on PetscRandom 142 143 Input Parameter: 144 . rnd - The random number generator context 145 146 Options Database: 147 + -random_seed <integer> - provide a seed to the random number generater 148 - -random_no_imaginary_part - makes the imaginary part of the random number zero, this is useful when you want the 149 same code to produce the same result when run with real numbers or complex numbers for regression testing purposes 150 151 Notes: 152 To see all options, run your program with the -help option. 153 Must be called after PetscRandomCreate() but before the rnd is used. 154 155 Level: beginner 156 157 .seealso: `PetscRandomCreate()`, `PetscRandomSetType()` 158 @*/ 159 PetscErrorCode PetscRandomSetFromOptions(PetscRandom rnd) 160 { 161 PetscBool set,noimaginary = PETSC_FALSE; 162 PetscInt seed; 163 164 PetscFunctionBegin; 165 PetscValidHeaderSpecific(rnd,PETSC_RANDOM_CLASSID,1); 166 167 PetscObjectOptionsBegin((PetscObject)rnd); 168 169 /* Handle PetscRandom type options */ 170 PetscCall(PetscRandomSetTypeFromOptions_Private(PetscOptionsObject,rnd)); 171 172 /* Handle specific random generator's options */ 173 if (rnd->ops->setfromoptions) { 174 PetscCall((*rnd->ops->setfromoptions)(PetscOptionsObject,rnd)); 175 } 176 PetscCall(PetscOptionsInt("-random_seed","Seed to use to generate random numbers","PetscRandomSetSeed",0,&seed,&set)); 177 if (set) { 178 PetscCall(PetscRandomSetSeed(rnd,(unsigned long int)seed)); 179 PetscCall(PetscRandomSeed(rnd)); 180 } 181 PetscCall(PetscOptionsBool("-random_no_imaginary_part","The imaginary part of the random number will be zero","PetscRandomSetInterval",noimaginary,&noimaginary,&set)); 182 #if defined(PETSC_HAVE_COMPLEX) 183 if (set) { 184 if (noimaginary) { 185 PetscScalar low,high; 186 PetscCall(PetscRandomGetInterval(rnd,&low,&high)); 187 low = low - PetscImaginaryPart(low); 188 high = high - PetscImaginaryPart(high); 189 PetscCall(PetscRandomSetInterval(rnd,low,high)); 190 } 191 } 192 #endif 193 PetscOptionsEnd(); 194 PetscCall(PetscRandomViewFromOptions(rnd,NULL, "-random_view")); 195 PetscFunctionReturn(0); 196 } 197 198 #if defined(PETSC_HAVE_SAWS) 199 #include <petscviewersaws.h> 200 #endif 201 202 /*@C 203 PetscRandomViewFromOptions - View from Options 204 205 Collective on PetscRandom 206 207 Input Parameters: 208 + A - the random number generator context 209 . obj - Optional object 210 - name - command line option 211 212 Level: intermediate 213 .seealso: `PetscRandom`, `PetscRandomView`, `PetscObjectViewFromOptions()`, `PetscRandomCreate()` 214 @*/ 215 PetscErrorCode PetscRandomViewFromOptions(PetscRandom A,PetscObject obj,const char name[]) 216 { 217 PetscFunctionBegin; 218 PetscValidHeaderSpecific(A,PETSC_RANDOM_CLASSID,1); 219 PetscCall(PetscObjectViewFromOptions((PetscObject)A,obj,name)); 220 PetscFunctionReturn(0); 221 } 222 223 /*@C 224 PetscRandomView - Views a random number generator object. 225 226 Collective on PetscRandom 227 228 Input Parameters: 229 + rnd - The random number generator context 230 - viewer - an optional visualization context 231 232 Notes: 233 The available visualization contexts include 234 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 235 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 236 output where only the first processor opens 237 the file. All other processors send their 238 data to the first processor to print. 239 240 You can change the format the vector is printed using the 241 option PetscViewerPushFormat(). 242 243 Level: beginner 244 245 .seealso: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()` 246 @*/ 247 PetscErrorCode PetscRandomView(PetscRandom rnd,PetscViewer viewer) 248 { 249 PetscBool iascii; 250 #if defined(PETSC_HAVE_SAWS) 251 PetscBool issaws; 252 #endif 253 254 PetscFunctionBegin; 255 PetscValidHeaderSpecific(rnd,PETSC_RANDOM_CLASSID,1); 256 PetscValidType(rnd,1); 257 if (!viewer) { 258 PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rnd),&viewer)); 259 } 260 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 261 PetscCheckSameComm(rnd,1,viewer,2); 262 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 263 #if defined(PETSC_HAVE_SAWS) 264 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws)); 265 #endif 266 if (iascii) { 267 PetscMPIInt rank; 268 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)rnd,viewer)); 269 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)rnd),&rank)); 270 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 271 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Random type %s, seed %lu\n",rank,((PetscObject)rnd)->type_name,rnd->seed)); 272 PetscCall(PetscViewerFlush(viewer)); 273 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 274 #if defined(PETSC_HAVE_SAWS) 275 } else if (issaws) { 276 PetscMPIInt rank; 277 const char *name; 278 279 PetscCall(PetscObjectGetName((PetscObject)rnd,&name)); 280 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 281 if (!((PetscObject)rnd)->amsmem && rank == 0) { 282 char dir[1024]; 283 284 PetscCall(PetscObjectViewSAWs((PetscObject)rnd,viewer)); 285 PetscCall(PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/Low",name)); 286 PetscStackCallSAWs(SAWs_Register,(dir,&rnd->low,1,SAWs_READ,SAWs_DOUBLE)); 287 } 288 #endif 289 } 290 PetscFunctionReturn(0); 291 } 292 293 /*@ 294 PetscRandomCreate - Creates a context for generating random numbers, 295 and initializes the random-number generator. 296 297 Collective 298 299 Input Parameters: 300 . comm - MPI communicator 301 302 Output Parameter: 303 . r - the random number generator context 304 305 Level: intermediate 306 307 Notes: 308 The random type has to be set by PetscRandomSetType(). 309 310 This is only a primative "parallel" random number generator, it should NOT 311 be used for sophisticated parallel Monte Carlo methods since it will very likely 312 not have the correct statistics across processors. You can provide your own 313 parallel generator using PetscRandomRegister(); 314 315 If you create a PetscRandom() using PETSC_COMM_SELF on several processors then 316 the SAME random numbers will be generated on all those processors. Use PETSC_COMM_WORLD 317 or the appropriate parallel communicator to eliminate this issue. 318 319 Use VecSetRandom() to set the elements of a vector to random numbers. 320 321 Example of Usage: 322 .vb 323 PetscRandomCreate(PETSC_COMM_SELF,&r); 324 PetscRandomSetType(r,PETSCRAND48); 325 PetscRandomGetValue(r,&value1); 326 PetscRandomGetValueReal(r,&value2); 327 PetscRandomDestroy(&r); 328 .ve 329 330 .seealso: `PetscRandomSetType()`, `PetscRandomGetValue()`, `PetscRandomGetValueReal()`, `PetscRandomSetInterval()`, 331 `PetscRandomDestroy()`, `VecSetRandom()`, `PetscRandomType` 332 @*/ 333 334 PetscErrorCode PetscRandomCreate(MPI_Comm comm,PetscRandom *r) 335 { 336 PetscRandom rr; 337 PetscMPIInt rank; 338 339 PetscFunctionBegin; 340 PetscValidPointer(r,2); 341 *r = NULL; 342 PetscCall(PetscRandomInitializePackage()); 343 344 PetscCall(PetscHeaderCreate(rr,PETSC_RANDOM_CLASSID,"PetscRandom","Random number generator","Sys",comm,PetscRandomDestroy,PetscRandomView)); 345 346 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 347 348 rr->data = NULL; 349 rr->low = 0.0; 350 rr->width = 1.0; 351 rr->iset = PETSC_FALSE; 352 rr->seed = 0x12345678 + 76543*rank; 353 PetscCall(PetscRandomSetType(rr,PETSCRANDER48)); 354 *r = rr; 355 PetscFunctionReturn(0); 356 } 357 358 /*@ 359 PetscRandomSeed - Seed the generator. 360 361 Not collective 362 363 Input Parameters: 364 . r - The random number generator context 365 366 Level: intermediate 367 368 Usage: 369 PetscRandomSetSeed(r,a positive integer); 370 PetscRandomSeed(r); PetscRandomGetValue() will now start with the new seed. 371 372 PetscRandomSeed(r) without a call to PetscRandomSetSeed() re-initializes 373 the seed. The random numbers generated will be the same as before. 374 375 .seealso: `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSetSeed()` 376 @*/ 377 PetscErrorCode PetscRandomSeed(PetscRandom r) 378 { 379 PetscFunctionBegin; 380 PetscValidHeaderSpecific(r,PETSC_RANDOM_CLASSID,1); 381 PetscValidType(r,1); 382 383 PetscCall((*r->ops->seed)(r)); 384 PetscCall(PetscObjectStateIncrease((PetscObject)r)); 385 PetscFunctionReturn(0); 386 } 387