NCEPLIBS-w3emc  2.11.0
mersenne_twister.f
Go to the documentation of this file.
1 
4 
90  module mersenne_twister
91  private
92 ! Public declarations
93  public random_stat
94  public random_seed
95  public random_setseed
96  public random_number
97  public random_number_f
98  public random_gauss
99  public random_gauss_f
100  public random_index
101  public random_index_f
102 ! Parameters
103  integer,parameter:: n=624
104  integer,parameter:: m=397
105  integer,parameter:: mata=-1727483681 ! constant vector a
106  integer,parameter:: umask=-2147483648 ! most significant w-r bits
107  integer,parameter:: lmask =2147483647 ! least significant r bits
108  integer,parameter:: tmaskb=-1658038656 ! tempering parameter
109  integer,parameter:: tmaskc=-272236544 ! tempering parameter
110  integer,parameter:: mag01(0:1)=(/0,mata/)
111  integer,parameter:: iseed=4357
112  integer,parameter:: nrest=n+6
113 ! Defined types
114  type random_stat
115  private
116  integer:: mti=n+1
117  integer:: mt(0:n-1)
118  integer:: iset
119  real:: gset
120  end type
121 ! Saved data
122  type(random_stat),save:: sstat
123 ! Overloaded interfaces
124  interface random_setseed
125  module procedure random_setseed_s
126  module procedure random_setseed_t
127  end interface
128  interface random_number
129  module procedure random_number_i
130  module procedure random_number_s
131  module procedure random_number_t
132  end interface
133  interface random_gauss
134  module procedure random_gauss_i
135  module procedure random_gauss_s
136  module procedure random_gauss_t
137  end interface
138  interface random_index
139  module procedure random_index_i
140  module procedure random_index_s
141  module procedure random_index_t
142  end interface
143 ! All the subprograms
144  contains
150  subroutine random_seed(size,put,get,stat)
151  implicit none
152  integer,intent(out),optional:: size
153  integer,intent(in),optional:: put(nrest)
154  integer,intent(out),optional:: get(nrest)
155  type(random_stat),intent(inout),optional:: stat
156  if(present(size)) then ! return size of seed array
157 ! if(present(put).or.present(get))&
158 ! call errmsg('RANDOM_SEED: more than one option set - some ignored')
159  size=nrest
160  elseif(present(put)) then ! restore from seed array
161 ! if(present(get))&
162 ! call errmsg('RANDOM_SEED: more than one option set - some ignored')
163  if(present(stat)) then
164  stat%mti=put(1)
165  stat%mt=put(2:n+1)
166  stat%iset=put(n+2)
167  stat%gset=transfer(put(n+3:nrest),stat%gset)
168  if(stat%mti.lt.0.or.stat%mti.gt.n.or.any(stat%mt.eq.0).or.
169  & stat%iset.lt.0.or.stat%iset.gt.1) then
170  call random_setseed_t(iseed,stat)
171 ! call errmsg('RANDOM_SEED: invalid seeds put - default seeds used')
172  endif
173  else
174  sstat%mti=put(1)
175  sstat%mt=put(2:n+1)
176  sstat%iset=put(n+2)
177  sstat%gset=transfer(put(n+3:nrest),sstat%gset)
178  if(sstat%mti.lt.0.or.sstat%mti.gt.n.or.any(sstat%mt.eq.0)
179  & .or.sstat%iset.lt.0.or.sstat%iset.gt.1) then
180  call random_setseed_t(iseed,sstat)
181 ! call errmsg('RANDOM_SEED: invalid seeds put - default seeds used')
182  endif
183  endif
184  elseif(present(get)) then ! save to seed array
185  if(present(stat)) then
186  if(stat%mti.eq.n+1) call random_setseed_t(iseed,stat)
187  get(1)=stat%mti
188  get(2:n+1)=stat%mt
189  get(n+2)=stat%iset
190  get(n+3:nrest)=transfer(stat%gset,get,nrest-(n+3)+1)
191  else
192  if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
193  get(1)=sstat%mti
194  get(2:n+1)=sstat%mt
195  get(n+2)=sstat%iset
196  get(n+3:nrest)=transfer(sstat%gset,get,nrest-(n+3)+1)
197  endif
198  else ! reset default seed
199  if(present(stat)) then
200  call random_setseed_t(iseed,stat)
201  else
202  call random_setseed_t(iseed,sstat)
203  endif
204  endif
205  end subroutine
208  subroutine random_setseed_s(inseed)
209  implicit none
210  integer,intent(in):: inseed
211  call random_setseed_t(inseed,sstat)
212  end subroutine
216  subroutine random_setseed_t(inseed,stat)
217  implicit none
218  integer,intent(in):: inseed
219  type(random_stat),intent(out):: stat
220  integer ii,mti
221  ii=inseed
222  if(ii.eq.0) ii=iseed
223  stat%mti=n
224  stat%mt(0)=iand(ii,-1)
225  do mti=1,n-1
226  stat%mt(mti)=iand(69069*stat%mt(mti-1),-1)
227  enddo
228  stat%iset=0
229  stat%gset=0.
230  end subroutine
233  function random_number_f() result(harvest)
234  implicit none
235  real:: harvest
236  real h(1)
237  if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
238  call random_number_t(h,sstat)
239  harvest=h(1)
240  end function
244  subroutine random_number_i(harvest,inseed)
245  implicit none
246  real,intent(out):: harvest(:)
247  integer,intent(in):: inseed
248  type(random_stat) stat
249  call random_setseed_t(inseed,stat)
250  call random_number_t(harvest,stat)
251  end subroutine
254  subroutine random_number_s(harvest)
255  implicit none
256  real,intent(out):: harvest(:)
257  if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
258  call random_number_t(harvest,sstat)
259  end subroutine
263  subroutine random_number_t(harvest,stat)
264  implicit none
265  real,intent(out):: harvest(:)
266  type(random_stat),intent(inout):: stat
267  integer j,kk,y
268  integer tshftu,tshfts,tshftt,tshftl
269  tshftu(y)=ishft(y,-11)
270  tshfts(y)=ishft(y,7)
271  tshftt(y)=ishft(y,15)
272  tshftl(y)=ishft(y,-18)
273  do j=1,size(harvest)
274  if(stat%mti.ge.n) then
275  do kk=0,n-m-1
276  y=ior(iand(stat%mt(kk),umask),iand(stat%mt(kk+1),lmask))
277  stat%mt(kk)=ieor(ieor(stat%mt(kk+m),ishft(y,-1)),
278  & mag01(iand(y,1)))
279  enddo
280  do kk=n-m,n-2
281  y=ior(iand(stat%mt(kk),umask),iand(stat%mt(kk+1),lmask))
282  stat%mt(kk)=ieor(ieor(stat%mt(kk+(m-n)),ishft(y,-1)),
283  & mag01(iand(y,1)))
284  enddo
285  y=ior(iand(stat%mt(n-1),umask),iand(stat%mt(0),lmask))
286  stat%mt(n-1)=ieor(ieor(stat%mt(m-1),ishft(y,-1)),
287  & mag01(iand(y,1)))
288  stat%mti=0
289  endif
290  y=stat%mt(stat%mti)
291  y=ieor(y,tshftu(y))
292  y=ieor(y,iand(tshfts(y),tmaskb))
293  y=ieor(y,iand(tshftt(y),tmaskc))
294  y=ieor(y,tshftl(y))
295  if(y.lt.0) then
296  harvest(j)=(real(y)+2.0**32)/(2.0**32-1.0)
297  else
298  harvest(j)=real(y)/(2.0**32-1.0)
299  endif
300  stat%mti=stat%mti+1
301  enddo
302  end subroutine
305  function random_gauss_f() result(harvest)
306  implicit none
307  real:: harvest
308  real h(1)
309  if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
310  call random_gauss_t(h,sstat)
311  harvest=h(1)
312  end function
316  subroutine random_gauss_i(harvest,inseed)
317  implicit none
318  real,intent(out):: harvest(:)
319  integer,intent(in):: inseed
320  type(random_stat) stat
321  call random_setseed_t(inseed,stat)
322  call random_gauss_t(harvest,stat)
323  end subroutine
326  subroutine random_gauss_s(harvest)
327  implicit none
328  real,intent(out):: harvest(:)
329  if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
330  call random_gauss_t(harvest,sstat)
331  end subroutine
335  subroutine random_gauss_t(harvest,stat)
336  implicit none
337  real,intent(out):: harvest(:)
338  type(random_stat),intent(inout):: stat
339  integer mx,my,mz,j
340  real r2(2),r,g1,g2
341  mz=size(harvest)
342  if(mz.le.0) return
343  mx=0
344  if(stat%iset.eq.1) then
345  mx=1
346  harvest(1)=stat%gset
347  stat%iset=0
348  endif
349  my=(mz-mx)/2*2+mx
350  do
351  call random_number_t(harvest(mx+1:my),stat)
352  do j=mx,my-2,2
353  call rgauss(harvest(j+1),harvest(j+2),r,g1,g2)
354  if(r.lt.1.) then
355  harvest(mx+1)=g1
356  harvest(mx+2)=g2
357  mx=mx+2
358  endif
359  enddo
360  if(mx.eq.my) exit
361  enddo
362  if(my.lt.mz) then
363  do
364  call random_number_t(r2,stat)
365  call rgauss(r2(1),r2(2),r,g1,g2)
366  if(r.lt.1.) exit
367  enddo
368  harvest(mz)=g1
369  stat%gset=g2
370  stat%iset=1
371  endif
372  contains
379  subroutine rgauss(r1,r2,r,g1,g2)
380  real,intent(in):: r1,r2
381  real,intent(out):: r,g1,g2
382  real v1,v2,fac
383  v1=2.*r1-1.
384  v2=2.*r2-1.
385  r=v1**2+v2**2
386  if(r.lt.1.) then
387  fac=sqrt(-2.*log(r)/r)
388  g1=v1*fac
389  g2=v2*fac
390  endif
391  end subroutine
392  end subroutine
396  function random_index_f(imax) result(iharvest)
397  implicit none
398  integer,intent(in):: imax
399  integer:: iharvest
400  integer ih(1)
401  if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
402  call random_index_t(imax,ih,sstat)
403  iharvest=ih(1)
404  end function
409  subroutine random_index_i(imax,iharvest,inseed)
410  implicit none
411  integer,intent(in):: imax
412  integer,intent(out):: iharvest(:)
413  integer,intent(in):: inseed
414  type(random_stat) stat
415  call random_setseed_t(inseed,stat)
416  call random_index_t(imax,iharvest,stat)
417  end subroutine
421  subroutine random_index_s(imax,iharvest)
422  implicit none
423  integer,intent(in):: imax
424  integer,intent(out):: iharvest(:)
425  if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
426  call random_index_t(imax,iharvest,sstat)
427  end subroutine
432  subroutine random_index_t(imax,iharvest,stat)
433  implicit none
434  integer,intent(in):: imax
435  integer,intent(out):: iharvest(:)
436  type(random_stat),intent(inout):: stat
437  integer,parameter:: mh=n
438  integer i1,i2,mz
439  real h(mh)
440  mz=size(iharvest)
441  do i1=1,mz,mh
442  i2=min((i1-1)+mh,mz)
443  call random_number_t(h(:i2-(i1-1)),stat)
444  iharvest(i1:i2)=max(ceiling(h(:i2-(i1-1))*imax),1)
445  enddo
446  end subroutine
447  end module
mersenne_twister
This module calculates random numbers using the Mersenne twister.
Definition: mersenne_twister.f:90
mersenne_twister::random_number_f
real function, public random_number_f()
Generates random numbers in functional mode.
Definition: mersenne_twister.f:234
mersenne_twister::random_seed
subroutine, public random_seed(size, put, get, stat)
Sets and gets state; overloads Fortran 90 standard.
Definition: mersenne_twister.f:151
mersenne_twister::random_index_f
integer function, public random_index_f(imax)
Generates random indices in functional mode.
Definition: mersenne_twister.f:397
mersenne_twister::random_gauss_f
real function, public random_gauss_f()
Generates Gaussian random numbers in functional mode.
Definition: mersenne_twister.f:306