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 Parameters: 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 PetscValidPointer(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 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 .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 PetscValidPointer(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 Parameters: 359 . r - The random number generator context 360 361 Level: intermediate 362 363 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