← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

Viewing changes to Template/NLO/Source/ranmar.f

  • Committer: Marco Zaro
  • Date: 2012-08-28 21:06:34 UTC
  • mto: (78.35.14 AutoMint)
  • mto: This revision was merged to the branch mainline in revision 249.
  • Revision ID: marco.zaro@gmail.com-20120828210634-5a06yvda3hplw8ur
doing some renaming:
 Template/FKS-born -> Template/NLO
 fks_born.py -> fks_base.py
 fks_born_helas_objects.py -> fks_helas_objects.py
 export_fks_born.py -> export_fks.py

also functions/classes and tests renamed
all unit tests ok, exporting ok

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/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, joffset
 
16
      integer     ij, kl, iseed1,iseed2
 
17
 
 
18
c
 
19
c     Global
 
20
c
 
21
c-------
 
22
c     18/6/2012 tjs promoted to integer*8 to avoid overflow for iseed > 60K
 
23
c------
 
24
      integer*8       iseed
 
25
      common /to_seed/iseed
 
26
c
 
27
c     Data
 
28
c
 
29
      data init /1/
 
30
      save ij, kl
 
31
c-----
 
32
c  Begin Code
 
33
c-----
 
34
      if (init .eq. 1) then
 
35
         init = 0
 
36
         call get_offset(ioffset)
 
37
         if (iseed .eq. 0) call get_base(iseed)
 
38
c
 
39
c     TJS 3/13/2008
 
40
c     Modified to allow for more sequences 
 
41
c     iseed can be between 0 and 30081*30081
 
42
c     before pattern repeats
 
43
c
 
44
c
 
45
c     TJS 12/3/2010
 
46
c     multipied iseed to give larger values more likely to make change
 
47
c     get offset for multiple runs of single process
 
48
c
 
49
c
 
50
c     TJS 18/6/2012
 
51
c     Updated to better divide iseed among ij and kl seeds
 
52
c     Note it may still be possible to get identical ij,kl for
 
53
c     different iseed, if have exactly compensating joffset, ioffset, jconfig
 
54
c
 
55
         call get_moffset(joffset)
 
56
         joffset = joffset * 3157
 
57
         iseed = iseed * 31300       
 
58
         ij=1802+jconfig + mod(iseed,30081)
 
59
         kl=9373+(iseed/30081)+ioffset + joffset     !Switched to 30081  20/6/12 to avoid dupes in range 30082-31328
 
60
         write(*,'(a,i6,a3,i6)') 'Using random seed offsets',jconfig," : ",ioffset
 
61
         write(*,*) ' with seed', iseed/31300
 
62
         do while (ij .gt. 31328)
 
63
            ij = ij - 31328
 
64
         enddo
 
65
         do while (kl .gt. 30081)
 
66
            kl = kl - 30081
 
67
         enddo
 
68
        call rmarin(ij,kl)         
 
69
      endif
 
70
      call ranmar(x)
 
71
      do while (x .lt. 1d-16)
 
72
         call ranmar(x)
 
73
      enddo
 
74
      x = a+x*(b-a)
 
75
      end
 
76
 
 
77
      subroutine get_base(iseed)
 
78
c-------------------------------------------------------
 
79
c     Looks for file iproc.dat to offset random number gen
 
80
c------------------------------------------------------
 
81
      implicit none
 
82
c
 
83
c     Constants
 
84
c
 
85
      integer    lun
 
86
      parameter (lun=22)
 
87
c
 
88
c     Arguments
 
89
c
 
90
      integer*8 iseed
 
91
c
 
92
c     Local
 
93
c
 
94
      character*60 fname
 
95
      logical done
 
96
      integer i,level
 
97
c-----
 
98
c  Begin Code
 
99
c-----
 
100
 
 
101
      fname = 'randinit'
 
102
      done = .false.
 
103
      level = 1
 
104
      do while(.not. done .and. level .lt. 5)
 
105
         open(unit=lun,file=fname,status='old',err=15)
 
106
         done = .true.
 
107
 15      level = level+1
 
108
         fname = '../' // fname
 
109
         i=index(fname,' ')
 
110
         if (i .gt. 0) fname=fname(1:i-1)
 
111
      enddo
 
112
      if (done) then
 
113
         read(lun,'(a)',end=24,err=24) fname
 
114
         i = index(fname,'=')
 
115
         if (i .gt. 0) fname=fname(i+1:)
 
