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