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 : !!***
|
