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