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