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