-
Notifications
You must be signed in to change notification settings - Fork 0
/
randoms.f90
199 lines (179 loc) · 6.02 KB
/
randoms.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
! Code converted using TO_F90 by Alan Miller
! Date: 2012-03-16 Time: 11:09:33
!> \file
!! Random numbers.
!!
!! \author Volker Blobel, University Hamburg, 2005-2009 (initial Fortran77 version)
!! \author Claus Kleinwort, DESY (maintenance and developement)
!!
!! \copyright
!! Copyright (c) 2009 - 2015 Deutsches Elektronen-Synchroton,
!! Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
!! This library is free software; you can redistribute it and/or modify
!! it under the terms of the GNU Library General Public License as
!! published by the Free Software Foundation; either version 2 of the
!! License, or (at your option) any later version. \n\n
!! This library is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU Library General Public License for more details. \n\n
!! You should have received a copy of the GNU Library General Public
!! License along with this program (see the file COPYING.LIB for more
!! details); if not, write to the Free Software Foundation, Inc.,
!! 675 Mass Ave, Cambridge, MA 02139, USA.
!!
!! Random number generators for Uniform and Normal distribution:
!!
!! URAN() for U(0,1)
!! GRAN() for N(0,1)
!> F.Gutbrod random number generator.
!!
!! Return N random numbers U(0,1) in array A(N).
!! Initialization by entry GBRVIN.
!!
!! \param[in] n number of requested random number
!! \param[out] a array of requested random number
SUBROUTINE gbrshi(n,a)
USE mpdef
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpi) :: ian
INTEGER(mpi) :: iboost
INTEGER(mpi) :: ic
INTEGER(mpi) :: idum
INTEGER(mpi) :: irotor
INTEGER(mpi) :: iseed
INTEGER(mpi) :: it
INTEGER(mpi) :: iwarm
INTEGER(mpi) :: j
INTEGER(mpi) :: jseed
INTEGER(mpi) :: jwarm
INTEGER(mpi) :: k
INTEGER(mpi) :: m
INTEGER(mpi) :: mbuff
INTEGER(mpi), INTENT(IN) :: n
REAL(mps), INTENT(OUT) :: a(*)
INTEGER(mpi), PARAMETER :: nb=511
INTEGER(mpi), PARAMETER :: ia=16807
INTEGER(mpi), PARAMETER :: im=2147483647
INTEGER(mpi), PARAMETER :: iq=127773
INTEGER(mpi), PARAMETER :: ir=2836
REAL(mps), PARAMETER :: aeps=1.0E-10
REAL(mps), PARAMETER :: scalin=4.6566125E-10
COMMON/ranbuf/mbuff(0:nb),ian,ic,iboost
INTEGER(mpi) :: istart
irotor(m,n)=IEOR(ior(ishft(m,17),ishft(m,-15)),n)
DATA istart/0/,iwarm/10/,iseed/4711/
IF(istart /= 0) GO TO 20
WRITE(*,*) ' Automatic GBRSHI initialization using:'
! initialize buffer
10 idum=iseed+9876543 ! prevent damage, if iseed=0
WRITE(*,*) ' ISEED=',iseed,' IWARM=',iwarm
DO j=0,nb+1 ! fill buffer
k=idum/iq ! minimal standard generator
idum=ia*(idum-k*iq)-ir*k ! with Schrages method
IF(idum < 0) idum=idum+im !
mbuff(j)=ishft(idum,1) ! fill in leading bit
END DO
ian=IAND(ian,nb) ! mask angle
ic=1 ! set pointer
iboost=0
DO j=1,iwarm*nb ! warm up a few times
it=mbuff(ian) ! hit ball angle
mbuff(ian)=irotor(it,ic) ! new spin
ic=it ! replace red spin
ian=IAND(it+iboost,nb) ! boost and mask angle
iboost=iboost+1 ! increment boost
END DO
IF(istart < 0) RETURN ! return for RBNVIN
istart=1 ! set done-flag
! generate array of r.n.
20 DO i=1,n
it=mbuff(ian) ! hit ball angle
mbuff(ian)=irotor(it,ic) ! new spin
ic=it ! replace red spin
ian=IAND(it+iboost,nb) ! boost and mask angle
a(i)=REAL(ishft(it,-1),mps)*scalin+aeps ! avoid zero output
iboost=iboost+1 ! increment boost
END DO
iboost=IAND(iboost,nb)
RETURN
ENTRY gbrvin(jseed,jwarm) ! initialize, but only once
IF(istart == 0) THEN
WRITE(*,*) ' Gbrshi initialization by GBRVIN-call using:'
iseed=jseed ! copy seed and
iwarm=jwarm ! warm-up parameter
istart=-1 ! start flag
GO TO 10
END IF
END SUBROUTINE gbrshi
!> GBRSHI initialization using TIME().
SUBROUTINE gbrtim
USE mpdef
IMPLICIT NONE
INTEGER(mpi) :: jseed
REAL(mps) :: time
LOGICAL :: done
DATA done/.FALSE./
IF(done) RETURN
jseed=time()
WRITE(*,*) ' Gbrshi initialialization using Time()'
CALL gbrvin(jseed,10)
done=.TRUE.
END SUBROUTINE gbrtim
!> Random number U(0,1) using RANSHI.
!!
!! \return random number U(0,1)
REAL(mps) FUNCTION uran() ! U(0,1)
USE mpdef
IMPLICIT NONE
INTEGER(mpi) :: indx
INTEGER(mpi) :: ndim
PARAMETER (ndim=100)
REAL(mps) :: buffer(ndim)
DATA indx/ndim/
SAVE indx,buffer
indx=MOD(indx,ndim)+1
IF(indx == 1) CALL gbrshi(ndim,buffer)
uran=buffer(indx)
END FUNCTION uran
!> Gauss random number.
!!
!! \return random number N(0,1)
REAL(mps) FUNCTION gran() ! N(0,1)
USE mpdef
IMPLICIT NONE
REAL(mps) :: al
REAL(mps) :: cs
INTEGER(mpi) :: indx
INTEGER(mpi) :: kn
INTEGER(mpi) :: ndim
REAL(mps) :: radsq
REAL(mps) :: rn1
REAL(mps) :: rn2
REAL(mps) :: sn
PARAMETER (ndim=100)
REAL(mps) :: buffer(ndim)
DATA indx/ndim/,kn/1/
SAVE indx,buffer,kn,cs,al
! ...
IF(kn <= 1) THEN
! two U(-1,+1) random numbers
10 indx=MOD(indx,ndim)+2
IF(indx == 2) CALL gbrshi(ndim,buffer)
rn1=buffer(indx-1)-1.0+buffer(indx-1)
rn2=buffer(indx )-1.0+buffer(indx)
radsq=rn1*rn1+rn2*rn2
IF(radsq > 1.0) GO TO 10 ! test point inside circle?
! sine and cosine for random phi
sn=rn1/SQRT(radsq)
cs=rn2/SQRT(radsq)
! transform to gaussians
al=SQRT(-2.0*LOG(radsq))
kn =2
gran=sn*al
ELSE
kn =1
gran=cs*al
END IF
END FUNCTION gran