-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathZBQ.f90
232 lines (225 loc) · 9.56 KB
/
ZBQ.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
!**********************************************************************
!******** FILE: randgen.f **************
!******** AUTHORS: Richard Chandler **************
!******** ([email protected]) **************
!******** Paul Northrop **************
!******** ([email protected]) **************
!******** LAST MODIFIED: 26/8/03 **************
!******** See file randgen.txt for details **************
!**********************************************************************
!* These subroutines have been extracted from the randgen.f file
!* obtained at http://www.ucl.ac.uk/~ucakarc/work/randgen.html
!* and translated from fixed form f77 to free form f90 language
!**********************************************************************
BLOCK DATA ZBQLBD01
!*
!* Initializes seed array etc. for random number generator.
!* The values below have themselves been generated using the
!* NAG generator.
!*
COMMON /ZBQL0001/ ZBQLIX,B,C
DOUBLE PRECISION ZBQLIX(43),B,C
INTEGER I
DATA (ZBQLIX(I),I=1,43) /8.001441D7,5.5321801D8, &
1.69570999D8,2.88589940D8,2.91581871D8,1.03842493D8,&
7.9952507D7,3.81202335D8,3.11575334D8,4.02878631D8, &
2.49757109D8,1.15192595D8,2.10629619D8,3.99952890D8,&
4.12280521D8,1.33873288D8,7.1345525D7,2.23467704D8, &
2.82934796D8,9.9756750D7,1.68564303D8,2.86817366D8, &
1.14310713D8,3.47045253D8,9.3762426D7 ,1.09670477D8,&
3.20029657D8,3.26369301D8,9.441177D6,3.53244738D8, &
2.44771580D8,1.59804337D8,2.07319904D8,3.37342907D8,&
3.75423178D8,7.0893571D7 ,4.26059785D8,3.95854390D8,&
2.0081010D7,5.9250059D7,1.62176640D8,3.20429173D8, &
2.63576576D8/
DATA B / 4.294967291D9 /
DATA C / 0.0D0 /
END
!**********************************************************************
!**********************************************************************
!**********************************************************************
SUBROUTINE ZBQLINI(SEED)
!**********************************************************************
!* To initialize the random number generator - either
!* repeatably or nonrepeatably. Need double precision
!* variables because integer storage can't handle the
!* numbers involved
!**********************************************************************
!* ARGUMENTS
!* =========
!* SEED (integer, input). User-input number which generates
!* elements of the array ZBQLIX, which is subsequently used
!* in the random number generation algorithm. If SEED=0,
!* the array is seeded using the system clock if the
!* FORTRAN implementation allows it.
!**********************************************************************
!* PARAMETERS
!* ==========
!* LFLNO (integer). Number of lowest file handle to try when
!* opening a temporary file to copy the system clock into.
!* Default is 80 to keep out of the way of any existing
!* open files (although the program keeps searching till
!* it finds an available handle). If this causes problems,
!* (which will only happen if handles 80 through 99 are
!* already in use), decrease the default value.
!**********************************************************************
INTEGER LFLNO
PARAMETER (LFLNO=80)
!**********************************************************************
!**********************************************************************
!* VARIABLES
!* =========
!* SEED See above
!* ZBQLIX Seed array for the random number generator. Defined
!* in ZBQLBD01
!* B,C Used in congruential initialisation of ZBQLIX
!* SS,MM,} System clock secs, mins, hours and days
!* HH,DD }
!* FILNO File handle used for temporary file
!* INIT Indicates whether generator has already been initialised
!*
INTEGER SEED,SS,MM,HH,DD,FILNO,I
INTEGER INIT
DOUBLE PRECISION ZBQLIX(43),B,C
DOUBLE PRECISION TMPVAR1,DSS,DMM,DHH,DDD
COMMON /ZBQL0001/ ZBQLIX,B,C
SAVE INIT
!*
!* Ensure we don't call this more than once in a program
!*
IF (INIT.GE.1) THEN
IF(INIT.EQ.1) THEN
WRITE(*,1)
INIT = 2
ENDIF
RETURN
ELSE
INIT = 1
ENDIF
!*
!* If SEED = 0, cat the contents of the clock into a file
!* and transform to obtain ZQBLIX(1), then use a congr.
!* algorithm to set remaining elements. Otherwise take
!* specified value of SEED.
!*
!*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
!*>>>>>>> NB FOR SYSTEMS WHICH DO NOT SUPPORT THE >>>>>>>>>>>>
!*>>>>>>> (NON-STANDARD) 'CALL SYSTEM' COMMAND, >>>>>>>>>>>>
!*>>>>>>> THIS WILL NOT WORK, AND THE FIRST CLAUSE >>>>>>>>>>>>
!*>>>>>>> OF THE FOLLOWING IF BLOCK SHOULD BE >>>>>>>>>>>>
!*>>>>>>> COMMENTED OUT. >>>>>>>>>>>>
!*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
IF (SEED.EQ.0) THEN
!*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
!*>>>>>>> COMMENT OUT FROM HERE IF YOU DON'T HAVE >>>>>>>>>>>>
!*>>>>>>> 'CALL SYSTEM' CAPABILITY ... >>>>>>>>>>>>
!*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
CALL SYSTEM(' date +%S%M%H%j > zbql1234.tmp')
!*
!* Try all file numbers for LFLNO to 999
!*
FILNO = LFLNO
10 OPEN(FILNO,FILE='zbql1234.tmp',ERR=11)
GOTO 12
11 FILNO = FILNO + 1
IF (FILNO.GT.999) THEN
WRITE(*,2)
RETURN
ENDIF
GOTO 10
12 READ(FILNO,'(3(I2),I3)') SS,MM,HH,DD
CLOSE(FILNO)
CALL SYSTEM('rm zbql1234.tmp')
DSS = DINT((DBLE(SS)/6.0D1) * B)
DMM = DINT((DBLE(MM)/6.0D1) * B)
DHH = DINT((DBLE(HH)/2.4D1) * B)
DDD = DINT((DBLE(DD)/3.65D2) * B)
TMPVAR1 = DMOD(DSS+DMM+DHH+DDD,B)
!*<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!*<<<<<<<< ... TO HERE (END OF COMMENTING OUT FOR <<<<<<<<<<<
!*<<<<<<<< USERS WITHOUT 'CALL SYSTEM' CAPABILITY <<<<<<<<<<<
!*<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
ELSE
TMPVAR1 = DMOD(DBLE(SEED),B)
ENDIF
ZBQLIX(1) = TMPVAR1
DO 100 I = 2,43
TMPVAR1 = ZBQLIX(I-1)*3.0269D4
TMPVAR1 = DMOD(TMPVAR1,B)
ZBQLIX(I) = TMPVAR1
100 CONTINUE
1 FORMAT(//5X,'****WARNING**** You have called routine ZBQLINI ',&
'more than',/5X,'once. I''m ignoring any subsequent calls.',//)
2 FORMAT(//5X,'**** ERROR **** In routine ZBQLINI, I couldn''t', &
' find an',/5X, &
'available file number. To rectify the problem, decrease the ',&
'value of',/5X, &
'the parameter LFLNO at the start of this routine (in file ', &
'randgen.f)',/5X, &
'and recompile. Any number less than 100 should work.')
END
!********************************************************************
FUNCTION ZBQLU01()
!********************************************************************
!* Returns a uniform random number between 0 & 1, using
!* a Marsaglia-Zaman type subtract-with-borrow generator.
!* Uses double precision, rather than integer, arithmetic
!* throughout because MZ's integer constants overflow
!* 32-bit integer storage (which goes from -2^31 to 2^31).
!* Ideally, we would explicitly truncate all integer
!* quantities at each stage to ensure that the double
!* precision representations do not accumulate approximation
!* error; however, on some machines the use of DNINT to
!* accomplish this is *seriously* slow (run-time increased
!* by a factor of about 3). This double precision version
!* has been tested against an integer implementation that
!* uses long integers (non-standard and, again, slow) -
!* the output was identical up to the 16th decimal place
!* after 10^10 calls, so we're probably OK ...
!*
DOUBLE PRECISION ZBQLU01,B,C,ZBQLIX(43),X,B2,BINV
INTEGER CURPOS,ID22,ID43
COMMON /ZBQL0001/ ZBQLIX,B,C
SAVE /ZBQL0001/
SAVE CURPOS,ID22,ID43
DATA CURPOS,ID22,ID43 /1,22,43/
B2 = B
BINV = 1.0D0/B
5 X = ZBQLIX(ID22) - ZBQLIX(ID43) - C
IF (X.LT.0.0D0) THEN
X = X + B
C = 1.0D0
ELSE
C = 0.0D0
ENDIF
ZBQLIX(ID43) = X
!*
!* Update array pointers. Do explicit check for bounds of each to
!* avoid expense of modular arithmetic. If one of them is 0 the
!* others
!* won't be
!*
CURPOS = CURPOS - 1
ID22 = ID22 - 1
ID43 = ID43 - 1
IF (CURPOS.EQ.0) THEN
CURPOS=43
ELSEIF (ID22.EQ.0) THEN
ID22 = 43
ELSEIF (ID43.EQ.0) THEN
ID43 = 43
ENDIF
!*
!* The integer arithmetic there can yield X=0, which can cause
!* problems in subsequent routines (e.g. ZBQLEXP). The problem
!* is simply that X is discrete whereas U is supposed to
!* be continuous - hence if X is 0, go back and generate another
!* X and return X/B^2 (etc.), which will be uniform on (0,1/B).
!*
IF (X.LT.BINV) THEN
B2 = B2*B
GOTO 5
ENDIF
ZBQLU01 = X/B2
END
!*************************************************