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