116
         read(fname,*,err=26,end=26) iseed
 
117
 24      close(lun)
 
118
c         write(*,*) 'Read iseed from randinit ',iseed
 
119
         return
 
120
 26      close(lun)
 
121
      endif
 
122
 25   iseed = 0
 
123
c      write(*,*) 'No base found using iseed=0'
 
124
      end
 
125
 
 
126
      subroutine get_offset(iseed)
 
127
c-------------------------------------------------------
 
128
c     Looks for file iproc.dat to offset random number gen
 
129
c------------------------------------------------------
 
130
      implicit none
 
131
c
 
132
c     Constants
 
133
c
 
134
      integer    lun
 
135
      parameter (lun=22)
 
136
c
 
137
c     Arguments
 
138
c
 
139
      integer iseed
 
140
c
 
141
c     Local
 
142
c
 
143
c-----
 
144
c  Begin Code
 
145
c-----
 
146
 
 
147
      open(unit=lun,file='./iproc.dat',status='old',err=15)
 
148
         read(lun,*,err=14) iseed
 
149
         close(lun)
 
150
         return
 
151
 14   close(lun)
 
152
 15   open(unit=lun,file='../iproc.dat',status='old',err=25)
 
153
         read(lun,*,err=24) iseed
 
154
         close(lun)
 
155
         return
 
156
 24   close(lun)
 
157
 25   iseed = 0
 
158
      end
 
159
 
 
160
      subroutine get_moffset(iseed)
 
161
c-------------------------------------------------------
 
162
c     Looks for file moffset.dat to offset random number gen
 
163
c------------------------------------------------------
 
164
      implicit none
 
165
c
 
166
c     Constants
 
167
c
 
168
      integer    lun
 
169
      parameter (lun=22)
 
170
c
 
171
c     Arguments
 
172
c
 
173
      integer iseed
 
174
c
 
175
c     Local
 
176
c
 
177
c-----
 
178
c  Begin Code
 
179
c-----
 
180
 
 
181
      open(unit=lun,file='./moffset.dat',status='old',err=25)
 
182
         read(lun,*,err=14) iseed
 
183
         write(*,*) "Got moffset",iseed
 
184
         close(lun)
 
185
         return
 
186
 14   close(lun)
 
187
 25   iseed = 0
 
188
      end
 
189
 
 
190
      subroutine ranmar(rvec)
 
191
*     -----------------
 
192
* universal random number generator proposed by marsaglia and zaman
 
193
* in report fsu-scri-87-50
 
194
* in this version rvec is a double precision variable.
 
195
      implicit real*8(a-h,o-z)
 
196
      common/ raset1 / ranu(97),ranc,rancd,rancm
 
197
      common/ raset2 / iranmr,jranmr
 
198
      save /raset1/,/raset2/
 
199
      uni = ranu(iranmr) - ranu(jranmr)
 
200
      if(uni .lt. 0d0) uni = uni + 1d0
 
201
      ranu(iranmr) = uni
 
202
      iranmr = iranmr - 1
 
203
      jranmr = jranmr - 1
 
204
      if(iranmr .eq. 0) iranmr = 97
 
205
      if(jranmr .eq. 0) jranmr = 97
 
206
      ranc = ranc - rancd
 
207
      if(ranc .lt. 0d0) ranc = ranc + rancm
 
208
      uni = uni - ranc
 
209
      if(uni .lt. 0d0) uni = uni + 1d0
 
210
      rvec = uni
 
211
      end
 
212
 
 
213
      subroutine rmarin(ij,kl)
 
214
*     -----------------
 
215
* initializing routine for ranmar, must be called before generating
 
216
* any pseudorandom numbers with ranmar. the input values should be in
 
217
* the ranges 0<=ij<=31328 ; 0<=kl<=30081
 
218
      implicit real*8(a-h,o-z)
 
219
      character*30 filename
 
220
      logical file_exists
 
221
      common/ raset1 / ranu(97),ranc,rancd,rancm
 
222
      common/ raset2 / iranmr,jranmr
 
223
      save /raset1/,/raset2/
 
224
* this shows correspondence between the simplified input seeds ij, kl
 
225
* and the original marsaglia-zaman seeds i,j,k,l.
 
226
* to get the standard values in the marsaglia-zaman paper (i=12,j=34
 
227
* k=56,l=78) put ij=1802, kl=9373
 
228
      write(*,*) "Ranmar initialization seeds",ij,kl
 
