← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.2-alpha0

« back to all changes in this revision

Viewing changes to Template/Source/ranmar.f

Added Template and HELAS into bzr

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/Source/ranmar.f
 
1
      subroutine ntuple(x,a,b,ii,jconfig)
 
2
c-------------------------------------------------------
 
3
c     Front to ranmar which allows user to easily
 
4
c     choose the seed.
 
5
c------------------------------------------------------
 
6
      implicit none
 
7
c
 
8
c     Arguments
 
9
c
 
10
      double precision x,a,b
 
11
      integer ii,jconfig
 
12
c
 
13
c     Local
 
14
c
 
15
      integer init, ioffset
 
16
      integer     ij, kl, iseed1,iseed2
 
17
 
 
18
c
 
19
c     Global
 
20
c
 
21
      integer         iseed
 
22
      common /to_seed/iseed
 
23
c
 
24
c     Data
 
25
c
 
26
      data init /1/
 
27
      save ij, kl
 
28
c-----
 
29
c  Begin Code
 
30
c-----
 
31
      if (init .eq. 1) then
 
32
         init = 0
 
33
         call get_offset(ioffset)
 
34
         if (iseed .eq. 0) call get_base(iseed)
 
35
c
 
36
c     TJS 3/13/2008
 
37
c     Modified to allow for more sequences 
 
38
c     iseed can be between 0 and 31328*30081
 
39
c     before pattern repeats
 
40
c
 
41
         ij=1802+jconfig + mod(iseed,30081)
 
42
         kl=9373+(iseed/31328)+ioffset 
 
43
         write(*,'(a,i6,a3,i6)') 'Using random seed offsets',jconfig," : ",ioffset
 
44
         write(*,*) ' with seed', iseed
 
45
         do while (ij .gt. 31328)
 
46
            ij = ij - 31328
 
47
         enddo
 
48
         do while (kl .gt. 30081)
 
49
            kl = kl - 30081
 
50
         enddo
 
51
        call rmarin(ij,kl)         
 
52
      endif
 
53
      call ranmar(x)
 
54
      do while (x .lt. 1d-16)
 
55
         call ranmar(x)
 
56
      enddo
 
57
      x = a+x*(b-a)
 
58
      end
 
59
 
 
60
      subroutine get_base(iseed)
 
61
c-------------------------------------------------------
 
62
c     Looks for file iproc.dat to offset random number gen
 
63
c------------------------------------------------------
 
64
      implicit none
 
65
c
 
66
c     Constants
 
67
c
 
68
      integer    lun
 
69
      parameter (lun=22)
 
70
c
 
71
c     Arguments
 
72
c
 
73
      integer iseed
 
74
c
 
75
c     Local
 
76
c
 
77
      character*60 fname
 
78
      logical done
 
79
      integer i,level
 
80
c-----
 
81
c  Begin Code
 
82
c-----
 
83
 
 
84
      fname = 'randinit'
 
85
      done = .false.
 
86
      level = 1
 
87
      do while(.not. done .and. level .lt. 5)
 
88
         open(unit=lun,file=fname,status='old',err=15)
 
89
         done = .true.
 
90
 15      level = level+1
 
91
         fname = '../' // fname
 
92
         i=index(fname,' ')
 
93
         if (i .gt. 0) fname=fname(1:i-1)
 
94
      enddo
 
95
      if (done) then
 
96
         read(lun,'(a)',end=24,err=24) fname
 
97
         i = index(fname,'=')
 
98
         if (i .gt. 0) fname=fname(i+1:)
 
99
         read(fname,*,err=26,end=26) iseed
 
100
 24      close(lun)
 
101
c         write(*,*) 'Read iseed from randinit ',iseed
 
102
         return
 
103
 26      close(lun)
 
104
      endif
 
105
 25   iseed = 0
 
106
c      write(*,*) 'No base found using iseed=0'
 
107
      end
 
108
 
 
109
      subroutine get_offset(iseed)
 
110
c-------------------------------------------------------
 
111
c     Looks for file iproc.dat to offset random number gen
 
