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) PetscCall((*rnd->ops->setfromoptions)(PetscOptionsObject,rnd)); 174 PetscCall(PetscOptionsInt("-random_seed","Seed to use to generate random numbers","PetscRandomSetSeed",0,&seed,&set)); 175 if (set) { 176 PetscCall(PetscRandomSetSeed(rnd,(unsigned long int)seed)); 177 PetscCall(PetscRandomSeed(rnd)); 178 } 179 PetscCall(PetscOptionsBool("-random_no_imaginary_part","The imaginary part of the random number will be zero","PetscRandomSetInterval",noimaginary,&noimaginary,&set)); 180 #if defined(PETSC_HAVE_COMPLEX) 181 if (set) { 182 if (noimaginary) { 183 PetscScalar low,high; 184 PetscCall(PetscRandomGetInterval(rnd,&low,&high)); 185 low = low - PetscImaginaryPart(low); 186 high = high - PetscImaginaryPart(high); 187 PetscCall(PetscRandomSetInterval(rnd,low,high)); 188 } 189 } 190 #endif 191 PetscOptionsEnd(); 192 PetscCall(PetscRandomViewFromOptions(rnd,NULL, "-random_view")); 193 PetscFunctionReturn(0); 194 } 195 196 #if defined(PETSC_HAVE_SAWS) 197 #include <petscviewersaws.h> 198 #endif 199 200 /*@C 201 PetscRandomViewFromOptions - View from Options 202 203 Collective on PetscRandom 204 205 Input Parameters: 206 + A - the random number generator context 207 . obj - Optional object 208 - name - command line option 209 210 Level: intermediate 211 .seealso: `PetscRandom`, `PetscRandomView`, `PetscObjectViewFromOptions()`, `PetscRandomCreate()` 212 @*/ 213 PetscErrorCode PetscRandomViewFromOptions(PetscRandom A,PetscObject obj,const char name[]) 214 { 215 PetscFunctionBegin; 216 PetscValidHeaderSpecific(A,PETSC_RANDOM_CLASSID,1); 217 PetscCall(PetscObjectViewFromOptions((PetscObject)A,obj,name)); 218 PetscFunctionReturn(0); 219 } 220 221 /*@C 222 PetscRandomView - Views a random number generator object. 223 224 Collective on PetscRandom 225 226 Input Parameters: 227 + rnd - The random number generator context 228 - viewer - an optional visualization context 229 230 Notes: 231 The available visualization contexts include 232 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 233 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 234 output where only the first processor opens 235 the file. All other processors send their 236 data to the first processor to print. 237 238 You can change the format the vector is printed using the 239 option PetscViewerPushFormat(). 240 241 Level: beginner 242 243 .seealso: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()` 244 @*/ 245 PetscErrorCode PetscRandomView(PetscRandom rnd,PetscViewer viewer) 246 { 247 PetscBool iascii; 248 #if defined(PETSC_HAVE_SAWS) 249 PetscBool issaws; 250 #endif 251 252 PetscFunctionBegin; 253 PetscValidHeaderSpecific(rnd,PETSC_RANDOM_CLASSID,1); 254 PetscValidType(rnd,1); 255 if (!viewer) { 256 PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rnd),&viewer)); 257 } 258 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 259 PetscCheckSameComm(rnd,1,viewer,2); 260 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 261 #if defined(PETSC_HAVE_SAWS) 262 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws)); 263 #endif 264 if (iascii) { 265 PetscMPIInt rank; 266 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)rnd,viewer)); 267 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)rnd),&rank)); 268 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 269 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Random type %s, seed %lu\n",rank,((PetscObject)rnd)->type_name,rnd->seed)); 270 PetscCall(PetscViewerFlush(viewer)); 271 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 272 #if defined(PETSC_HAVE_SAWS) 273 } else if (issaws) { 274 PetscMPIInt rank; 275 const char *name; 276 277 PetscCall(PetscObjectGetName((PetscObject)rnd,&name)); 278 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 279 if (!((PetscObject)rnd)->amsmem && rank == 0) { 280 char dir[1024]; 281 282 PetscCall(PetscObjectViewSAWs((PetscObject)rnd,viewer)); 283 PetscCall(PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/Low",name)); 284 PetscStackCallSAWs(SAWs_Register,(dir,&rnd->low,1,SAWs_READ,SAWs_DOUBLE)); 285 } 286 #endif 287 } 288 PetscFunctionReturn(0); 289 } 290 291 /*@ 292 PetscRandomCreate - Creates a context for generating random numbers, 293 and initializes the random-number generator. 294 295 Collective 296 297 Input Parameters: 298 . comm - MPI communicator 299 300 Output Parameter: 301 . r - the random number generator context 302 303 Level: intermediate 304 305 Notes: 306 The random type has to be set by PetscRandomSetType(). 307 308 This is only a primitive "parallel" random number generator, it should NOT 309 be used for sophisticated parallel Monte Carlo methods since it will very likely 310 not have the correct statistics across processors. You can provide your own 311 parallel generator using PetscRandomRegister(); 312 313 If you create a PetscRandom() using PETSC_COMM_SELF on several processors then 314 the SAME random numbers will be generated on all those processors. Use PETSC_COMM_WORLD 315 or the appropriate parallel communicator to eliminate this issue. 316 317 Use VecSetRandom() to set the elements of a vector to random numbers. 318 319 Example of Usage: 320 .vb 321 PetscRandomCreate(PETSC_COMM_SELF,&r); 322 PetscRandomSetType(r,PETSCRAND48); 323 PetscRandomGetValue(r,&value1); 324 PetscRandomGetValueReal(r,&value2); 325 PetscRandomDestroy(&r); 326 .ve 327 328 .seealso: `PetscRandomSetType()`, `PetscRandomGetValue()`, `PetscRandomGetValueReal()`, `PetscRandomSetInterval()`, 329 `PetscRandomDestroy()`, `VecSetRandom()`, `PetscRandomType` 330 @*/ 331 332 PetscErrorCode PetscRandomCreate(MPI_Comm comm,PetscRandom *r) 333 { 334 PetscRandom rr; 335 PetscMPIInt rank; 336 337 PetscFunctionBegin; 338 PetscValidPointer(r,2); 339 *r = NULL; 340 PetscCall(PetscRandomInitializePackage()); 341 342 PetscCall(PetscHeaderCreate(rr,PETSC_RANDOM_CLASSID,"PetscRandom","Random number generator","Sys",comm,PetscRandomDestroy,PetscRandomView)); 343 344 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 345 346 rr->data = NULL; 347 rr->low = 0.0; 348 rr->width = 1.0; 349 rr->iset = PETSC_FALSE; 350 rr->seed = 0x12345678 + 76543*rank; 351 PetscCall(PetscRandomSetType(rr,PETSCRANDER48)); 352 *r = rr; 353 PetscFunctionReturn(0); 354 } 355 356 /*@ 357 PetscRandomSeed - Seed the generator. 358 359 Not collective 360 361 Input Parameters: 362 . r - The random number generator context 363 364 Level: intermediate 365 366 Usage: 367 PetscRandomSetSeed(r,a positive integer); 368 PetscRandomSeed(r); PetscRandomGetValue() will now start with the new seed. 369 370 PetscRandomSeed(r) without a call to PetscRandomSetSeed() re-initializes 371 the seed. The random numbers generated will be the same as before. 372 373 .seealso: `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSetSeed()` 374 @*/ 375 PetscErrorCode PetscRandomSeed(PetscRandom r) 376 { 377 PetscFunctionBegin; 378 PetscValidHeaderSpecific(r,PETSC_RANDOM_CLASSID,1); 379 PetscValidType(r,1); 380 381 PetscCall((*r->ops->seed)(r)); 382 PetscCall(PetscObjectStateIncrease((PetscObject)r)); 383 PetscFunctionReturn(0); 384 } 385