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