1 | MODULE mt19937ar |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE mt19937ar *** |
---|
4 | !! Mersen Twister random generator |
---|
5 | !!====================================================================== |
---|
6 | !! |
---|
7 | !! ** Purpose : |
---|
8 | !! ** Method : |
---|
9 | !! References : This a a Fortran90 translation of mt19937ar.c: |
---|
10 | !! A C-program for MT19937, with initialization improved 2002/1/26. |
---|
11 | !! Coded by Takuji Nishimura and Makoto Matsumoto. |
---|
12 | !! |
---|
13 | !! Before using, initialize the state by using init_genrand(seed) |
---|
14 | !! or init_by_array(init_key, key_length). |
---|
15 | !! |
---|
16 | !! Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, |
---|
17 | !! All rights reserved. |
---|
18 | !! |
---|
19 | !! Redistribution and use in source and binary forms, with or without |
---|
20 | !! modification, are permitted provided that the following conditions |
---|
21 | !! are met: |
---|
22 | !! |
---|
23 | !! 1. Redistributions of source code must retain the above copyright |
---|
24 | !! notice, this list of conditions and the following disclaimer. |
---|
25 | !! |
---|
26 | !! 2. Redistributions in binary form must reproduce the above copyright |
---|
27 | !! notice, this list of conditions and the following disclaimer in the |
---|
28 | !! documentation and/or other materials provided with the distribution. |
---|
29 | !! |
---|
30 | !! 3. The names of its contributors may not be used to endorse or promote |
---|
31 | !! products derived from this software without specific prior written |
---|
32 | !! permission. |
---|
33 | !! |
---|
34 | !! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
---|
35 | !! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
---|
36 | !! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
---|
37 | !! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR |
---|
38 | !! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
---|
39 | !! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
---|
40 | !! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
---|
41 | !! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
---|
42 | !! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
---|
43 | !! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
---|
44 | !! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
---|
45 | !! |
---|
46 | !! |
---|
47 | !! Any feedback is very welcome. |
---|
48 | !! http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html |
---|
49 | !! email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) |
---|
50 | !! |
---|
51 | !! Until 2001/4/6, MT had been distributed under GNU Public License, |
---|
52 | !! but after 2001/4/6, we decided to let MT be used for any purpose, |
---|
53 | !! including commercial use. 2002-versions mt19937ar.c, mt19937ar-cok.c |
---|
54 | !! are considered to be usable freely. |
---|
55 | !! |
---|
56 | !! History: |
---|
57 | !! ! 97 - 02 (Makoto Matsumoto/Takuji Nishimura) original c version |
---|
58 | !! ! 10-03 (A. Vidard) F90 translation and NEMOTAM adaptation |
---|
59 | !-----------------------------------------------------------------------! |
---|
60 | USE par_kind, ONLY : wp |
---|
61 | USE dom_oce , ONLY : nproc |
---|
62 | USE lib_mpp , ONLY : get_unit |
---|
63 | |
---|
64 | PRIVATE |
---|
65 | INTEGER, PARAMETER :: jpn = 624 ! period parameters |
---|
66 | INTEGER, PARAMETER :: jpm = 397 |
---|
67 | |
---|
68 | INTEGER*8, PARAMETER :: jpa = Z'9908b0df' ! constant vector A |
---|
69 | INTEGER*8, PARAMETER :: jpmsb = Z'80000000' ! most significant w-r bit mask |
---|
70 | INTEGER*8, PARAMETER :: jplsb = Z'7fffffff' ! least significant r bit mask |
---|
71 | ! |
---|
72 | INTEGER*8, PARAMETER :: jpall = Z'ffffffff' ! all bits mask |
---|
73 | |
---|
74 | INTEGER*8, PARAMETER :: jptmp1= Z'9d2c5680' ! Tampering paramters |
---|
75 | INTEGER*8, PARAMETER :: jptmp2= Z'efc60000' |
---|
76 | |
---|
77 | INTEGER :: mti |
---|
78 | INTEGER*8, DIMENSION(0:jpn-1) :: mt |
---|
79 | LOGICAL :: l_mtinit = .FALSE. |
---|
80 | |
---|
81 | PUBLIC & |
---|
82 | & l_mtinit, & |
---|
83 | & mtrand_int31, & |
---|
84 | & mtrand_int32, & |
---|
85 | & mtrand_real1, & |
---|
86 | & mtrand_real2, & |
---|
87 | & mtrand_real3, & |
---|
88 | & mtrand_real53,& |
---|
89 | & init_by_array,& |
---|
90 | & init_mtrand, & |
---|
91 | & mtrand_seeddump, & |
---|
92 | & mtrand_seedread |
---|
93 | CONTAINS |
---|
94 | SUBROUTINE init_mtrand(ks) |
---|
95 | !----------------------------------------------------------------------- |
---|
96 | !! initialize mt(0:N-1) with a seed |
---|
97 | !----------------------------------------------------------------------- |
---|
98 | INTEGER :: ks |
---|
99 | ! |
---|
100 | IF (.NOT.l_mtinit) THEN |
---|
101 | mt(0) = IAND(INT(ks,kind=8),jpall) |
---|
102 | DO mti = 1, jpn-1 |
---|
103 | mt(mti) = 1812433253 * IEOR(mt(mti-1),ISHFT(mt(mti-1),-30)) + mti |
---|
104 | ! See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. |
---|
105 | ! In the previous versions, MSBs of the seed affect |
---|
106 | ! only MSBs of the array mt[]. |
---|
107 | ! modified by Makoto Matsumoto |
---|
108 | mt(mti) = IAND(mt(mti),jpall) |
---|
109 | ! for >32 bit machines |
---|
110 | END DO |
---|
111 | l_mtinit = .TRUE. |
---|
112 | ENDIF |
---|
113 | END SUBROUTINE init_mtrand |
---|
114 | |
---|
115 | SUBROUTINE init_by_array(kinit_key,key_length) |
---|
116 | !----------------------------------------------------------------------- |
---|
117 | ! initialize by an array with array-length |
---|
118 | ! init_key is the array for initializing keys |
---|
119 | ! key_length is its length |
---|
120 | !----------------------------------------------------------------------- |
---|
121 | INTEGER :: kinit_key(0:*) |
---|
122 | INTEGER :: key_length |
---|
123 | INTEGER :: ji,jj,jk |
---|
124 | ! |
---|
125 | CALL init_mtrand(19650218) |
---|
126 | ji = 1 |
---|
127 | jj = 0 |
---|
128 | DO jk = MAX(jpn,key_length),1,-1 |
---|
129 | mt(ji) = IEOR(mt(ji),IEOR(mt(ji-1),ISHFT(mt(ji-1),-30)) * 1664525) & |
---|
130 | & + kinit_key(jj) + jj |
---|
131 | mt(ji) = IAND(mt(ji),jpall) |
---|
132 | ji = ji + 1 |
---|
133 | jj = jj + 1 |
---|
134 | IF (ji .GE. jpn) THEN |
---|
135 | mt(0) = mt(jpn-1) |
---|
136 | ji = 1 |
---|
137 | ENDIF |
---|
138 | IF (jj .GE. key_length) THEN |
---|
139 | jj = 0 |
---|
140 | ENDIF |
---|
141 | END DO |
---|
142 | DO jk = jpn-1,1,-1 |
---|
143 | mt(ji) = IEOR(mt(ji),IEOR(mt(ji-1),ISHFT(mt(ji-1),-30)) * 1566083941) - ji ! non linear |
---|
144 | mt(ji) = IAND(mt(ji),jpall) ! for WORDSIZE > 32 machines |
---|
145 | ji = ji + 1 |
---|
146 | IF (ji .GE. jpn) THEN |
---|
147 | mt(0) = mt(jpn-1) |
---|
148 | ji = 1 |
---|
149 | ENDIF |
---|
150 | END DO |
---|
151 | mt(0) = jpmsb ! MSB is 1; assuring non-zero initial array |
---|
152 | ! |
---|
153 | END SUBROUTINE init_by_array |
---|
154 | FUNCTION mtrand_int32() |
---|
155 | !----------------------------------------------------------------------- |
---|
156 | ! generates a random number on [0,0xffffffff]-interval |
---|
157 | !----------------------------------------------------------------------- |
---|
158 | INTEGER :: mtrand_int32 |
---|
159 | INTEGER :: jk |
---|
160 | INTEGER*8 :: iy |
---|
161 | INTEGER*8, DIMENSION(0:1) :: mag01 |
---|
162 | ! |
---|
163 | mag01(0) = 0 |
---|
164 | mag01(1) = jpa ! mag01[x] = x * vector A for x=0,1 |
---|
165 | ! |
---|
166 | IF(.NOT. l_mtinit)THEN ! if init_mtrand() has not been called |
---|
167 | CALL init_mtrand(5489) ! a default initial seed is used |
---|
168 | ENDIF |
---|
169 | ! |
---|
170 | IF (mti .GE. jpn) THEN |
---|
171 | DO jk = 0,jpn-jpm-1 |
---|
172 | iy = IOR(IAND(mt(jk),jpmsb),IAND(mt(jk+1),jplsb)) |
---|
173 | mt(jk) = IEOR(IEOR(mt(jk+jpm),ISHFT(iy,-1)),mag01(IAND(iy,INT(1,kind=8)))) |
---|
174 | END DO |
---|
175 | DO jk = jpn-jpm,jpn-1-1 |
---|
176 | iy = IOR(IAND(mt(jk),jpmsb),IAND(mt(jk+1),jplsb)) |
---|
177 | mt(jk) = IEOR(IEOR(mt(jk+(jpm-jpn)),ISHFT(iy,-1)),mag01(IAND(iy,INT(1,kind=8)))) |
---|
178 | END DO |
---|
179 | iy = IOR(IAND(mt(jpn-1),jpmsb),IAND(mt(0),jplsb)) |
---|
180 | mt(jk) = IEOR(IEOR(mt(jpm-1),ISHFT(iy,-1)),mag01(IAND(iy,INT(1,kind=8)))) |
---|
181 | mti = 0 |
---|
182 | ENDIF |
---|
183 | ! |
---|
184 | iy = mt(mti) |
---|
185 | mti = mti + 1 |
---|
186 | ! |
---|
187 | iy = IEOR(iy,ISHFT(iy,-11)) |
---|
188 | iy = IEOR(iy,IAND(ISHFT(iy,7),jptmp1)) |
---|
189 | iy = IEOR(iy,IAND(ISHFT(iy,15),jptmp2)) |
---|
190 | iy = IEOR(iy,ISHFT(iy,-18)) |
---|
191 | ! |
---|
192 | mtrand_int32 = iy |
---|
193 | END FUNCTION mtrand_int32 |
---|
194 | FUNCTION mtrand_int31() |
---|
195 | !----------------------------------------------------------------------- |
---|
196 | ! generates a random number on [0,0x7fffffff]-interval |
---|
197 | !----------------------------------------------------------------------- |
---|
198 | INTEGER :: mtrand_int31 |
---|
199 | mtrand_int31 = INT(ISHFT(mtrand_int32(),-1)) |
---|
200 | END FUNCTION mtrand_int31 |
---|
201 | FUNCTION mtrand_real1() |
---|
202 | !----------------------------------------------------------------------- |
---|
203 | ! generates a random number on [0,1]-real-interval |
---|
204 | !----------------------------------------------------------------------- |
---|
205 | REAL(wp) :: mtrand_real1, zr |
---|
206 | zr = REAL(mtrand_int32(),wp) |
---|
207 | IF (zr.LT.0._wp) zr = zr + 2._wp**32 |
---|
208 | mtrand_real1 = zr / 4294967295._wp |
---|
209 | END FUNCTION mtrand_real1 |
---|
210 | FUNCTION mtrand_real2() |
---|
211 | !----------------------------------------------------------------------- |
---|
212 | ! generates a random number on [0,1)-real-interval |
---|
213 | !----------------------------------------------------------------------- |
---|
214 | REAL(wp) :: mtrand_real2,zr |
---|
215 | zr = REAL(mtrand_int32(),wp) |
---|
216 | IF (zr.LT.0._wp) zr = zr + 2._wp**32 |
---|
217 | mtrand_real2=zr / 4294967296._wp |
---|
218 | END FUNCTION mtrand_real2 |
---|
219 | FUNCTION mtrand_real3() |
---|
220 | !----------------------------------------------------------------------- |
---|
221 | ! generates a random number on (0,1)-real-interval |
---|
222 | !----------------------------------------------------------------------- |
---|
223 | REAL(wp) :: mtrand_real3, zr |
---|
224 | zr = REAL(mtrand_int32(),wp) |
---|
225 | IF (zr.LT.0._wp) zr = zr + 2._wp**32 |
---|
226 | mtrand_real3 = ( zr + 0.5_wp ) / 4294967296._wp |
---|
227 | END FUNCTION mtrand_real3 |
---|
228 | FUNCTION mtrand_res53() |
---|
229 | !----------------------------------------------------------------------- |
---|
230 | ! generates a random number on [0,1) with 53-bit resolution |
---|
231 | !----------------------------------------------------------------------- |
---|
232 | REAL(wp) :: mtrand_res53 |
---|
233 | REAL(wp) :: za,zb |
---|
234 | za = ISHFT(mtrand_int32(),-5) |
---|
235 | zb = ISHFT(mtrand_int32(),-6) |
---|
236 | IF (za.LT.0._wp) za = za + 2._wp**32 |
---|
237 | IF (zb.LT.0._wp) zb = zb + 2._wp**32 |
---|
238 | mtrand_res53 = ( za * 67108864._wp + zb ) / 9007199254740992._wp |
---|
239 | END FUNCTION mtrand_res53 |
---|
240 | |
---|
241 | |
---|
242 | SUBROUTINE mtrand_seeddump(cdfile) |
---|
243 | !----------------------------------------------------------------------- |
---|
244 | ! Dumps the seed in a restart file. |
---|
245 | !----------------------------------------------------------------------- |
---|
246 | CHARACTER(len=*) :: cdfile |
---|
247 | INTEGER :: iunit |
---|
248 | |
---|
249 | iunit = get_unit() |
---|
250 | OPEN( unit=iunit, file = TRIM(cdfile), form = 'unformatted', & |
---|
251 | & access = 'sequential', status = 'new' ) |
---|
252 | WRITE( iunit ) mt |
---|
253 | WRITE( iunit ) mti |
---|
254 | CLOSE( iunit ) |
---|
255 | |
---|
256 | END SUBROUTINE mtrand_seeddump |
---|
257 | |
---|
258 | SUBROUTINE mtrand_seedread(cdfile) |
---|
259 | !----------------------------------------------------------------------- |
---|
260 | ! Reads the seed from a restart file. |
---|
261 | !----------------------------------------------------------------------- |
---|
262 | CHARACTER(len=*) :: cdfile |
---|
263 | INTEGER :: iunit |
---|
264 | |
---|
265 | iunit = get_unit() |
---|
266 | OPEN( unit=iunit, file = TRIM(cdfile), form = 'unformatted', & |
---|
267 | & access = 'sequential', status = 'old' ) |
---|
268 | READ( iunit ) mt |
---|
269 | READ( iunit ) mti |
---|
270 | CLOSE( iunit ) |
---|
271 | l_mtinit = .TRUE. |
---|
272 | |
---|
273 | END SUBROUTINE mtrand_seedre |
---|
274 | |
---|
275 | END MODULE mt19937ar |
---|