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