229
c
 
230
c    18/6/2012 TJS  Added check to ensure ij and kl are in range
 
231
c      
 
232
      if (ij .lt. 0 .or. ij .gt. 31328 .or.
 
233
     $     kl .lt. 0 .or. kl .gt. 30081) then
 
234
         filename='../../error'
 
235
         INQUIRE(FILE="../../RunWeb", EXIST=file_exists)
 
236
         if(.not.file_exists) filename = '../' // filename
 
237
         open(unit=26,file=filename,status='unknown')
 
238
         if (ij .lt. 0 .or. ij .gt. 31328) then
 
239
            write(26,*) 'Bad initialization value of ij in rmarin ', ij
 
240
            write(*,*) 'Bad initialization value of ij in rmarin ', ij
 
241
         elseif (kl .lt. 0 .or. kl .gt. 30081) then
 
242
            write(26,*) 'Bad initialization value of kl in rmarin ', kl
 
243
            write(*,*) 'Bad initialization value of kl in rmarin ', kl
 
244
         endif
 
245
         stop
 
246
      endif
 
247
 
 
248
      i = mod( ij/177 , 177 ) + 2
 
249
      j = mod( ij     , 177 ) + 2
 
250
      k = mod( kl/169 , 178 ) + 1
 
251
      l = mod( kl     , 169 )
 
252
      do 300 ii = 1 , 97
 
253
        s =  0d0
 
254
        t = .5d0
 
255
        do 200 jj = 1 , 24
 
256
          m = mod( mod(i*j,179)*k , 179 )
 
257
          i = j
 
258
          j = k
 
259
          k = m
 
260
          l = mod( 53*l+1 , 169 )
 
261
          if(mod(l*m,64) .ge. 32) s = s + t
 
262
          t = .5d0*t
 
263
  200   continue
 
264
        ranu(ii) = s
 
265
  300 continue
 
266
      ranc  =   362436d0 / 16777216d0
 
267
      rancd =  7654321d0 / 16777216d0
 
268
      rancm = 16777213d0 / 16777216d0
 
269
      iranmr = 97
 
270
      jranmr = 33
 
271
      end
 
272
 
 
273
 
 
274
 
 
275
 
 
276
 
 
277
 
 
278
 
 
279
 
 
280
 
 
281
 
 
282
      function ran2()
 
283
c-------------------------------------------------------
 
284
c     Front to ranmar which allows user to easily
 
285
c     choose the seed.
 
286
c------------------------------------------------------
 
287
      implicit none
 
288
c
 
289
c     Arguments
 
290
c
 
291
      real*8 ran2
 
292
c
 
293
c     Local
 
294
c
 
295
      real*8 x
 
296
      integer init,ioffset
 
297
      integer ij,kl,iseed1,iseed2
 
298
 
 
299
c
 
300
c     Global
 
301
c
 
302
      integer*8       iseed
 
303
      common /to_seed/iseed
 
304
 
 
305
      integer           iconfig
 
306
      common/to_configs/iconfig
 
307
c
 
308
c     Data
 
309
c
 
310
      data init /1/
 
311
      save ij, kl
 
312
c-----
 
313
c  Begin Code
 
314
c-----
 
315
      if (init .eq. 1) then
 
316
         init = 0
 
317
         call get_offset(ioffset)
 
318
         if (iseed .eq. 0) call get_base(iseed)
 
319
c
 
320
c     TJS 3/13/2008
 
321
c     Modified to allow for more sequences 
 
322
c     iseed can be between 0 and 31328*30081
 
323
c     before pattern repeats
 
324
c
 
325
         ij=1802+iconfig + mod(iseed,30081)
 
326
         kl=9373+(iseed/31328)+ioffset 
 
327
         write(*,'(a,i6,a3,i6)') 'Using random seed offsets',iconfig,
 
328
     &        " : ",ioffset
 
329
         write(*,*) ' with seed', iseed
 
330
         do while (ij .gt. 31328)
 
331
            ij = ij - 31328
 
332
         enddo
 
333
         do while (kl .gt. 30081)
 
334
            kl = kl - 30081
 
335
         enddo
 
336
        call rmarin(ij,kl)         
 
337
      endif
 
338
      call ranmar(x)
 
339
      do while (x .lt. 1d-16)
 
340
         call ranmar(x)
 
341
      enddo
 
342
      ran2=x
 
343
      end