Branch data Line data Source code
1 : : !!****m* ABINIT/m_random_zbq
2 : : !! NAME
3 : : !! m_random_zbq
4 : : !!
5 : : !! FUNCTION
6 : : !! This library contains routines for generating uniform random numbers
7 : : !! between 0 & 1, using a Marsaglia-Zaman type subtract-with-borrow generator.
8 : : !! This random number generator is non deterministic (it is not "pseudo-random").
9 : : !!
10 : : !! These routines have been downloaded from Richard Chandler homepage,
11 : : !! then converted into Fortran90 standard.
12 : : !! Authors:
13 : : !! Richard Chandler (email: richard@stats.ucl.ac.uk)
14 : : !! Paul Northrop (email: paul@stats.ucl.ac.uk)
15 : : !!
16 : : !! COPYRIGHT
17 : : !! Copyright (C) Richard Chandler/Paul Northrop
18 : : !! Please feel free to use and adapt these routines for your own use. We
19 : : !! do, however, request that if the routines are used in work which is later
20 : : !! presented in public, or published, the source of the code (ie. us!) is
21 : : !! acknowledged. If for no other reason, this will allow other people to test
22 : : !! the routines to see if there's anything horribly wrong with them that will
23 : : !! invalidate your results.
24 : : !!
25 : : !! See http://www.ucl.ac.uk/~ucakarc/work/software/randgen.txt
26 : : !! for further information
27 : : !!
28 : : !! SOURCE
29 : :
30 : : #if defined HAVE_CONFIG_H
31 : : #include "config.h"
32 : : #endif
33 : :
34 : : #include "abi_common.h"
35 : :
36 : :
37 : : module m_random_zbq
38 : :
39 : : use m_abicore
40 : :
41 : : use defs_basis, only : sp, dp, std_out
42 : :
43 : : implicit none
44 : :
45 : : private
46 : : !!***
47 : :
48 : : public :: ZBQLU01
49 : : public :: ZBQLINI
50 : :
51 : : REAL(DP),SAVE,PRIVATE :: B=4.294967291D9,C=0.0D0
52 : : REAL(DP),SAVE,PRIVATE :: ZBQLIX(43)= &
53 : : & (/8.001441D7,5.5321801D8,&
54 : : & 1.69570999D8,2.88589940D8,2.91581871D8,1.03842493D8,&
55 : : & 7.9952507D7,3.81202335D8,3.11575334D8,4.02878631D8, &
56 : : & 2.49757109D8,1.15192595D8,2.10629619D8,3.99952890D8,&
57 : : & 4.12280521D8,1.33873288D8,7.1345525D7,2.23467704D8, &
58 : : & 2.82934796D8,9.9756750D7,1.68564303D8,2.86817366D8, &
59 : : & 1.14310713D8,3.47045253D8,9.3762426D7 ,1.09670477D8,&
60 : : & 3.20029657D8,3.26369301D8,9.441177D6,3.53244738D8, &
61 : : & 2.44771580D8,1.59804337D8,2.07319904D8,3.37342907D8,&
62 : : & 3.75423178D8,7.0893571D7 ,4.26059785D8,3.95854390D8,&
63 : : & 2.0081010D7,5.9250059D7,1.62176640D8,3.20429173D8, &
64 : : & 2.63576576D8/)
65 : :
66 : :
67 : : contains
68 : :
69 : : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
70 : :
71 : : !*******************************************************************
72 : 0 : REAL(DP) FUNCTION ZBQLU01(DUMMY)
73 : : !*******************************************************************
74 : : ! Returns a uniform random number between 0 & 1, using
75 : : ! a Marsaglia-Zaman type subtract-with-borrow generator.
76 : : ! Uses double precision, rather than integer, arithmetic
77 : : ! throughout because MZ's integer constants overflow
78 : : ! 32-bit integer storage (which goes from -2^31 to 2^31).
79 : : ! Ideally, we would explicitly truncate all integer
80 : : ! quantities at each stage to ensure that the double
81 : : ! precision representations do not accumulate approximation
82 : : ! error; however, on some machines the use of DNINT to
83 : : ! accomplish this is *seriously* slow (run-time increased
84 : : ! by a factor of about 3). This double precision version
85 : : ! has been tested against an integer implementation that
86 : : ! uses long integers (non-standard and, again, slow) -
87 : : ! the output was identical up to the 16th decimal place
88 : : ! after 10^10 calls, so we're probably OK ...
89 : : !
90 : : REAL(DP) :: DUMMY
91 : :
92 : : INTEGER,SAVE :: CURPOS=1,ID22=22,ID43=43
93 : : REAL(DP) :: X,B2,BINV
94 : :
95 : 0 : B2 = B
96 : 0 : BINV = 1.0D0/B
97 : 0 : 5 X = ZBQLIX(ID22) - ZBQLIX(ID43) - C
98 [ # # ]: 0 : IF (X<0.0D0) THEN
99 : 0 : X = X + B
100 : 0 : C = 1.0D0
101 : : ELSE
102 : 0 : C = 0.0D0
103 : : ENDIF
104 : 0 : ZBQLIX(ID43) = X
105 : : !
106 : : ! Update array pointers. Do explicit check for bounds of each to
107 : : ! avoid expense of modular arithmetic. If one of them is 0 the others
108 : : ! won't be
109 : : !
110 : 0 : CURPOS = CURPOS - 1
111 : 0 : ID22 = ID22 - 1
112 : 0 : ID43 = ID43 - 1
113 [ # # ]: 0 : IF (CURPOS==0) THEN
114 : 0 : CURPOS=43
115 [ # # ]: 0 : ELSEIF (ID22==0) THEN
116 : 0 : ID22 = 43
117 [ # # ]: 0 : ELSEIF (ID43==0) THEN
118 : 0 : ID43 = 43
119 : : ENDIF
120 : : !
121 : : ! The integer arithmetic there can yield X=0, which can cause
122 : : ! problems in subsequent routines (e.g. ZBQLEXP). The problem
123 : : ! is simply that X is discrete whereas U is supposed to
124 : : ! be continuous - hence if X is 0, go back and generate another
125 : : ! X and return X/B^2 (etc.), which will be uniform on (0,1/B).
126 [ # # ]: 0 : IF (X<BINV) THEN
127 : 0 : B2 = B2*B
128 : 0 : GOTO 5
129 : : ENDIF
130 : :
131 : 0 : ZBQLU01 = X/B2
132 : :
133 : : RETURN;X=DUMMY
134 : :
135 : : END FUNCTION ZBQLU01
136 : :
137 : : !*********************************************************************
138 : 0 : SUBROUTINE ZBQLINI(SEED)
139 : : !*********************************************************************
140 : : ! To initialize the random number generator - either
141 : : ! repeatably or nonrepeatably. Need double precision
142 : : ! variables because integer storage can't handle the
143 : : ! numbers involved
144 : : !*********************************************************************
145 : : ! ARGUMENTS
146 : : ! =========
147 : : ! SEED (integer, input). User-input number which generates
148 : : ! elements of the array ZBQLIX, which is subsequently used
149 : : ! in the random number generation algorithm. If SEED=0,
150 : : ! the array is seeded using the system clock if the
151 : : ! FORTRAN implementation allows it.
152 : : !*********************************************************************
153 : : INTEGER :: SEED
154 : : !*********************************************************************
155 : : ! PARAMETERS
156 : : ! ==========
157 : : ! LFLNO (integer). Number of lowest file handle to try when
158 : : ! opening a temporary file to copy the system clock into.
159 : : ! Default is 80 to keep out of the way of any existing
160 : : ! open files (although the program keeps searching till
161 : : ! it finds an available handle). If this causes problems,
162 : : ! (which will only happen if handles 80 through 99 are
163 : : ! already in use), decrease the default value.
164 : : !*********************************************************************
165 : : ! INTEGER,PARAMETER :: LFLNO=80
166 : : !*********************************************************************
167 : : ! VARIABLES
168 : : ! =========
169 : : ! SEED See above
170 : : ! ZBQLIX Seed array for the random number generator. Defined
171 : : ! in ZBQLBD01
172 : : ! B,C Used in congruential initialisation of ZBQLIX
173 : : ! SS,MM,} System clock secs, mins, hours and days
174 : : ! HH,DD }
175 : : ! FILNO File handle used for temporary file
176 : : ! INIT Indicates whether generator has already been initialised
177 : :
178 : : INTEGER,SAVE :: INIT
179 : : INTEGER :: SS,MM,HH,DD,I
180 : : REAL(DP) :: TMPVAR1,DSS,DMM,DHH,DDD
181 : : ! Variable used by date_and_time function:
182 : : CHARACTER(LEN=8) :: DATE
183 : : CHARACTER(LEN=10) :: TIME
184 : : CHARACTER(LEN=5) :: ZONE
185 : : INTEGER :: VALUES(8)
186 : :
187 : : ! Ensure we don't call this more than once in a program
188 [ # # ]: 0 : IF (INIT>=1) THEN
189 [ # # ]: 0 : IF(INIT==1) THEN
190 : 0 : WRITE(std_out,1)
191 : : 1 FORMAT(//5X,'****WARNING**** You have called routine ZBQLINI ',&
192 : : & 'more than',/5X,'once. I''m ignoring any subsequent calls.',//)
193 : 0 : INIT = 2
194 : : END IF
195 : 0 : RETURN
196 : : ELSE
197 : 0 : INIT = 1
198 : : ENDIF
199 : :
200 : : ! If SEED = 0, cat the contents of the clock into a file
201 : : ! and transform to obtain ZQBLIX(1), then use a congr.
202 : : ! algorithm to set remaining elements. Otherwise take
203 : : ! specified value of SEED.
204 : : !
205 [ # # ]: 0 : IF (SEED==0) THEN
206 : : ! Initial coding:
207 : : !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
208 : : !>>>>>>> COMMENT OUT FROM HERE IF YOU DON'T HAVE >>>>>>>>>>>>
209 : : !>>>>>>> 'CALL SYSTEM' CAPABILITY ... >>>>>>>>>>>>
210 : : !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
211 : : ! CALL SYSTEM(' date +%S%M%H%j > zbql1234.tmp')
212 : : ! ! Try all file numbers for LFLNO to 999
213 : : ! FILNO = LFLNO
214 : : ! 10 OPEN(FILNO,FILE='zbql1234.tmp',ERR=11)
215 : : ! GOTO 12
216 : : ! 11 FILNO = FILNO + 1
217 : : ! IF (FILNO.GT.999) THEN
218 : : ! WRITE(std_out,*) 'ZBQLINI: Error on temporary file !'
219 : : ! RETURN
220 : : ! ENDIF
221 : : ! GOTO 10
222 : : ! 12 READ(FILNO,'(3(I2),I3)') SS,MM,HH,DD
223 : : ! CLOSE(FILNO)
224 : : ! CALL SYSTEM('rm zbql1234.tmp')
225 : :
226 : : ! New F90 coding:
227 : 0 : CALL DATE_AND_TIME(DATE,TIME,ZONE,VALUES)
228 : 0 : HH=VALUES(4)
229 : 0 : MM=VALUES(5)
230 : 0 : SS=VALUES(6)
231 : 0 : DD=(VALUES(2)-1)*30+VALUES(3) ! Approx. for the day index
232 : :
233 : 0 : DSS = DINT((DBLE(SS)/6.0D1) * B)
234 : 0 : DMM = DINT((DBLE(MM)/6.0D1) * B)
235 : 0 : DHH = DINT((DBLE(HH)/2.4D1) * B)
236 : 0 : DDD = DINT((DBLE(DD)/3.65D2) * B)
237 : 0 : TMPVAR1 = DMOD(DSS+DMM+DHH+DDD,B)
238 : : !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
239 : : !<<<<<<<< ... TO HERE (END OF COMMENTING OUT FOR <<<<<<<<<<<
240 : : !<<<<<<<< USERS WITHOUT 'CALL SYSTEM' CAPABILITY <<<<<<<<<<<
241 : : !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
242 : : ELSE
243 : 0 : TMPVAR1 = DMOD(DBLE(SEED),B)
244 : : ENDIF
245 : :
246 : 0 : ZBQLIX(1) = TMPVAR1
247 [ # # ]: 0 : DO I = 2,43
248 : 0 : TMPVAR1 = ZBQLIX(I-1)*3.0269D4
249 : 0 : TMPVAR1 = DMOD(TMPVAR1,B)
250 : 0 : ZBQLIX(I) = TMPVAR1
251 : : END DO
252 : :
253 : : END SUBROUTINE ZBQLINI
254 : :
255 : : end module m_random_zbq
256 : : !!***
|
