Actual source code: random.c
1: /*$Id: random.c,v 1.60 2001/08/07 03:02:06 balay Exp $*/
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:
7: Numbers can be shuffled for use as random tuples
9: Multiple random number generators may be used
11: We are still not sure what interface we want here. There should be
12: one to reinitialize and set the seed.
13: */
15: #include petsc.h
16: #include petscsys.h
17: #include <stdlib.h>
19: int PETSC_RANDOM_COOKIE;
21: /* Private data */
22: struct _p_PetscRandom {
23: PETSCHEADER(int)
24: unsigned long seed;
25: PetscScalar low,width; /* lower bound and width of the interval over
26: which the random numbers are distributed */
27: PetscTruth iset; /* if true, indicates that the user has set the interval */
28: /* array for shuffling ??? */
29: };
31: #undef __FUNCT__
33: /*@C
34: PetscRandomDestroy - Destroys a context that has been formed by
35: PetscRandomCreate().
37: Collective on PetscRandom
39: Intput Parameter:
40: . r - the random number generator context
42: Level: intermediate
44: .seealso: PetscRandomGetValue(), PetscRandomCreate(), VecSetRandom()
45: @*/
46: int PetscRandomDestroy(PetscRandom r)
47: {
50: if (--r->refct > 0) return(0);
52: PetscLogObjectDestroy((PetscObject)r);
53: PetscHeaderDestroy((PetscObject)r);
54: return(0);
55: }
57: #undef __FUNCT__
59: /*@C
60: PetscRandomSetInterval - Sets the interval over which the random numbers
61: will be randomly distributed. By default, this interval is [0,1).
63: Collective on PetscRandom
65: Input Parameters:
66: . r - the random number generator context
68: Example of Usage:
69: .vb
70: PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&r);
71: PetscRandomSetInterval(RANDOM_DEFAULT,&r);
72: PetscRandomGetValue(r,&value1);
73: PetscRandomGetValue(r,&value2);
74: PetscRandomDestroy(r);
75: .ve
77: Level: intermediate
79: Concepts: random numbers^range
81: .seealso: PetscRandomCreate()
82: @*/
83: int PetscRandomSetInterval(PetscRandom r,PetscScalar low,PetscScalar high)
84: {
87: #if defined(PETSC_USE_COMPLEX)
88: if (PetscRealPart(low) >= PetscRealPart(high)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"only low < high");
89: if (PetscImaginaryPart(low) >= PetscImaginaryPart(high)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"only low < high");
90: #else
91: if (low >= high) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"only low < high: Instead %g %g",low,high);
92: #endif
93: r->low = low;
94: r->width = high-low;
95: r->iset = PETSC_TRUE;
96: return(0);
97: }
99: /*
100: For now we have set up using the DRAND48() generater. We need to deal
101: with other variants of random number generators. We should also add
102: a routine to enable restarts [seed48()]
103: */
104: #if defined(PETSC_HAVE_DRAND48)
105: EXTERN_C_BEGIN
106: extern double drand48();
107: extern void srand48(long);
108: EXTERN_C_END
110: #undef __FUNCT__
112: /*@C
113: PetscRandomCreate - Creates a context for generating random numbers,
114: and initializes the random-number generator.
116: Collective on MPI_Comm
118: Input Parameters:
119: + comm - MPI communicator
120: - type - the type of random numbers to be generated, usually RANDOM_DEFAULT
122: Output Parameter:
123: . r - the random number generator context
125: Level: intermediate
127: Notes:
128: By default, we generate random numbers via srand48()/drand48() that
129: are uniformly distributed over [0,1). The user can shift and stretch
130: this interval by calling PetscRandomSetInterval().
131:
132: Currently three types of random numbers are supported. These types
133: are equivalent when working with real numbers.
134: . RANDOM_DEFAULT - both real and imaginary components are random
135: . RANDOM_DEFAULT_REAL - real component is random; imaginary component is 0
136: . RANDOM_DEFAULT_IMAGINARY - imaginary component is random; real component is 0
138: Use VecSetRandom() to set the elements of a vector to random numbers.
140: Example of Usage:
141: .vb
142: PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&r);
143: PetscRandomGetValue(r,&value1);
144: PetscRandomGetValue(r,&value2);
145: PetscRandomGetValue(r,&value3);
146: PetscRandomDestroy(r);
147: .ve
149: Concepts: random numbers^creating
151: .seealso: PetscRandomGetValue(), PetscRandomSetInterval(), PetscRandomDestroy(), VecSetRandom()
152: @*/
153: int PetscRandomCreate(MPI_Comm comm,PetscRandomType type,PetscRandom *r)
154: {
155: PetscRandom rr;
156: int ierr,rank;
159: *r = 0;
160: if (type != RANDOM_DEFAULT && type != RANDOM_DEFAULT_REAL && type != RANDOM_DEFAULT_IMAGINARY){
161: SETERRQ(PETSC_ERR_SUP,"Not for this random number type");
162: }
163: PetscHeaderCreate(rr,_p_PetscRandom,int,PETSC_RANDOM_COOKIE,type,"random",comm,PetscRandomDestroy,0);
164: PetscLogObjectCreate(rr);
165: rr->low = 0.0;
166: rr->width = 1.0;
167: rr->iset = PETSC_FALSE;
168: rr->seed = 0;
169: MPI_Comm_rank(comm,&rank);
170: srand48(0x12345678+rank);
171: *r = rr;
172: return(0);
173: }
175: #undef __FUNCT__
177: /*@C
178: PetscRandomGetValue - Generates a random number. Call this after first calling
179: PetscRandomCreate().
181: Not Collective
183: Intput Parameter:
184: . r - the random number generator context
186: Output Parameter:
187: . val - the value
189: Level: intermediate
191: Notes:
192: Use VecSetRandom() to set the elements of a vector to random numbers.
194: Example of Usage:
195: .vb
196: PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&r);
197: PetscRandomGetValue(r,&value1);
198: PetscRandomGetValue(r,&value2);
199: PetscRandomGetValue(r,&value3);
200: PetscRandomDestroy(r);
201: .ve
203: Concepts: random numbers^getting
205: .seealso: PetscRandomCreate(), PetscRandomDestroy(), VecSetRandom()
206: @*/
207: int PetscRandomGetValue(PetscRandom r,PetscScalar *val)
208: {
211: #if defined(PETSC_USE_COMPLEX)
212: if (r->type == RANDOM_DEFAULT) {
213: if (r->iset == PETSC_TRUE) {
214: *val = PetscRealPart(r->width)*drand48() + PetscRealPart(r->low) +
215: (PetscImaginaryPart(r->width)*drand48() + PetscImaginaryPart(r->low)) * PETSC_i;
216: }
217: else *val = drand48() + drand48()*PETSC_i;
218: } else if (r->type == RANDOM_DEFAULT_REAL) {
219: if (r->iset == PETSC_TRUE) *val = PetscRealPart(r->width)*drand48() + PetscRealPart(r->low);
220: else *val = drand48();
221: } else if (r->type == RANDOM_DEFAULT_IMAGINARY) {
222: if (r->iset == PETSC_TRUE) *val = (PetscImaginaryPart(r->width)*drand48()+PetscImaginaryPart(r->low))*PETSC_i;
223: else *val = drand48()*PETSC_i;
224: } else {
225: SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid random number type");
226: }
227: #else
228: if (r->iset == PETSC_TRUE) *val = r->width * drand48() + r->low;
229: else *val = drand48();
230: #endif
231: return(0);
232: }
234: #elif defined(PETSC_HAVE_RAND)
236: #undef __FUNCT__
238: int PetscRandomCreate(MPI_Comm comm,PetscRandomType type,PetscRandom *r)
239: {
240: PetscRandom rr;
241: int rank,ierr;
244: PetscLogInfo(0,"PetscRandomCreate: using rand(). not as efficinet as dran48n");
245: *r = 0;
246: if (type != RANDOM_DEFAULT && type != RANDOM_DEFAULT_REAL && type != RANDOM_DEFAULT_IMAGINARY) {
247: SETERRQ(PETSC_ERR_SUP,"Not for this random number type");
248: }
249: PetscHeaderCreate(rr,_p_PetscRandom,int,PETSC_RANDOM_COOKIE,type,"random",comm,PetscRandomDestroy,0);
250: PetscLogObjectCreate(rr);
251: rr->low = 0.0;
252: rr->width = 1.0;
253: rr->iset = PETSC_FALSE;
254: rr->seed = 0;
255: MPI_Comm_rank(comm,&rank);
256: srand(0x12345678+rank);
257: *r = rr;
258: return(0);
259: }
261: #define RAND_WRAP() (rand()/(double)((unsigned int)RAND_MAX+1))
262: #undef __FUNCT__
264: int PetscRandomGetValue(PetscRandom r,PetscScalar *val)
265: {
268: #if defined(PETSC_USE_COMPLEX)
269: if (r->type == RANDOM_DEFAULT) {
270: if (r->iset == PETSC_TRUE)
271: *val = PetscRealPart(r->width)*RAND_WRAP() + PetscRealPart(r->low) +
272: (PetscImaginaryPart(r->width)*RAND_WRAP() + PetscImaginaryPart(r->low)) * PETSC_i;
273: else *val = RAND_WRAP() + RAND_WRAP()*PETSC_i;
274: } else if (r->type == RANDOM_DEFAULT_REAL) {
275: if (r->iset == PETSC_TRUE) *val = PetscRealPart(r->width)*RAND_WRAP() + PetscRealPart(r->low);
276: else *val = RAND_WRAP();
277: } else if (r->type == RANDOM_DEFAULT_IMAGINARY) {
278: if (r->iset == PETSC_TRUE) *val = (PetscImaginaryPart(r->width)*RAND_WRAP()+PetscImaginaryPart(r->low))*PETSC_i;
279: else *val = RAND_WRAP()*PETSC_i;
280: } else SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid random number type");
281: #else
282: if (r->iset == PETSC_TRUE) *val = r->width * RAND_WRAP() + r->low;
283: else *val = RAND_WRAP();
284: #endif
285: return(0);
286: }
288: #else
289: /* Should put a simple, portable random number generator here? */
291: extern double drand48();
293: #undef __FUNCT__
295: int PetscRandomCreate(MPI_Comm comm,PetscRandomType type,PetscRandom *r)
296: {
297: PetscRandom rr;
298: char arch[10];
299: int ierr;
302: *r = 0;
303: if (type != RANDOM_DEFAULT) SETERRQ(PETSC_ERR_SUP,"Not for this random number type");
304: PetscHeaderCreate(rr,_p_PetscRandom,int,PETSC_RANDOM_COOKIE,type,"random",comm,PetscRandomDestroy,0);
305: PetscLogObjectCreate(rr);
306: *r = rr;
307: PetscGetArchType(arch,10);
308: PetscPrintf(comm,"PetscRandomCreate: Warning: Random number generator not set for machine %s; using fake random numbers.n",arch);
309: return(0);
310: }
312: #undef __FUNCT__
314: int PetscRandomGetValue(PetscRandom r,PetscScalar *val)
315: {
318: #if defined(PETSC_USE_COMPLEX)
319: *val = (0.5,0.5);
320: #else
321: *val = 0.5;
322: #endif
323: return(0);
324: }
325: #endif