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 PetscFunctionBegin; 36 if (!*r) PetscFunctionReturn(0); 37 PetscValidHeaderSpecific(*r, PETSC_RANDOM_CLASSID, 1); 38 if (--((PetscObject)(*r))->refct > 0) { 39 *r = NULL; 40 PetscFunctionReturn(0); 41 } 42 if ((*r)->ops->destroy) { PetscCall((*(*r)->ops->destroy)(*r)); } 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 PetscFunctionBegin; 64 PetscValidHeaderSpecific(r, PETSC_RANDOM_CLASSID, 1); 65 if (seed) { 66 PetscValidPointer(seed, 2); 67 *seed = r->seed; 68 } 69 PetscFunctionReturn(0); 70 } 71 72 /*@C 73 PetscRandomSetSeed - Sets the random seed. You MUST call PetscRandomSeed() after this call to have the new seed used. 74 75 Not collective 76 77 Input Parameters: 78 + r - The random number generator context 79 - seed - The random seed 80 81 Level: intermediate 82 83 Usage: 84 PetscRandomSetSeed(r,a positive integer); 85 PetscRandomSeed(r); PetscRandomGetValue() will now start with the new seed. 86 87 PetscRandomSeed(r) without a call to PetscRandomSetSeed() re-initializes 88 the seed. The random numbers generated will be the same as before. 89 90 .seealso: `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSeed()` 91 @*/ 92 PetscErrorCode PetscRandomSetSeed(PetscRandom r, unsigned long seed) { 93 PetscFunctionBegin; 94 PetscValidHeaderSpecific(r, PETSC_RANDOM_CLASSID, 1); 95 r->seed = seed; 96 PetscCall(PetscInfo(NULL, "Setting seed to %d\n", (int)seed)); 97 PetscFunctionReturn(0); 98 } 99 100 /* ------------------------------------------------------------------- */ 101 /* 102 PetscRandomSetTypeFromOptions_Private - Sets the type of random generator from user options. Defaults to type PETSCRAND48 or PETSCRAND. 103 104 Collective on PetscRandom 105 106 Input Parameter: 107 . rnd - The random number generator context 108 109 Level: intermediate 110 111 .seealso: `PetscRandomSetFromOptions()`, `PetscRandomSetType()` 112 */ 113 static PetscErrorCode PetscRandomSetTypeFromOptions_Private(PetscRandom rnd, PetscOptionItems *PetscOptionsObject) { 114 PetscBool opt; 115 const char *defaultType; 116 char typeName[256]; 117 118 PetscFunctionBegin; 119 if (((PetscObject)rnd)->type_name) { 120 defaultType = ((PetscObject)rnd)->type_name; 121 } else { 122 defaultType = PETSCRANDER48; 123 } 124 125 PetscCall(PetscRandomRegisterAll()); 126 PetscCall(PetscOptionsFList("-random_type", "PetscRandom type", "PetscRandomSetType", PetscRandomList, defaultType, typeName, 256, &opt)); 127 if (opt) { 128 PetscCall(PetscRandomSetType(rnd, typeName)); 129 } else { 130 PetscCall(PetscRandomSetType(rnd, defaultType)); 131 } 132 PetscFunctionReturn(0); 133 } 134 135 /*@ 136 PetscRandomSetFromOptions - Configures the random number generator from the options database. 137 138 Collective on PetscRandom 139 140 Input Parameter: 141 . rnd - The random number generator context 142 143 Options Database: 144 + -random_seed <integer> - provide a seed to the random number generater 145 - -random_no_imaginary_part - makes the imaginary part of the random number zero, this is useful when you want the 146 same code to produce the same result when run with real numbers or complex numbers for regression testing purposes 147 148 Notes: 149 To see all options, run your program with the -help option. 150 Must be called after PetscRandomCreate() but before the rnd is used. 151 152 Level: beginner 153 154 .seealso: `PetscRandomCreate()`, `PetscRandomSetType()` 155 @*/ 156 PetscErrorCode PetscRandomSetFromOptions(PetscRandom rnd) { 157 PetscBool set, noimaginary = PETSC_FALSE; 158 PetscInt seed; 159 160 PetscFunctionBegin; 161 PetscValidHeaderSpecific(rnd, PETSC_RANDOM_CLASSID, 1); 162 163 PetscObjectOptionsBegin((PetscObject)rnd); 164 165 /* Handle PetscRandom type options */ 166 PetscCall(PetscRandomSetTypeFromOptions_Private(rnd, PetscOptionsObject)); 167 168 /* Handle specific random generator's options */ 169 PetscTryTypeMethod(rnd, setfromoptions, PetscOptionsObject); 170 PetscCall(PetscOptionsInt("-random_seed", "Seed to use to generate random numbers", "PetscRandomSetSeed", 0, &seed, &set)); 171 if (set) { 172 PetscCall(PetscRandomSetSeed(rnd, (unsigned long int)seed)); 173 PetscCall(PetscRandomSeed(rnd)); 174 } 175 PetscCall(PetscOptionsBool("-random_no_imaginary_part", "The imaginary part of the random number will be zero", "PetscRandomSetInterval", noimaginary, &noimaginary, &set)); 176 #if defined(PETSC_HAVE_COMPLEX) 177 if (set) { 178 if (noimaginary) { 179 PetscScalar low, high; 180 PetscCall(PetscRandomGetInterval(rnd, &low, &high)); 181 low = low - PetscImaginaryPart(low); 182 high = high - PetscImaginaryPart(high); 183 PetscCall(PetscRandomSetInterval(rnd, low, high)); 184 } 185 } 186 #endif 187 PetscOptionsEnd(); 188 PetscCall(PetscRandomViewFromOptions(rnd, NULL, "-random_view")); 189 PetscFunctionReturn(0); 190 } 191 192 #if defined(PETSC_HAVE_SAWS) 193 #include <petscviewersaws.h> 194 #endif 195 196 /*@C 197 PetscRandomViewFromOptions - View from Options 198 199 Collective on PetscRandom 200 201 Input Parameters: 202 + A - the random number generator context 203 . obj - Optional object 204 - name - command line option 205 206 Level: intermediate 207 .seealso: `PetscRandom`, `PetscRandomView`, `PetscObjectViewFromOptions()`, `PetscRandomCreate()` 208 @*/ 209 PetscErrorCode PetscRandomViewFromOptions(PetscRandom A, PetscObject obj, const char name[]) { 210 PetscFunctionBegin; 211 PetscValidHeaderSpecific(A, PETSC_RANDOM_CLASSID, 1); 212 PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name)); 213 PetscFunctionReturn(0); 214 } 215 216 /*@C 217 PetscRandomView - Views a random number generator object. 218 219 Collective on PetscRandom 220 221 Input Parameters: 222 + rnd - The random number generator context 223 - viewer - an optional visualization context 224 225 Notes: 226 The available visualization contexts include 227 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 228 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 229 output where only the first processor opens 230 the file. All other processors send their 231 data to the first processor to print. 232 233 You can change the format the vector is printed using the 234 option PetscViewerPushFormat(). 235 236 Level: beginner 237 238 .seealso: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()` 239 @*/ 240 PetscErrorCode PetscRandomView(PetscRandom rnd, PetscViewer viewer) { 241 PetscBool iascii; 242 #if defined(PETSC_HAVE_SAWS) 243 PetscBool issaws; 244 #endif 245 246 PetscFunctionBegin; 247 PetscValidHeaderSpecific(rnd, PETSC_RANDOM_CLASSID, 1); 248 PetscValidType(rnd, 1); 249 if (!viewer) { PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rnd), &viewer)); } 250 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 251 PetscCheckSameComm(rnd, 1, viewer, 2); 252 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 253 #if defined(PETSC_HAVE_SAWS) 254 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws)); 255 #endif 256 if (iascii) { 257 PetscMPIInt rank; 258 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)rnd, viewer)); 259 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)rnd), &rank)); 260 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 261 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Random type %s, seed %lu\n", rank, ((PetscObject)rnd)->type_name, rnd->seed)); 262 PetscCall(PetscViewerFlush(viewer)); 263 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 264 #if defined(PETSC_HAVE_SAWS) 265 } else if (issaws) { 266 PetscMPIInt rank; 267 const char *name; 268 269 PetscCall(PetscObjectGetName((PetscObject)rnd, &name)); 270 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 271 if (!((PetscObject)rnd)->amsmem && rank == 0) { 272 char dir[1024]; 273 274 PetscCall(PetscObjectViewSAWs((PetscObject)rnd, viewer)); 275 PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/Low", name)); 276 PetscCallSAWs(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 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 primitive "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 .seealso: `PetscRandomSetType()`, `PetscRandomGetValue()`, `PetscRandomGetValueReal()`, `PetscRandomSetInterval()`, 321 `PetscRandomDestroy()`, `VecSetRandom()`, `PetscRandomType` 322 @*/ 323 324 PetscErrorCode PetscRandomCreate(MPI_Comm comm, PetscRandom *r) { 325 PetscRandom rr; 326 PetscMPIInt rank; 327 328 PetscFunctionBegin; 329 PetscValidPointer(r, 2); 330 *r = NULL; 331 PetscCall(PetscRandomInitializePackage()); 332 333 PetscCall(PetscHeaderCreate(rr, PETSC_RANDOM_CLASSID, "PetscRandom", "Random number generator", "Sys", comm, PetscRandomDestroy, PetscRandomView)); 334 335 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 336 337 rr->data = NULL; 338 rr->low = 0.0; 339 rr->width = 1.0; 340 rr->iset = PETSC_FALSE; 341 rr->seed = 0x12345678 + 76543 * rank; 342 PetscCall(PetscRandomSetType(rr, PETSCRANDER48)); 343 *r = rr; 344 PetscFunctionReturn(0); 345 } 346 347 /*@ 348 PetscRandomSeed - Seed the generator. 349 350 Not collective 351 352 Input Parameters: 353 . r - The random number generator context 354 355 Level: intermediate 356 357 Usage: 358 PetscRandomSetSeed(r,a positive integer); 359 PetscRandomSeed(r); PetscRandomGetValue() will now start with the new seed. 360 361 PetscRandomSeed(r) without a call to PetscRandomSetSeed() re-initializes 362 the seed. The random numbers generated will be the same as before. 363 364 .seealso: `PetscRandomCreate()`, `PetscRandomGetSeed()`, `PetscRandomSetSeed()` 365 @*/ 366 PetscErrorCode PetscRandomSeed(PetscRandom r) { 367 PetscFunctionBegin; 368 PetscValidHeaderSpecific(r, PETSC_RANDOM_CLASSID, 1); 369 PetscValidType(r, 1); 370 371 PetscUseTypeMethod(r, seed); 372 PetscCall(PetscObjectStateIncrease((PetscObject)r)); 373 PetscFunctionReturn(0); 374 } 375