112
c------------------------------------------------------
 
113
      implicit none
 
114
c
 
115
c     Constants
 
116
c
 
117
      integer    lun
 
118
      parameter (lun=22)
 
119
c
 
120
c     Arguments
 
121
c
 
122
      integer iseed
 
123
c
 
124
c     Local
 
125
c
 
126
c-----
 
127
c  Begin Code
 
128
c-----
 
129
 
 
130
      open(unit=lun,file='./iproc.dat',status='old',err=15)
 
131
         read(lun,*,err=14) iseed
 
132
         close(lun)
 
133
         return
 
134
 14   close(lun)
 
135
 15   open(unit=lun,file='../iproc.dat',status='old',err=25)
 
136
         read(lun,*,err=24) iseed
 
137
         close(lun)
 
138
         return
 
139
 24   close(lun)
 
140
 25   iseed = 0
 
141
      end
 
142
 
 
143
      subroutine ranmar(rvec)
 
144
*     -----------------
 
145
* universal random number generator proposed by marsaglia and zaman
 
146
* in report fsu-scri-87-50
 
147
* in this version rvec is a double precision variable.
 
148
      implicit real*8(a-h,o-z)
 
149
      common/ raset1 / ranu(97),ranc,rancd,rancm
 
150
      common/ raset2 / iranmr,jranmr
 
151
      save /raset1/,/raset2/
 
152
      uni = ranu(iranmr) - ranu(jranmr)
 
153
      if(uni .lt. 0d0) uni = uni + 1d0
 
154
      ranu(iranmr) = uni
 
155
      iranmr = iranmr - 1
 
156
      jranmr = jranmr - 1
 
157
      if(iranmr .eq. 0) iranmr = 97
 
158
      if(jranmr .eq. 0) jranmr = 97
 
159
      ranc = ranc - rancd
 
160
      if(ranc .lt. 0d0) ranc = ranc + rancm
 
161
      uni = uni - ranc
 
162
      if(uni .lt. 0d0) uni = uni + 1d0
 
163
      rvec = uni
 
164
      end
 
165
 
 
166
      subroutine rmarin(ij,kl)
 
167
*     -----------------
 
168
* initializing routine for ranmar, must be called before generating
 
169
* any pseudorandom numbers with ranmar. the input values should be in
 
170
* the ranges 0<=ij<=31328 ; 0<=kl<=30081
 
171
      implicit real*8(a-h,o-z)
 
172
      common/ raset1 / ranu(97),ranc,rancd,rancm
 
173
      common/ raset2 / iranmr,jranmr
 
174
      save /raset1/,/raset2/
 
175
* this shows correspondence between the simplified input seeds ij, kl
 
176
* and the original marsaglia-zaman seeds i,j,k,l.
 
177
* to get the standard values in the marsaglia-zaman paper (i=12,j=34
 
178
* k=56,l=78) put ij=1802, kl=9373
 
179
      write(*,*) "Ranmar initialization seeds",ij,kl
 
180
      i = mod( ij/177 , 177 ) + 2
 
181
      j = mod( ij     , 177 ) + 2
 
182
      k = mod( kl/169 , 178 ) + 1
 
183
      l = mod( kl     , 169 )
 
184
      do 300 ii = 1 , 97
 
185
        s =  0d0
 
186
        t = .5d0
 
187
        do 200 jj = 1 , 24
 
188
          m = mod( mod(i*j,179)*k , 179 )
 
189
          i = j
 
190
          j = k
 
191
          k = m
 
192
          l = mod( 53*l+1 , 169 )
 
193
          if(mod(l*m,64) .ge. 32) s = s + t
 
194
          t = .5d0*t
 
195
  200   continue
 
196
        ranu(ii) = s
 
197
  300 continue
 
198
      ranc  =   362436d0 / 16777216d0
 
199
      rancd =  7654321d0 / 16777216d0
 
200
      rancm = 16777213d0 / 16777216d0
 
201
      iranmr = 97
 
202
      jranmr = 33
 
203
      end