1
c *************************************************************************
3
c This file contains the routines used by the MC@NLO to integrate the
4
c NLO cross sections, and to generate weighted/unweighted events.
5
c Most of these routines are identical to those of the original
6
c package BASES/SPRING. They have been modified in order to treat
7
c the case of functions with arbitrary sign. The user should be
8
c only interested in the two wrappers
11
c which perform the integration and generate unweighted events respectively.
12
c See the comments at the beginning of the bodies of these routines.
14
c ------------------------------ SIDE EFFECTS -----------------------------
16
c When the user is interested in generating weighted events, he can access
17
c the Vegas weight by inserting in the integrand function the following
19
c COMMON /CBSWGT/BSEWGT
20
c This common block is filled by BASES on event-by-event basis. BSEWGT
21
c is actually the Vegas weight, divided by a normalization factor equal
22
c to the number of Vegas calls per iteration, times the number of iterations
23
c in the integration step; BSEWGT can thus be used to fill histograms,
24
c with no further normalization.
26
c It is possible to generate unweighted events with an integrand function
27
c of arbitrary sign. Suppose the integrand function is FF0; then one has
28
c to define FF=ABS(FF0), and give FF as input to BASES/SPRING. The
29
c information on the sign of FF0 must be passed by the user to BASES/SPRING
30
c by inserting the common block
31
c COMMON /CBSSGN/BSFSGN
32
c in the body of FF, where BSFSGN is such that FF0=BSFSGN*FF (thus,
33
c BSFSGN is either +1 or -1). Tipically, the definition of BSFSGN is
34
c the last line of code in the body of the function.
36
c ------------------------------ EVENT RECORD -----------------------------
38
c After SPRING successfully generates an event, it calls (from RUN_OF_SPRING)
41
c The user must use this routine to finalize his computation: selecting
42
c (on statistical basis) a given flavour/color/kinematics configuration,
43
c if the event corresponds to several configurations, and writing the
44
c event on file. It must be stressed that SPRFIN is called
45
c immediately after SPRING has gone through the input function FF in
46
c the generation of an event: thus, FF must pass all the relevant
47
c quantities to SPRFIN through common blocks.
49
c *************************************************************************
54
subroutine run_bases(ff,prefix,ndim,nwild,ncall,it1,it2,ac1,ac2,
55
# res,sig,resneg,signeg,ctime,itd1,itd2,iseed,
57
c Master routine to integrate a real*8 function using bases. The entries are:
60
c prefix: bases will write the relevant information on prefix_bs.data
61
c ndim: the dimension of the space that contains the support of ff
62
c nwild: the number of wild variables
63
c ncall: the number of sampling points per iteration (user's input)
64
c it1: the maximum number of iterations in the grid optimization step
65
c it2: the maximum number of iterations in the integration step
66
c ac1: the relative accuracy required in the grid optimization step
67
c ac2: the relative accuracy required in the integration step
68
c res: the result of the integration of ff
69
c sig: its standard deviation
70
c resneg: the result of the integration of ff0, where ff=abs(ff0)
71
c signeg: its standard deviation
72
c ctime: the computing time used in the integration step
73
c itd1: iterations performed in the grid optimization step
74
c itd2: iterations performed in the integration step
75
c iseed: the initial seed for random number generation
76
c iwrite: if equal to zero, bases does NOT write prefix_bs.data
77
c ircall: the actual number of sampling points per iteration (bases output)
79
c ******************************** WARNING ********************************
80
c Resneg and signeg return sensible results only if the common block
81
c COMMON /CBSSGN/BSFSGN
82
c is inserted in the body of the integrand function ff
83
c ******************************** WARNING ********************************
84
c Bases has been modified in order to output weighted events, thus
85
c avoiding the use of spring. For this option to work, the common block
86
c COMMON /CBSWGT/BSEWGT
87
c must be intserted in the body of the integrand function ff
89
c In the case of the wild variables, bases saves the probability
90
c information for each hypercube associated to those variables;
91
c this is crucial when spring is attached. If a variable is not
92
c wild, the information on the shape of the integrand is lost
93
c during the phase of event generation.
94
c Bases also leaves to the user the choice of whether a variable
95
c is optimized (i.e., the grid is iteratively adjusted), or not.
96
c This option is controlled by the value of IG(i) (1=optimization,
97
c 0=no optimization). In the present interface, we always set IG(i)=1
99
real*8 ff,ac1,ac2,res,sig,resneg,signeg,ctime
100
integer ndim,nwild,ncall,it1,it2,itd1,itd2,iseed,iwrite,ircall
101
integer ing,inpg,insp
102
common/brcall/ing,inpg,insp
103
character * 60 prefix,fname
105
parameter (isunit=23)
107
data verbose/.false./
110
call init_of_bases(ndim,nwild,ncall,it1,it2,iseed,ac1,ac2)
111
call bases(ff,res,sig,resneg,signeg,ctime,itd1,itd2)
114
call fk88strcat(prefix,'_bs.data',fname)
115
open(unit=isunit,file=fname,
116
# form='unformatted',status='unknown')
121
write(6,*)'integral: ',res,' +- ',sig
122
write(6,*)'elapsed time: ',ctime
123
write(6,*)'itns: ',itd1,itd2
129
subroutine run_spring(ff,prefix,mxevts,mxtrls,nevts,ntrls,
131
c Master routine to generate events using spring. The entries are:
133
c ff: the function previously integrated by bases
134
c prefix: spring will read the relevant information from prefix_bs.data
135
c mxevts: the required number of generated events
136
c mxtrls: the maximum number of trials to generate mxevts
137
c nevts: the number of generated events
138
c ntrls: the number of trials
139
c ndim: the dimension of the space that contains the support of ff
140
c nwild: the number of wild variables
141
c iseed: the initial seed for random number generation
143
c ******************************** WARNING ********************************
144
c For the event to be finalised (for example, written on file),
147
c must be written by the user; the relevant information is passed
148
c from FF to SPRFIN by means of common blocks
151
integer mxevts,mxtrls,nevts,ntrls,ndim,nwild,iseed,idum
152
character * 60 prefix
153
parameter (dum=-1.d0)
157
if(mxevts.eq.0)return
158
call init_of_spring(prefix)
159
call init_of_bases(ndim,nwild,idum,idum,idum,iseed,dum,dum)
161
write(6,*)'no events will be generated'
164
if(mxtrls.lt.mxevts)then
166
# 'SPRING warning: the required number of events cannot be'
168
# 'generated: mxevts, mxtrls=',mxevts,mxtrls
170
call run_of_spring(ff,mxevts,mxtrls,nevts,ntrls)
171
if(nevts.ne.mxevts.and.ntrls.ne.mxtrls)then
172
write(6,*)'SPRING: fatal error'
173
write(6,*)mxevts,nevts
174
write(6,*)mxtrls,ntrls
182
c Begin of BASES routines
185
subroutine init_of_bases(mdim,mwild,mcall,it1,it2,iseed,ac1,ac2)
186
c This subroutine is called in order to initialize bases parameters;
187
c it has been obtained by getting together bsdims.f, bsgrid.f, and
188
c bsparm.f of the original package (BASES50/SPRING50). When the last
189
c four entries are set to negative values, default values (defined
190
c in bsinit) are used
192
integer mdim,mwild,mcall,it1,it2,iseed
195
parameter (mxdim = 50)
196
parameter (mxwl = 15)
197
* NDIM : The number of dimension of integral *
198
* NWILD : The number of wild variables *
199
* NCALL : The number of sample points per iteration. *
200
* This actual number is to be determined by taking the *
201
* number of dimensions into account. *
202
* XL(i) : The lower value of the i-th integral variable *
203
* XU(i) : The upper value of the i-th integral variable *
204
* IG(i) : The flag switches whether the grid of i-th variable *
205
* is to be optimized ( 1 ) or kept uniform ( 0 ). *
206
real * 8 xl(mxdim),xu(mxdim)
207
integer ndim,nwild,ncall,ig(mxdim)
208
common/bparm1/xl,xu,ndim,nwild,ig,ncall
209
* ACC1 : The required accuracy at the grid optimization step *
210
* ACC2 : The required accuracy at the integration step. *
211
* ITMX1 : The max. number of iteration at the grid opt. step. *
212
* ITMX2 : The max. number of iteration at the integration step. *
215
common/bparm2/acc1,acc2,itmx1,itmx2
217
common/base0/jflag,ibases
220
common/bfxlhisto/lhisto
223
if(mdim.gt.mxdim.or.mwild.gt.mxwl)then
224
write(6,*)'just kidding....'
238
if(it1.ge.0)itmx1=it1
239
if(it2.ge.0)itmx2=it2
240
if(ac1.ge.0.d0)acc1=ac1
241
if(ac2.ge.0.d0)acc2=ac2
246
************************************************************************
247
* ==================================================== *
248
SUBROUTINE BASES(FXN,S,SIGMA,SNEG,SIGMANEG,CTIME,IT1,IT2 )
249
* ==================================================== *
250
* Subroutine BASES for the Numerical integration. *
251
* In terms of this program Integration can be done, furthermore *
252
* a probability distribution can be made for the event generation.*
253
* The event with weight one is generated by program SPRING. *
255
* from the arguement *
256
* FXN : Name of function program *
257
* from the labeled common /BASE1/ *
258
* XL(50) : Lower limits of the integration variabels *
259
* XU(50) : upper limits of the integration variabels *
260
* NDIM : Dimension of the integration *
261
* NCALL : Number of sampling points per iteration *
262
* from the lebeled common /BASE2/ *
263
* ITMX* : Number of iteration *
264
* ACC* : Required accuracies *
266
* S : Estimate of the integral *
267
* SIGMA : Standard deviation of the estimate *
268
* SNEG : Estimate of the integral with sign (see run_bases) *
269
* SIGMANEG: Standard deviation of the estimate for sneg *
270
* CTIME : Computing time required for integration *
271
* IT1 : Number of iterations for the grid defining step *
272
* IT2 : Number of iterations for the integration step *
274
C* Coded by S.Kawabata April '94 *
276
C* Modified by S. Frixione (sneg, sigmaneg, and common block brsneg *
277
C* have been added) *
278
C***********************************************************************
281
IMPLICIT REAL*8 (A-H,O-Z)
283
PARAMETER (MXDIM = 50)
285
* JFLAG = 0 : First trial for defining grids.
286
* JFLAG = 1 : First trial for event accumulation.
287
* JFLAG = 2 : Second or more trial for defining grids.
288
* JFLAG = 3 : Second or more trial for accumulation.
290
COMMON /BASE0/ JFLAG,IBASES
291
COMMON /BASE1/ XL(MXDIM),XU(MXDIM),NDIM,NWILD,
293
COMMON /BASE2/ ACC1,ACC2,ITMX1,ITMX2
295
COMMON /BSRSLT/AVGI,SD,CHI2A,STIME,ITG,ITF
296
COMMON /BSRNEG/AVGINEG,SDNEG
298
COMMON /BWARN1/ NERROR
299
COMMON /BWARN2/ ERROR(3,3)
300
* INTV = ( 0 / 1 / any ) = ( Batch / Batch(Unix) / Interactive )
301
* IPNT = ( 0 / any ) = ( IBM Type / Ascii printer )
302
COMMON /BSCNTL/ INTV, IPNT, NLOOP, MLOOP
304
COMMON/NINFO/ NODEID, NUMNOD
305
COMMON /BDATE/ IDATE(3),ITIME(2)
306
* IDATE(1) : year ITIME(1) : hour
307
* IDATE(2) : month ITIME(2) : minute
309
REAL*4 TIMEBS,TIMINT,TIMESP,TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
310
COMMON /BTIME1/ TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
311
COMMON /BTIME2/ TIMEBS(0:2),TIMINT,TIMESP(0:2)
313
*-------------------------------------------------
314
* Check the parameters defined by user
315
*------------------------------------------------------
319
* ---------------------------------------------------------------
321
* ---------------------------------------------------------------
327
IF( INTV .GT. 1 ) THEN
328
CALL BSPRNT( LU, 1, IDUM1, IDUM2 )
331
C -----------------------------------------------------
333
C -----------------------------------------------------
341
IF( INTV .GT. 1 ) THEN
342
CALL BSPRNT( LU, 4, IDUM1, IDUM2 )
347
* ===================
349
* =================== For a parallel computer
350
C CALL BSCAST( JFLAG, 1 )
352
* ----------------------------------------------------
353
* Accumulation to make probability distribution
354
* ----------------------------------------------------
355
* ===================
357
* =================== For a parallel computer
358
C CALL BSCAST( JFLAG, 1 )
370
IF( NERROR .GT. 0 ) THEN
372
9900 FORMAT(1X,'****************************************',
373
. '***************************************',
374
. /1X,'* (((( Warning in the integration step ',
380
WRITE(6,9901) ERROR(I,J)
386
. /1X,'*(( Suggestion )) ',
388
. /1X,'* (1) Try integration again with larger ',
389
. 'number of sample points than this job.*',
392
. /1X,'* (2) The integral variables are not sui',
393
. 'ted for the function. *',
394
. /1X,'* Take another integral variables !!',
398
. /1X,'****************************************',
399
. '***************************************')
402
IF( INTV .GT. 1 ) THEN
403
CALL BSPRNT( LU, 2, IDUM1, IDUM2 )
410
************************************************************************
411
* =================== *
413
* =================== *
415
* To check user's initialization parameters. *
417
* Coded by S.Kawabata Oct. '85 *
419
************************************************************************
421
IMPLICIT REAL*8 (A-H,O-Z)
422
PARAMETER ( MXDIM = 50)
423
COMMON /BPARM1/ XL(MXDIM),XU(MXDIM),NDIM,NWILD,
425
COMMON /BPARM2/ ACC1,ACC2,ITMX1,ITMX2
427
COMMON /BASE0/ JFLAG,IBASES
428
COMMON /BASE1/ XLT(MXDIM),XUT(MXDIM),NDIMT,NWILDT,
430
COMMON /BASE2/ ACC1T,ACC2T,ITMX1T,ITMX2T
431
COMMON /BSCNTL/ INTV, IPNT, NLOOP, MLOOP
436
IF( IBASES .NE. 1 ) THEN
439
. 5X,'*************************************************',
441
. /5X,'* BSINIT was not called before calling BASES *',
443
. /5X,'* Process was terminated due to this error. *',
445
. /5X,'*************************************************')
450
IF( NDIM .LT. 1) THEN
453
. 5X,'*************************************************',
455
. /5X,'* NDIM was not set before calling BASES. *',
457
. /5X,'* Process was terminated due to this error. *',
459
. /5X,'*************************************************')
466
IF( XU(I) .LE. -1.0D37) THEN
469
. 5X,'*************************************************',
471
. /5X,'* XL(',I6,' ). XU(',I6,' ) were not set *',
472
. /5X,'* before calling BASES. *',
473
. /5X,'* Process was terminated due to this error. *',
475
. /5X,'*************************************************')
485
C Change the maximum number of the wild variables
487
IF( NWILD .LT. 0) THEN
488
NWILD = MIN( NDIM, 15)
491
. 5X,'*************************************************',
493
. /5X,'* NWILD was not set before calling BASES. *',
495
. /5X,'* NWILD is set equal to the value(',I6,' ). *',
497
. /5X,'*************************************************')
499
.IF( NWILD .GT. 15) THEN
501
NWILD = MIN( NDIM, 15)
502
WRITE(6,9400) NWILDO, NWILD
504
. 5X,'*************************************************',
506
. /5X,'* NWILD(',I6,' ) was too large number. *',
508
. /5X,'* NWILD is set equal to the value(',I6,' ). *',
510
. /5X,'*************************************************')
525
C***********************************************************************
527
C*======================== *
528
C* SUBROUTINE BSETGU *
529
C*======================== *
531
C* Initialization of Bases progam *
532
C* This is called only when IFLAG=0. *
533
C* ( IFLAG = 0 ; First Trial of Defining Grid step ) *
535
C* Changed by S.Kawabata Aug. 1984 at Nagoya Univ. *
536
C* Last update Oct. 1985 at KEK *
538
C* Modified by S. Frixione (common block brcall has been added in *
539
C* order to pass npg, ng, and nsp without using base4) *
540
C***********************************************************************
544
IMPLICIT REAL*8 (A-H,O-Z)
545
PARAMETER (MXDIM = 50, NDMX = 50, LENG = 32768)
546
COMMON /BASE1/ XL(MXDIM),XU(MXDIM),NDIM,NWILD,
548
COMMON /BASE4/ XI(NDMX,MXDIM),DX(MXDIM),DXD(LENG),DXP(LENG),
549
. ND,NG,NPG,MA(MXDIM)
550
COMMON /BASE6/ D(NDMX,MXDIM),
551
. ALPH,XSAVE(NDMX,MXDIM),XTI,XTSI,XACC,ITSX
552
COMMON /BRCALL/ING,INPG,INSP
557
C---------------------------------------------------------------
558
C Define the number of grids and sub-regions
559
C---------------------------------------------------------------
560
C==> Determine NG : Number of grids
561
NG = (NCALL/2.)**(1./NWILD)
562
IF(NG .GT. 25) NG = 25
563
100 IF(NG .LT. 2) NG = 1
564
IF(NG**NWILD .GT. LENG) THEN
569
C==> Determine ND : Number of sub-regions
573
C==> Determine NPG: Number of sampling points per subhypercube
585
IF( NDIM .GT. 1 ) THEN
589
IF( J .LE. NWILD ) THEN
595
C---------------------------------------------------------------
596
C Set size of subregions uniform
597
C---------------------------------------------------------------
609
145 IF(RC .GT. DR) GO TO 140
612
XIN(I)= XN-(XN-XO)*DR
613
IF(I .LT. NDM) GO TO 145
619
********************************************* Updated Feb.08 '94
620
IF( ITSX .GT. 0 ) THEN
636
C***********************************************************************
638
C*======================== *
639
C* SUBROUTINE BSETGV( IFLAG ) *
640
C*======================== *
642
C* Refine the grid sizes *
644
C* Coded by S.Kawabata Aug. 1984 at Nagoya Univ. *
645
C* Last update Oct. 1985 at KEK *
647
C***********************************************************************
649
SUBROUTINE BSETGV( IFLAG )
651
IMPLICIT REAL*8 (A-H,O-Z)
652
PARAMETER (MXDIM = 50, NDMX = 50, LENG = 32768)
653
COMMON /BASE1/ XL(MXDIM),XU(MXDIM),NDIM,NWILD,
655
COMMON /BASE4/ XI(NDMX,MXDIM),DX(MXDIM),DXD(LENG),DXP(LENG),
656
. ND,NG,NPG,MA(MXDIM)
657
COMMON /BASE3/ SCALLS,WGT,TI,TSI,TACC,IT
658
COMMON /BASE6/ D(NDMX,MXDIM),
659
. ALPH,XSAVE(NDMX,MXDIM),XTI,XTSI,XACC,ITSX
661
COMMON /BSRSLT/AVGI,SD,CHI2A,STIME,ITG,ITF
664
DIMENSION XIN(NDMX),R(NDMX),DT(MXDIM),DDX(NDMX)
665
DATA ONE/1.0D0/,ZERO/0.0D0/,N0/0/,N1/1/
667
*========= Save the grid information for the best accuracy ===========
669
IF( ITSX .GT. 0 ) THEN
670
IF( IFLAG .EQ. 0 ) THEN
672
IF( ( TI .GT. AVGI+SD) .AND. TSI .LT. XTSI ) THEN
684
IF( ( XTI .GT. TI) .AND. XTSI .LT. TSI ) THEN
696
C======= SMOOTHING THE FUNCTION D(I,J)
698
CLOGE = 1.0D0/LOG(10.0D0)
702
IF( IG(J) .EQ. 1 ) THEN
703
DDX(1)= 0.5D0*(D(1,J) + D(2,J))
705
DDX(I)= (D(I+1,J) + D(I,J) + D(I-1,J))/3.D0
707
DDX(ND) = 0.5D0*(D(NDM,J) + D(ND,J))
714
C=========== REDEFINE THE GRID
722
IF(D(I,J) .GT. ZERO) THEN
724
IF( DT10 - CLOGE*DILOG .LE. 70.0D0 ) THEN
726
R(I) = ((XO-ONE)/(XO*(DTLOG-DILOG)))**ALPH
729
R(I) = (DTLOG-DILOG)**(-ALPH)
743
750 IF(RC.GT.DR)GO TO 740
746
XIN(I)= XN-(XN-XO)*DR/R(K)
747
IF(I.LT.NDM)GO TO 750
759
C***********************************************************************
761
C*======================== *
762
C* SUBROUTINE BSGETW( WEIGHT ) *
763
C*======================== *
767
C* Coded by T.Ishikawa Jun. 1995 at KEK *
768
C* Last update Jun. 1995 at KEK *
770
C***********************************************************************
772
SUBROUTINE BSGETW( WEIGHT )
774
IMPLICIT REAL*8 (A-H,O-Z)
775
COMMON /BASE3/ SCALLS,WGT,TI,TSI,TACC,IT
777
*========= Save the grid information for the best accuracy ===========
785
***********************************************************************
786
*============================ *
787
SUBROUTINE BSINFO( LU )
788
*============================ *
790
* Print the information for *
791
* (1) BASES parameters *
792
* (2) Computer time information *
793
* (3) Convergency behavior of the Grid optimization step *
794
* (4) Convergency behavior of the integration step *
796
* LU : Logical unit number of printer *
798
* by S.Kawabata March 1994 at KEK
800
***********************************************************************
802
IMPLICIT REAL*8 (A-H,O-Z)
804
COMMON /BSRSLT/AVGI,SD,CHI2A,STIME,ITG,ITF
808
CALL BSPRNT( LU, 1, IDUM1, IDUM2 )
810
* Print Bases parameters
812
CALL BSPRNT( LU, 4, IDUM1, IDUM2 )
814
* Print Computing time information
816
CALL BSPRNT( LU, 3, IDUM1, IDUM2 )
818
* Print Convergency Behaviors
822
IF( ISTEP .EQ. 1 ) ITX = ITF
824
IF( ITX .GT. 0 ) THEN
826
CALL BSPRNT( LU, 8, ITX, ISTEP )
835
************************************************************************
836
* =================== *
838
* =================== *
840
* Initialization of BASE50/SPRING50. *
841
* Function of this routine is *
842
* (0) Set the size of histogram and scatter plot buffers *
843
* (1) Set the parameters INTV and IPNT *
844
* INTV = ( 0 / 1 / any ) *
845
* = ( Batch / Batch(Unix) / Interactive ) *
846
* IPNT = ( 0 / any ) *
847
* = ( IBM Type / Ascii printer ) *
848
* (2) Set the acceleration factor ALPHA by 1.5 *
849
* The range of this value is from 0.0 to 2.0. *
850
* ALPHA = 0.0 results in no grid-optimization. *
851
* (3) Set the grid-optimization flag IGOPT ( Default value 0 ) *
852
* IGOPT = 0 : The grid is optimized by VEGAS algorithm *
853
* IGOPT = 1 : The grid is optimized so that the accuracy *
854
* of each iteration be minimized. *
855
* (4) Set Node-ID number NODEID and the number of nodes NUMNOD *
856
* (5) Set seed of radom number *
857
* (6) Set the values of BASES paremeters with default ones. *
858
* (7) Set the values of parameters with non-sense values, *
859
* which should be set again with the true values by User *
860
* before running BASES. *
862
* Coded by S.Kawabata March '94 *
864
* Modified by S. Frixione in order to exclude the built-in *
865
* histogramming package. The common block bfxlhisto has been added *
866
************************************************************************
868
IMPLICIT REAL*8 (A-H,O-Z)
869
PARAMETER (MXDIM = 50, NDMX = 50 )
870
COMMON /BPARM1/ XL(MXDIM),XU(MXDIM),NDIM,NWILD,
872
COMMON /BPARM2/ ACC1,ACC2,ITMX1,ITMX2
874
COMMON /BASE0/ JFLAG,IBASES
875
COMMON /BASE6/ D(NDMX,MXDIM),
876
. ALPH,XSAVE(NDMX,MXDIM),XTI,XTSI,XACC,IGOPT
877
COMMON /BSCNTL/ INTV, IPNT, NLOOP, MLOOP
878
COMMON/NINFO/ NODEID, NUMNOD
879
COMMON /BDATE/ IDATE(3),ITIME(2)
880
* IDATE(1) : year ITIME(1) : hour
881
* IDATE(2) : month ITIME(2) : minute
883
REAL*4 TIMEBS,TIMINT,TIMESP,TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
884
COMMON /BTIME1/ TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
885
COMMON /BTIME2/ TIMEBS(0:2),TIMINT,TIMESP(0:2)
887
COMMON/BFXLHISTO/LHISTO
888
*=========================================================
889
* (0) Initialization of timer and Histogram buffer
890
* Timer initialization
891
CALL BSTIME( TIME0, 0 )
895
* Histogram buffer initialization
901
*=========================================================
903
* (1) Set the parameters INTV and IPNT
906
* (2) Set the acceleration factor ALPHA by 1.5
908
* (3) Set the grid-optimization flag IGOPT
910
* (4) Set Node-ID number NODEID and the number of nodes NUMNOD
911
* IF( INTV .EQ. 0 ) THEN
919
C---------------------------------------------------------------
920
C (5) Set initial seeds of random number generator
921
c Modified by SF: this is now in init_of_bases
922
C---------------------------------------------------------------
925
c CALL DRNSET( ISEED )
926
C ---------------------------------------------------------------
927
C (6),(7) Set BASES parameters equal to default values
928
C ---------------------------------------------------------------
942
* Initialization of computing time table of BASES
947
*-------------------------------------------
948
* Don't change IBASES from this value
949
*-------------------------------------------
956
***********************************************************************
958
* ========================== *
959
SUBROUTINE BSINTG( FXN )
960
* ========================== *
962
* Subroutine performs N-dimensional Monte Carlo integration *
963
* for four vector generation of simulated events *
965
* JFLAG = 0 ; First Trial of Defining Grid *
966
* JFLAG = 1 ; First Trial of Data Accumulation *
967
* JFLAG = 2 ; Second Trial of Defining Grid *
968
* JFLAG = 3 ; Second Trial of Data Accumulation *
970
* Coded by S.Kawabata July 1980 at DESY, Hamburg *
971
* Last update March 1994 *
973
* Modified by S. Frixione in order to exclude the built-in *
974
* histogramming package. The common block bfxlhisto has been added *
975
* The common block bsrneg and CBSSGN have been added to treat *
976
* functions with sign; CBSSGN must mandatorily be inserted also in *
977
* in the integrand function. The commond block CBSWGT has been *
978
* added in order to deal with weighted events; it must be inserted *
979
* also in the integrand function. *
980
***********************************************************************
982
IMPLICIT REAL*8 (A-H,O-Z)
985
PARAMETER (MXDIM = 50, NDMX = 50, LENG = 32768)
986
COMMON /BASE0/ JFLAG,IBASES
987
COMMON /BASE1/ XL(MXDIM),XU(MXDIM),NDIM,NWILD,
989
COMMON /BASE2/ ACC1,ACC2,ITMX1,ITMX2
990
COMMON /BASE3/ SCALLS,WGT,TI,TSI,TACC,IT
991
COMMON /BASE4/ XI(NDMX,MXDIM),DX(MXDIM),DXD(LENG),DXP(LENG),
992
. ND,NG,NPG,MA(MXDIM)
994
REAL*4 TIME, EFF, WRONG, TRSLT, TSTD, PCNT
995
COMMON /BASE5/ ITRAT(ITM,0:1),TIME(ITM,0:2),EFF(ITM,0:1),
996
. WRONG(ITM,0:1),RESLT(ITM,0:1),ACSTD(ITM,0:1),
997
. TRSLT(ITM,0:1),TSTD(ITM,0:1),PCNT(ITM,0:1)
998
COMMON /BASE6/ D(NDMX,MXDIM),
999
. ALPH,XSAVE(NDMX,MXDIM),XTI,XTSI,XACC,ITSX
1001
COMMON /BSRSLT/AVGI,SD,CHI2A,STIME,ITG,ITF
1002
COMMON /BSRNEG/AVGINEG,SDNEG
1003
COMMON /BRCALL/ING,INPG,INSP
1004
COMMON /CBSSGN/BSFSGN
1005
COMMON /CBSWGT/BSEWGT
1007
COMMON /BWARN1/ NERROR
1008
COMMON /BWARN2/ ERROR(3,3)
1010
* INTV = ( 0 / 1 / any ) = ( Batch / Batch(Unix) / Interactive )
1011
* IPNT = ( 0 / any ) = ( IBM Type / Ascii printer )
1012
COMMON /BSCNTL/ INTV, IPNT, NLOOP, MLOOP
1015
INTEGER KG(MXDIM),IA(MXDIM)
1017
COMMON/NINFO/ NODEID, NUMNOD
1018
REAL*4 TIMEBS,TIMINT,TIMESP,TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
1019
COMMON /BTIME1/ TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
1020
COMMON /BTIME2/ TIMEBS(0:2),TIMINT,TIMESP(0:2)
1022
INTEGER NCNODE(2,512),NPNODE(2,512)
1025
COMMON/BFXLHISTO/LHISTO
1027
* Parameters for checking convergency
1029
DATA ACLMT,FC / 25.0D0, 5.0D0 /
1032
DATA ONE/ 1.0D0/, ZERO/0.0D0/, LU / 6/
1033
DATA N0/0/, N1/1/, HUNDRT/100.0D0/
1035
************************************************************************
1036
* Initialization Part
1037
************************************************************************
1038
*=======================================================================
1039
* Determine the number of hypercubes NSP
1040
*=======================================================================
1050
DV2G = DXG**(2*NWILD)/NPG/NPG/(NPG-1)
1053
IF( NSP .EQ. 1 ) THEN
1054
*=======================================================================
1055
* Determination of the number of sampling points
1056
* per node in the single hypercube case
1057
*=======================================================================
1058
MEX = MOD(NPG,NUMNOD)
1061
DO 12 NODEX = 1,NUMNOD
1063
NPGT = NPGT + NPERCP
1064
IF( NODEX .LE. MEX ) NPGT = NPGT + 1
1067
NPNODE(1,NODEX) = NPGS
1068
NPNODE(2,NODEX) = NPGT
1071
*=======================================================================
1072
* Determination of the number of hypercubes
1073
* per node in many hypercubes case
1074
*=======================================================================
1075
MEX = MOD(NSP,NUMNOD)
1078
DO 15 NODEX = 1,NUMNOD
1080
NSPT = NSPT + NPERCP
1081
IF( NODEX .LE. MEX ) NSPT = NSPT + 1
1082
NCNODE(1,NODEX) = NSPS
1083
NCNODE(2,NODEX) = NSPT
1085
NPNODE(2,NODEX) = NPG
1088
*=======================================================================
1098
IF(JFLAG .EQ. N0 .OR. JFLAG .EQ. N1 ) THEN
1099
*-----------------------------------------------------------------------
1100
* JFLAG = 0 : The first trial of the grid optim. step
1101
* JFLAG = 1 : The first trial of the integration step
1102
*-----------------------------------------------------------------------
1127
*-----------------------------------------------------------------------
1128
* JFLAG = 2 : The continuation of the grid optim. step
1129
* JFLAG = 3 : The continuation of the integration step
1130
*-----------------------------------------------------------------------
1131
C IF( JFLAG .EQ. 2 ) THEN
1136
C . IF( JFLAG .EQ. 3 ) THEN
1149
*------- Set the expected accuracy and the max. iteration number -------
1153
IF( ISTEP .EQ. N1 ) THEN
1158
*-------- Print the title of the convergency behavior table -----------
1159
* in the interactive mode
1160
IF( INTV .GT. 1 ) THEN
1161
* -----------------------------------
1162
CALL BSPRNT( LU, 5, ISTEP, IDUM2 )
1163
* -----------------------------------
1167
* =====================
1169
* =====================
1171
*********************************************************************
1172
* Main Integration Loop
1173
*********************************************************************
1175
DO 500 IT = IT1,ITMX
1177
*=======================================================================
1178
* Initialization for the iteration
1179
*=======================================================================
1181
SCALLS = SCALLS + CALLS
1189
IF( ISTEP .EQ. N0 ) THEN
1197
IF( NODEID .EQ. 0 ) NODEX = NUMNOD
1199
*---------------------------------------------------------------------
1200
* Distributing hyper cubes to NumNode nodes
1201
* NCNODE(1,NODEX) : 1st cube number for the node NODEX
1202
* NCNODE(2,NODEX) : Last cube number for the node NODEX
1203
* NODEX : node number 1 => NumNode(=0)
1204
* NODEX : node number 1 => NumNode(=0)
1205
*---------------------------------------------------------------------
1207
NSP1 = NCNODE(1,NODEX)
1208
NSP2 = NCNODE(2,NODEX)
1209
* Dummy loopfor a parallel processor
1210
C IF( NSP1 .GT. 1 ) THEN
1211
C CALL DRLOOP( NDIM*NPG*(NSP1-1) )
1214
*=====================================================================
1215
* Loop for hypercube from NSP1 to NSP2 in the NodeX-th node
1216
*=====================================================================
1218
DO 400 NCB = NSP1, NSP2
1225
IF( NWILD .GT. 1 ) THEN
1226
DO 210 J = 1,NWILD-1
1227
NUM = MOD(NP,MA(J+1))
1228
KG(J) = NUM/MA(J) + 1
1231
KG(NWILD) = NP/MA(NWILD) + 1
1233
*---------------------------------------------------------------------
1234
* If number of hypercubes is only one,
1235
* Distributing sampling points to NumNode nodes
1236
* NPNODE(1,NODEX) : 1st sample point for the node NODEX
1237
* NPNODE(2,NODEX) : Last sample point for the node NODEX
1238
* NODEX : node number 1 => NumNode(=0)
1239
*---------------------------------------------------------------------
1241
NPG1 = NPNODE(1,NODEX)
1242
NPG2 = NPNODE(2,NODEX)
1243
* Dummy loop for a parallel processor
1244
C IF( NPG1 .GT. 1 ) THEN
1245
C CALL DRLOOP( NDIM*(NPG1-1) )
1248
*=====================================================================
1249
* Loop for sampling points from NPG1 to NPG2
1250
* in the single hypercube case
1251
*=====================================================================
1253
DO 300 NTY = NPG1,NPG2
1255
*---------------------------------------------------------------------
1256
* Determine the integration variables by random numbers
1257
*---------------------------------------------------------------------
1261
IF( J .LE. NWILD ) THEN
1262
XN = (KG(J)-DRN(IDUMY))*DXG+1.D0
1264
XN = ND*DRN(IDUMY)+1.D0
1268
IF( IAJ .EQ. 1) THEN
1272
XO = XI(IAJ,J)-XI(IAJ-1,J)
1273
RC = XI(IAJ-1,J)+(XN-IAJ)*XO
1275
X(J) = XL(J)+RC*DX(J)
1278
*-----------------------------------------------------------------------
1280
c The Vegas weight passed to the integrand function is zero in
1281
c the grid optimisation step
1283
IF(JFLAG.EQ.1)BSEWGT = WGT/(ITMX2*INSP*INPG)
1287
*-----------------------------------------------------------------------
1288
* Check the value of the integrand
1289
*-----------------------------------------------------------------------
1291
IF( FXG .NE. 0.0 ) THEN
1293
IF( ISTEP .EQ. 1 ) THEN
1294
DXD(NCB) = DXD(NCB) + FXG
1295
IF( FXG .GT. DXP(NCB) ) DXP(NCB) = FXG
1297
IF( FXG .LT. 0.0 ) THEN
1299
IF( NEGFLG .EQ. 0 ) THEN
1300
WRITE(6,9200) IT,NODEID
1302
. '******* WARNING FROM BASES ********',
1304
. /1X,'* Negative FUNCTION at IT =',I3,1X,
1305
. ', node = ',I3,1X,'*',
1306
. /1X,'***********************************',
1313
*-----------------------------------------------------------------------
1314
* Accumulation of FXG and FXG*FXG
1315
*-----------------------------------------------------------------------
1320
F2NEG = FXGNEG*FXGNEG
1321
FBNEG = FBNEG + FXGNEG
1322
F2BNEG = F2BNEG + F2NEG
1323
IF( ISTEP .EQ. 0 ) THEN
1325
D(IA(J),J)= D(IA(J),J)+F2
1331
*------------------------------------------- for a parallel processor
1332
* Dummy loop for a parallel processor
1333
C IF( NPG2 .LT. NPG ) THEN
1334
C CALL DRLOOP(NDIM*(NPG-NPG1))
1336
* Global sum of FB and F2B
1337
C IF( NSP .EQ. 1 ) THEN
1338
C CALL BSDSUM( FB, 1 )
1339
C CALL BSDSUM( F2B, 1 )
1341
*-----------------------------------------------------------------------
1343
*-----------------------------------------------------------------------
1344
* Calculate the estimate and variance in the hypercube
1345
*-----------------------------------------------------------------------
1347
F2B = DSQRT(F2B*NPG)
1348
F2S = (F2B-FB)*(F2B+FB)
1351
F2BNEG = DSQRT(F2BNEG*NPG)
1352
F2SNEG = (F2BNEG-FBNEG)*(F2BNEG+FBNEG)
1354
TSINEG = TSINEG + F2SNEG
1358
*------------------------------------------- for a parallel processor
1360
C IF( NSP2 .LT. NSP ) THEN
1361
C CALL DRLOOP(NDIM*NPG*(NSP-NSP2))
1364
* Global sum of efficiency and frequency
1365
* of negative valued function
1368
C CALL BSISUM( NEFF, 2 )
1372
C IF( NSP .EQ. 1 ) THEN
1373
C CALL BSDSUM( TX, 2 )
1376
* Global sum of grid information
1377
C IF( ISTEP .EQ. 0 ) THEN
1378
C NOWORK = NDMX*NDIM
1379
C CALL BSDSUM( D, NOWORK )
1382
*=====================================================================
1383
* Compute Result of this Iteration
1384
*=====================================================================
1385
*--------------------------------------------------------------------
1386
* Accumulate the histogram entries
1387
*--------------------------------------------------------------------
1393
*--------------------------------------------------------------------
1403
TSINEG = TSINEG*DV2G
1405
IF( TSI .LE. 1.0D-37 ) TSI = 1.0D-37
1406
IF( TSINEG .LE. 1.0D-37 ) TSINEG = 1.0D-37
1409
TI2NEG = TINEG*TINEG
1411
IF( NGOOD .LE. 10 ) THEN
1412
* --------------------------------
1413
CALL BSPRNT( LU, 9, IDUM1, IDUM2 )
1414
* --------------------------------
1421
*--------------------------------------------------------------------
1422
* Calculate the cumulative result
1423
*--------------------------------------------------------------------
1430
IF(IT .GT. N1 ) CHI2A = (SCHI - SI*AVGI)/(IT-.999D0)
1431
SD = DSQRT(ONE/SWGT)
1433
SINEG = SINEG+TINEG*WGTNEG
1434
SWGTNEG = SWGTNEG+WGTNEG
1435
SCHINEG = SCHINEG+TI2NEG*WGTNEG
1436
AVGINEG = SINEG/SWGTNEG
1437
SDNEG = DSQRT(ONE/SWGTNEG)
1439
*---------------------------------------------------------------------
1440
* Save the results in the buffer
1441
*---------------------------------------------------------------------
1445
IF( ITX .EQ. 0 ) ITX = ITM
1446
ITRAT(ITX,ISTEP) = IT
1447
EFF (ITX,ISTEP) = NGOOD/CALLS*HUNDRT
1448
WRONG(ITX,ISTEP) = NEGTIV/CALLS*HUNDRT
1449
RESLT(ITX,ISTEP) = AVGI
1450
ACSTD(ITX,ISTEP) = SD
1451
TRSLT(ITX,ISTEP) = TI
1452
TACC = ABS(TSI/TI*HUNDRT)
1453
TSTD (ITX,ISTEP) = TACC
1454
PCNT (ITX,ISTEP) = ABS(SD/AVGI*HUNDRT)
1456
*----------------------------------------------------------------------
1457
* Check cumulative accuracy
1458
*----------------------------------------------------------------------
1460
IF( NODEID .EQ. 0 ) THEN
1462
*------------------- Check cumulative accuracy -----------------------
1465
IF((ABS(SDAV) .LE. ACC)) NEND = N1
1467
IF( ISTEP .EQ. N1 ) THEN
1468
IF( TACC .GT. ACLMT ) THEN
1469
IF( NER1 .EQ. 0 ) THEN
1471
WRITE(ERROR(1,NERROR),9900) NERROR,IT,ACLMT
1472
9900 FORMAT('* (',I1,') Temp. accuracy of it-#',
1473
. I3,' is too large comparing to',
1474
. F6.2,' percent.',6X,'*')
1475
WRITE(ERROR(2,NERROR),9901) TACC,ACLMT
1476
9901 FORMAT('*',8X,'Temp. accuracy (',
1478
. F7.4,' % )',23X,'*')
1479
WRITE(ERROR(3,NERROR),9902)
1480
9902 FORMAT('*',77X,'*')
1484
IF( IT .GT. 1 ) THEN
1485
IF(( TI .GT. AVTI+FDEVI ) .OR.
1486
. ( TI .LT. AVTI-FDEVI ) ) THEN
1487
IF( NER2 .EQ. 0 ) THEN
1489
WRITE(ERROR(1,NERROR),9910) NERROR,IT,FC
1490
9910 FORMAT('* (',I1,') Temp. estimate of ',
1491
. 'it-#',I3,' fluctuates more than ',
1492
. F4.1,'*average-sigma.',6X,'*')
1494
*patch TI:1995/08/25
1496
*old CALL BSORDR( RE, FX2, ORDER, IORDR )
1497
CALL BSORDR( ARE, FX2, ORDER, IORDR )
1502
*patch TI:1995/08/25
1505
IF( ARE1 .GE. AAC ) THEN
1506
CALL BSORDR( ARE1, FX2, ORDR1, IORDR1)
1508
CALL BSORDR( AAC, FX2, ORDR1, IORDR1)
1510
* IF( RE1 .GE. AC ) THEN
1511
* CALL BSORDR( RE1, FX2, ORDR1, IORDR1)
1513
* CALL BSORDR( AC, FX2, ORDR1, IORDR1)
1518
WRITE(ERROR(2,NERROR),9911) RE,IORDR,
1520
9911 FORMAT('* Temp. Estimate (',
1521
. F10.6,' E',I3,') > (',F10.6,'+',F8.6,
1522
. ' ) E',I3,', or',1X,'*')
1523
WRITE(ERROR(3,NERROR),9912) RE,IORDR,
1525
9912 FORMAT('* Temp. Estimate (',
1526
. F10.6,' E',I3,') < (',F10.6,'-',F8.6,
1531
IF( TSI .GT. FDEVI ) THEN
1532
IF( NER3 .EQ. 0 ) THEN
1534
WRITE(ERROR(1,NERROR),9920) NERROR,IT,FC
1535
9920 FORMAT('* (',I1,') Error of it-#',
1536
. I3,' fluctuates more than',F4.1,
1537
. '*average-sigma.',16X,'*')
1539
*patch TI:1995/08/25
1541
* CALL BSORDR( RE1, FX2, ORDER, IORDR)
1542
CALL BSORDR( ARE1, FX2, ORDER, IORDR)
1546
*patch TI:1995/08/25
1548
* CALL BSORDR( AC, FX2, ORDR1, IORDR1)
1549
CALL BSORDR( AAC, FX2, ORDR1, IORDR1)
1552
WRITE(ERROR(2,NERROR),9921) RE1,IORDR,
1554
9921 FORMAT('* Temp. Error (',
1555
. F10.6,' E',I3,') > (',F10.6,
1556
. ' E',I3,')',18X,'*')
1557
WRITE(ERROR(3,NERROR),9902)
1562
SUMTSI = SUMTSI + TSI
1564
AVTSI = SUMTSI/FLOAT(IT)
1565
AVTI = SUMTI/FLOAT(IT)
1570
*------------------------------------------- for a parallel processor
1573
C CALL BSCAST( NEND, 1 )
1575
*----------------------------------------------------------------------
1576
* Smoothing the Distribution D(I,J) and refine the grids
1577
*----------------------------------------------------------------------
1579
IF( ISTEP .LE. N0 ) THEN
1580
IF( IT .EQ. ITMX ) NEND = N1
1581
* ---------------------
1583
* ---------------------
1585
* ==========================
1586
CALL BSUTIM( 0, ISTEP )
1587
* ==========================
1589
TIME (ITX,ISTEP) = TIMINT
1592
*---- Print the convergency behavior table in the interactive mode ----
1593
IF( INTV .GT. 1 ) THEN
1594
* ---------------------------------
1595
CALL BSPRNT ( LU, 6, ISTEP, IDUM2 )
1596
* ---------------------------------
1599
IF( NEND .EQ. N1 ) GO TO 600
1601
* ======================
1603
* ======================
1610
***********************************************************************
1611
* Termination of BASES
1612
***********************************************************************
1616
*---------------------------------------------- For a parallel computer
1618
* Global sum of histograms
1620
* Global sum of probabilities
1621
C CALL BSDSUM( DXD, NSP )
1622
* Global sum of the max.value in each HC
1623
C CALL BSDSUM( DXP, NSP )
1626
*======================= End of the step ? ============================
1628
IF( NEND .EQ. N1 ) THEN
1629
IF( INTV .GT. 1 ) THEN
1630
* ---------------------------------
1631
CALL BSPRNT ( LU, 7, IDUM1, IDUM2 )
1632
* ---------------------------------
1634
IF( ISTEP .EQ. N0) THEN
1642
* ======================
1644
* ======================
1650
***********************************************************************
1651
* =================================== *
1652
SUBROUTINE BSLIST( LU, I, ISTEP )
1653
* =================================== *
1655
* Print out results of each iteration and cumulative result *
1658
* LU : Logical unit number for the printer *
1659
* I : Address in the arrays of common /BASE5/ *
1660
* ISTEP : The Set-Identifier *
1661
* ISTEP = ( 0 / 1 ) = ( Grid opt. / Integration step ) *
1663
* S. Kawabata March '94 *
1664
***********************************************************************
1666
IMPLICIT REAL*8 (A-H,O-Z)
1667
PARAMETER (ITM = 50)
1668
REAL*4 TIME, EFF, WRONG, TRSLT, TSTD, PCNT
1669
COMMON /BASE5/ ITRAT(ITM,0:1),TIME(ITM,0:2),EFF(ITM,0:1),
1670
. WRONG(ITM,0:1),RESLT(ITM,0:1),ACSTD(ITM,0:1),
1671
. TRSLT(ITM,0:1),TSTD(ITM,0:1),PCNT(ITM,0:1)
1673
CALL BSTCNV( TIME(I,ISTEP), IH, MN, IS1, IS2 )
1676
AC = ABS(ACSTD(I,ISTEP))
1678
IF( ARE .GE. AC) THEN
1679
CALL BSORDR( ARE, F2, ORDER, IORDR)
1681
CALL BSORDR( AC, F2, ORDER, IORDR )
1686
WRITE(LU,9631) ITRAT(I,ISTEP),IEFF,WRONG(I,ISTEP),
1687
. TRSLT(I,ISTEP),TSTD(I,ISTEP),
1688
. RE,AC,IORDR,PCNT(I,ISTEP),IH,MN,IS1,IS2
1689
9631 FORMAT(I4,I4,F6.2,1P,E11.3, 0P,1X,F6.3,
1690
. F10.6,'(+-',F8.6,')E',I3.2,1X,F6.3,
1691
. 1X,I3,':',I2,':',I2,'.',I2.2)
1698
C***********************************************************************
1700
C*============================================= *
1701
C* SUBROUTINE BSORDR( VAL, F2, ORDER, IORDR) *
1702
C*============================================= *
1704
C* To resolve the real number VAL into mantester and exponent parts.*
1705
C* When VAL = 1230.0 is given, output are *
1706
C* F2 = 1.2 and ORDER = 4.0. *
1708
C* VAL : Real*8 value *
1710
C* F2 : The upper two digits is given *
1711
C* ORDER: Order is given *
1712
C* IORDR: Exponent is given *
1716
C***********************************************************************
1718
SUBROUTINE BSORDR(VAL, F2, ORDER, IORDR)
1719
IMPLICIT REAL*8 (A-H,O-Z)
1721
IF( VAL .NE. 0.0 ) THEN
1722
ORDER = LOG10( VAL )
1723
IORDR = INT( ORDER )
1724
IF( ORDER .LT. 0.0D0 ) IORDR = IORDR - 1
1725
ORDER = 10.D0**IORDR
1737
***********************************************************************
1738
* ======================================= *
1739
SUBROUTINE BSPRNT( LU, ID, IP1, IP2 )
1740
* ======================================= *
1742
* Print out routine of BASES. *
1744
* ID : Identity number of printouts. *
1745
* IP1... IP2 : Integer *
1747
* S. Kawabata May 1992 *
1748
* Last update March 1994 *
1749
***********************************************************************
1751
IMPLICIT REAL*8 (A-H,O-Z)
1752
PARAMETER (MXDIM = 50, NDMX = 50, LENG = 32768)
1753
COMMON /BASE0/ JFLAG,IBASES
1754
COMMON /BASE1/ XL(MXDIM),XU(MXDIM),NDIM,NWILD,
1756
COMMON /BASE2/ ACC1,ACC2,ITMX1,ITMX2
1757
COMMON /BASE3/ SCALLS,WGT,TI,TSI,TACC,IT
1758
COMMON /BASE4/ XI(NDMX,MXDIM),DX(MXDIM),DXD(LENG),DXP(LENG),
1759
. ND,NG,NPG,MA(MXDIM)
1760
PARAMETER (ITM = 50)
1761
REAL*4 TIME, EFF, WRONG, TRSLT, TSTD, PCNT
1762
COMMON /BASE5/ ITRAT(ITM,0:1),TIME(ITM,0:2),EFF(ITM,0:1),
1763
. WRONG(ITM,0:1),RESLT(ITM,0:1),ACSTD(ITM,0:1),
1764
. TRSLT(ITM,0:1),TSTD(ITM,0:1),PCNT(ITM,0:1)
1766
COMMON /BSRSLT/AVGI,SD,CHI2A,STIME,IT1,ITF
1767
CHARACTER*51 ICH(0:1)
1769
* INTV = ( 0 / 1 / any ) = ( Batch / Batch(Unix) / Interactive )
1770
* IPNT = ( 0 / any ) = ( IBM Type / Ascii printer )
1771
COMMON /BSCNTL/ INTV, IPNT, NLOOP, MLOOP
1773
COMMON /BDATE/ IDATE(3),ITIME(2)
1774
* IDATE(1) : year ITIME(1) : hour
1775
* IDATE(2) : month ITIME(2) : minute
1777
REAL*4 TIMEBS,TIMINT,TIMESP,TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
1778
COMMON /BTIME1/ TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
1779
COMMON /BTIME2/ TIMEBS(0:2),TIMINT,TIMESP(0:2)
1782
COMMON/NINFO/ NODEID, NUMNOD
1784
DATA ICH / 'Convergency Behavior for the Grid Optimization Step',
1785
. 'Convergency Behavior for the Integration Step '/
1787
IF( NODEID .NE. 0 ) RETURN
1790
GO TO ( 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000 ), ID
1791
C----------------------------------------------------------- BSMAIN
1793
100 IF( IPNT .EQ. 0 ) THEN
1795
9600 FORMAT(/1H1,/1H )
1800
WRITE(LU,9620) (IDATE(I),I=1,3),(ITIME(J),J=1,2)
1801
9620 FORMAT(55X,'Date: ',I2,'/',I2,'/',I2,2X,I2.2,':',I2.2)
1804
. 8X,'**********************************************************',
1806
./8X,'* BBBBBBB AAAA SSSSSS EEEEEE SSSSSS *',
1807
./8X,'* BB BB AA AA SS SS EE SS SS *',
1808
./8X,'* BB BB AA AA SS EE SS *',
1809
./8X,'* BBBBBBB AAAAAAAA SSSSSS EEEEEE SSSSSS *',
1810
./8X,'* BB BB AA AA SS EE SS *',
1811
./8X,'* BB BB AA AA SS SS EE SS SS *',
1812
./8X,'* BBBB BB AA AA SSSSSS EEEEEE SSSSSS *',
1814
./8X,'* BASES Version 5.1 *',
1815
./8X,'* coded by S.Kawabata KEK, March 1994 *',
1816
./8X,'**********************************************************')
1819
C----------------------------------------------------------- BSMAIN
1821
200 IF( IPNT .EQ. 0 ) THEN
1828
. '****** END OF BASES *********')
1830
C----------------------------------------------------------- BSMAIN
1835
.//5X,'<< Computing Time Information >>')
1837
* WRITE(LU,9310) (IDATE(I),I=1,3),(ITIME(J),J=1,2)
1838
*9310 FORMAT(/15X,'Start at: ',I2,'/',I2,'/',I2,2X,I2.2,':',I2.2)
1840
* WRITE(LU,9320) (IDATE(I),I=1,3),(ITIME(J),J=1,2)
1841
*9320 FORMAT(15X,'End at: ',I2,'/',I2,'/',I2,2X,I2.2,':',I2.2)
1843
9330 FORMAT(/15X,'(1) For BASES H: M: Sec')
1844
CALL BSTCNV(TIMEBS(2),IH,MN,IS1,IS2)
1845
WRITE(LU,9340) IH, MN, IS1, IS2
1846
9340 FORMAT(19X,'Overhead : ',I3,':',I2,':',I2,'.',I2.2)
1847
CALL BSTCNV(TIMEBS(0),IH,MN,IS1,IS2)
1848
WRITE(LU,9350) IH, MN, IS1, IS2
1849
9350 FORMAT(19X,'Grid Optim. Step : ',I3,':',I2,':',I2,'.',I2.2)
1850
CALL BSTCNV(TIMEBS(1),IH,MN,IS1,IS2)
1851
WRITE(LU,9360) IH, MN, IS1, IS2
1852
9360 FORMAT(19X,'Integration Step : ',I3,':',I2,':',I2,'.',I2.2)
1853
XTIME = TIMEB2 - TIMEB1
1854
CALL BSTCNV(XTIME,IH,MN,IS1,IS2)
1855
WRITE(LU,9365) IH, MN, IS1, IS2
1856
9365 FORMAT(19X,'Go time for all : ',I3,':',I2,':',I2,'.',I2.2)
1857
EXTIM = TIMEBS(1)*1000.0/SCALLS/0.7
1859
9375 FORMAT(/15X,'(2) Expected event generation time')
1860
WRITE(LU,9376) EXTIM
1861
9376 FORMAT(19X,'Expected time for 1000 events :',F10.2,' Sec')
1864
C----------------------------------------------------------- BASES
1868
WRITE(LU,9400) NDIM,NWILD,MCALL,NCALL,ND,NG,NSP
1870
.//5X,'<< Parameters for BASES >>',
1871
.//5X,' (1) Dimensions of integration etc.',
1872
. /5X,' # of dimensions : Ndim =',I9,3X,'( 50 at max.)',
1873
. /5X,' # of Wilds : Nwild =',I9,3X,'( 15 at max.)',
1874
. /5X,' # of sample points : Ncall =',I9,'(real)',
1876
. /5X,' # of subregions : Ng =',I9,' / variable',
1877
. /5X,' # of regions : Nregion =',I9,' / variable',
1878
. /5X,' # of Hypercubes : Ncube =',I9,
1879
.//5X,' (2) About the integration variables')
1881
9405 FORMAT(10X,'------',2('+---------------'),'+-------+-------')
1883
9410 FORMAT(10X,' i XL(i) XU(i) ',
1887
IF( I .LE. NWILD ) THEN
1888
WRITE(LU,9420) I,XL(I),XU(I),IG(I)
1889
9420 FORMAT(10X,I5,1P,2(' ',E14.6),' ',3X,0P,I1,3X,
1892
WRITE(LU,9421) I,XL(I),XU(I),IG(I)
1893
9421 FORMAT(10X,I5,1P,2(' ',E14.6),' ',3X,0P,I1,3X,
1898
WRITE(LU,9450) ITMX1,ACC1,ITMX2,ACC2
1900
. /5X,' (3) Parameters for the grid optimization step',
1901
. /5X,' Max.# of iterations: ITMX1 =',I9,
1902
. /5X,' Expected accuracy : Acc1 =',F9.4,' %',
1903
.//5X,' (4) Parameters for the integration step',
1904
. /5X,' Max.# of iterations: ITMX2 =',I9,
1905
. /5X,' Expected accuracy : Acc2 =',F9.4,' %')
1908
C----------------------------------------------------------- BASES
1910
500 IF( INTV .LE. 1 ) RETURN
1912
IF( IPNT .EQ. 0 ) THEN
1917
WRITE(LU,9620) (IDATE(I),I=1,3),(ITIME(J),J=1,2)
1918
WRITE(LU,9500) ICH(ISTEP)
1922
9550 FORMAT(1X,'<- Result of each iteration ->',
1923
. 2X,'<- Cumulative Result ->',
1924
. 1X,'< CPU time >',
1925
. /1X,' IT Eff R_Neg Estimate Acc %',
1926
. 2X,'Estimate(+- Error )order Acc %',
1927
. 1X,'( H: M: Sec )')
1929
9570 FORMAT(1X,7('----------'),'--------')
1932
C----------------------------------------------------------- BASES
1934
600 IF( INTV .LE. 1 ) RETURN
1937
IF( ITX .EQ. 0 ) ITX = ITM
1939
CALL BSLIST( LU, ITX, ISTEP )
1943
700 IF( INTV .LE. 1 ) RETURN
1947
C----------------------------------------------------------- BASES
1951
ITX = MOD( ITJ, ITM )
1952
IF( ITX .EQ. 0 ) ITX = ITM
1954
IF( ITRAT(1,ISTEP) .EQ. 1 ) THEN
1961
IF( ITRAT(I,ISTEP) .LT. ITMN ) THEN
1963
ITMN = ITRAT(I,ISTEP)
1966
IF( ITST .EQ. 1 ) NDEV = 1
1969
IF( IPNT .EQ. 0 ) THEN
1974
WRITE(LU,9620) (IDATE(I),I=1,3),(ITIME(J),J=1,2)
1975
WRITE(LU,9500) ICH(ISTEP)
1980
625 IF( NDEV .EQ. 1 ) THEN
1985
DO 650 I = ITST, ITFN
1987
CALL BSLIST( LU, I, ISTEP )
1991
IF( NDEV .GT. 0 ) GO TO 625
1996
C----------------------------------------------------------- BASES
1999
9950 FORMAT(1X,'******** FATAL ERROR IN BASES **************',
2000
. /1X,'There are no enough good points in this iteration.',
2001
. /1X,'Process was terminated due to this error.')
2005
C-----------------------------------------------------------------
2007
IF( IP2 .NE. 0 ) THEN
2008
IF( IPNT .EQ. 0 ) THEN
2013
WRITE(LU,9620) (IDATE(I),I=1,3),(ITIME(J),J=1,2)
2016
. 20X,'Results of Integration',
2017
. /10X,5('----------'),'------',
2018
. /10X,' Loop# Estimate(+- Error )order',
2019
. ' It1 It2 ( H: M: Sec )',
2020
. /10X,5('----------'),'------')
2026
IF( ARE .GE. AC) THEN
2027
CALL BSORDR( ARE, F2, ORDER, IORDR)
2029
CALL BSORDR( AC, F2, ORDER, IORDR )
2033
CALL BSTCNV( STIME, IH, MN, IS1, IS2)
2034
WRITE(LU,9660) LOOP,RE,AC,IORDR,IT1,IT,IH,MN,IS1,IS2
2035
9660 FORMAT(10X,I6,F10.6,'(+-',F8.6,')E',I3.2,2I5,
2036
. 1X,I3,':',I2,':',I2,'.',I2.2,
2037
. /10X,5('----------'),'------')
2043
C***********************************************************************
2045
C*======================== *
2046
C* SUBROUTINE BSPUTW( WEIGHT ) *
2047
C*======================== *
2051
C* Coded by T.Ishikawa Jun. 1995 at KEK *
2052
C* Last update Jun. 1995 at KEK *
2054
C***********************************************************************
2056
SUBROUTINE BSPUTW( WEIGHT )
2058
IMPLICIT REAL*8 (A-H,O-Z)
2059
COMMON /BASE3/ SCALLS,WGT,TI,TSI,TACC,IT
2061
*========= Save the grid information for the best accuracy ===========
2069
************************************************************************
2071
* ========================== *
2072
SUBROUTINE BSREAD( LUN )
2073
* ========================== *
2075
* Read temporary result from the logocal unit LUN *
2077
* S.Kawabata June '90 at KEK *
2079
************************************************************************
2082
IMPLICIT REAL*8 (A-H,O-Z)
2083
PARAMETER (MXDIM = 50, NDMX = 50, LENG = 32768)
2084
COMMON /BASE1/ ND1(5*MXDIM+3)
2085
* COMMON /BASE1/ XL(MXDIM),XU(MXDIM),NDIM,NWILD,
2087
C COMMON /BASE2/ ND2(6)
2088
* COMMON /BASE2/ ACC1,ACC2,ITMX1,ITMX2
2089
COMMON /BASE3/ ND3(11)
2090
* COMMON /BASE3/ SCALLS,WGT,TI,TSI,TACC,IT
2091
COMMON /BASE4/ ND4(2*MXDIM*(NDMX+1)+4*LENG+MXDIM+3)
2092
* COMMON /BASE4/ XI(NDMX,MXDIM),DX(MXDIM),DXD(LENG),DXP(LENG),
2093
* . ND,NG,NPG,MA(MXDIM)
2094
PARAMETER (ITM = 50 )
2095
* COMMON /BASE5/ ND5(22*ITM)
2096
COMMON /BASE5/ ND5(23*ITM)
2097
* REAL*4 TIME, EFF, WRONG, TRSLT, TSTD, PCNT
2098
* COMMON /BASE5/ ITRAT(ITM,0:1),TIME(ITM,0:2),EFF(ITM,0:1),
2099
* . WRONG(ITM,0:1),RESLT(ITM,0:1),ACSTD(ITM,0:1),
2100
* . TRSLT(ITM,0:1),TSTD(ITM,0:1),PCNT(ITM,0:1)
2101
COMMON /RANDM/ ND6(45)
2103
PARAMETER ( NHS = 50, NSC = 50 )
2104
COMMON /PLOTH/ NPH(18*(NHS+NSC)+29),NW
2105
COMMON /PLOTB/ IBUF( 281*NHS + 2527*NSC )
2106
* INTEGER*4 XHASH,DHASH,NHIST,MAPL,IFBASE,NSCAT,MAPD
2107
* COMMON/PLOTH/ XHASH(NHS+1,13),DHASH(NSC+1,14),IFBASE(NHS),
2108
* . NHIST, MAPL(4,NHS),
2109
* . NSCAT, MAPD(4,NSC),
2112
COMMON/NINFO/ NODEID, NUMNOD
2114
IF( NODEID .NE. 0 ) RETURN
2117
READ(LUN) ND1,ND3,ND4,ND5,ND6,NPH
2118
C READ(LUN) ND1,ND2,ND3,ND4,ND5,ND6,NPH
2120
READ(LUN) NW,(IBUF(I),I=1,NW)
2126
************************************************************************
2127
*=================================================
2128
SUBROUTINE BSTCNV( TIME, IH, MN, IS1, IS2 )
2129
*=================================================
2131
* Resolve TIME in second into IH, MN, IS1, IS2
2133
* TIME : in the unit of second
2140
* S.Kawabata 1992 June 15
2141
************************************************************************
2143
IMPLICIT REAL*8 (A-H,O-Z)
2146
DATA HOUR, MINUT, N100/ 360000, 6000, 100 /
2151
IF( ISEC .GE. MINUT ) THEN
2153
IF( ISEC .GE. HOUR ) THEN
2160
ISEC = ISEC - MN*MINUT
2163
IS2 = MOD( ISEC, N100)
2169
SUBROUTINE BSUTIM( JOB, ID )
2171
C COMMON/NINFO/ NODEID, NUMNOD
2172
COMMON /BTIME1/ TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
2173
COMMON /BTIME2/ TIMEBS(0:2),TIMINT,TIMESP(0:2)
2175
* Prior to call thisroutine, BSTIME( TIME0, 1 ) should be called
2176
* for initialize the time offset TIME0.
2178
* print *,'bsutim .. job, id ',job,id
2179
CALL BSTIME( RTIME, 1)
2180
DTIME = RTIME - TIME0
2182
IF( JOB .EQ. 0 ) THEN
2183
* For BASES computing time
2184
* ID = 0 : Grid defining step
2185
* 1 : Integration step
2188
TIMEBS(ID) = TIMEBS(ID) + DTIME
2190
IF( ID .LE. 1 ) THEN
2191
TIMINT = TIMINT + DTIME
2194
* For SPRING computing time
2195
* ID = 0 : Event generation
2199
TIMESP(ID) = TIMESP(ID) + DTIME
2209
************************************************************************
2211
* ========================== *
2212
SUBROUTINE BSWRIT( LUN )
2213
* ===================== *
2215
* Read temporary result from disk file. *
2217
* S.Kawabata June '90 at KEK *
2219
************************************************************************
2221
IMPLICIT REAL*8 (A-H,O-Z)
2222
PARAMETER (MXDIM = 50, NDMX = 50, LENG = 32768)
2223
COMMON /BASE1/ ND1(5*MXDIM+3)
2224
* COMMON /BASE1/ XL(MXDIM),XU(MXDIM),NDIM,NWILD,
2226
C COMMON /BASE2/ ND2(6)
2227
* COMMON /BASE2/ ACC1,ACC2,ITMX1,ITMX2
2228
COMMON /BASE3/ ND3(11)
2229
* COMMON /BASE3/ SCALLS,WGT,TI,TSI,TACC,IT
2230
COMMON /BASE4/ ND4(2*MXDIM*(NDMX+1)+4*LENG+MXDIM+3)
2231
* COMMON /BASE4/ XI(NDMX,MXDIM),DX(MXDIM),DXD(LENG),DXP(LENG),
2232
* . ND,NG,NPG,MA(MXDIM)
2233
PARAMETER (ITM = 50 )
2234
* COMMON /BASE5/ ND5(22*ITM)
2235
COMMON /BASE5/ ND5(23*ITM)
2236
* REAL*4 TIME, EFF, WRONG, TRSLT, TSTD, PCNT
2237
* COMMON /BASE5/ ITRAT(ITM,0:1),TIME(ITM,0:2),EFF(ITM,0:1),
2238
* . WRONG(ITM,0:1),RESLT(ITM,0:1),ACSTD(ITM,0:1),
2239
* . TRSLT(ITM,0:1),TSTD(ITM,0:1),PCNT(ITM,0:1)
2240
COMMON /RANDM/ ND6(45)
2242
PARAMETER ( NHS = 50, NSC = 50 )
2243
COMMON /PLOTH/ NPH(18*(NHS+NSC)+29),NW
2244
COMMON /PLOTB/ IBUF( 281*NHS + 2527*NSC )
2245
* INTEGER*4 XHASH,DHASH,NHIST,MAPL,IFBASE,NSCAT,MAPD
2246
* COMMON/PLOTH/ XHASH(NHS+1,13),DHASH(NSC+1,14),IFBASE(NHS),
2247
* . NHIST, MAPL(4,NHS),
2248
* . NSCAT, MAPD(4,NSC),
2251
COMMON/NINFO/ NODEID, NUMNOD
2253
IF( NODEID .NE. 0 ) RETURN
2256
WRITE(LUN) ND1,ND3,ND4,ND5,ND6,NPH
2257
C WRITE(LUN) ND1,ND2,ND3,ND4,ND5,ND6,NPH
2258
IF(NW .EQ. 0 ) NW = 281
2259
WRITE(LUN) NW,(IBUF(I),I=1,NW)
2265
C**********************************************************************
2266
C*====================== *
2267
C* FUNCTION DRN( ISEED) *
2268
C*====================== *
2269
C* Machine-independent Random number generator *
2270
C* General purpose Version, OK as long as >= 32 bits *
2274
C**********************************************************************
2276
* REAL FUNCTION DRN*8(ISEED)
2277
DOUBLE PRECISION FUNCTION DRN(ISEED)
2279
COMMON/RANDM/RDM(31),RM1,RM2,IA1,IC1,M1,IX1,
2283
C Generate Next number in sequence
2285
IX1 = MOD( IA1*IX1+IC1, M1 )
2286
IX2 = MOD( IA2*IX2+IC2, M2 )
2287
IX3 = MOD( IA3*IX3+IC3, M3 )
2290
RDM(J) = ( FLOAT(IX1)+FLOAT(IX2)*RM2 )*RM1
2292
C Omit following statement if function arguement passed by value:
2299
C**********************************************************************
2300
C*============================ *
2301
C* Subroutine DRNSET( ISEED ) *
2302
C*============================ *
2304
C* Initialization routine of *
2305
C* Machine-independent Random number generator *
2306
C* General purpose Version, OK as long as >= 32 bits *
2310
C**********************************************************************
2312
SUBROUTINE DRNSET( ISEED )
2314
COMMON/RANDM/RDM(31),RM1,RM2,IA1,IC1,M1,IX1,
2330
IX1 = MOD( ISEED, M1 )
2331
IX1 = MOD( IA1*IX1+IC1, M1 )
2332
IX2 = MOD( IX1, M2 )
2333
IX1 = MOD( IA1*IX1+IC1, M1 )
2338
IX1 = MOD( IA1*IX1+IC1, M1 )
2339
IX2 = MOD( IA2*IX2+IC2, M2 )
2340
RDM(J)= ( FLOAT(IX1)+FLOAT(IX2)*RM2 )*RM1
2347
c End of BASES routines
2352
c Begin of SPRING routines
2355
subroutine init_of_spring(prefix)
2356
c This subroutine is called in order to initialize spring parameters
2359
common/sprng0/ispring
2360
character * 60 prefix,fname
2362
parameter (isunit=23)
2364
call fk88strcat(prefix,'_bs.data',fname)
2365
open(unit=isunit,file=fname,
2366
# form='unformatted',status='old')
2374
subroutine run_of_spring(ff,mxevts,mxtrls,nevts,ntrls)
2377
integer mxevts,mxtrls,nevts,ntrls
2379
parameter (mxtry=50)
2381
parameter (iounit=6)
2382
integer mxtryp,nevent,ntrial,miss
2383
common/sprng2/mxtryp,nevent,ntrial,miss
2388
dowhile(ntrial.lt.mxtrls.and.nevent.lt.mxevts)
2389
call spring(ff,mxtry,mxevts)
2399
c It appears that the main unweighting routine is sprgen; spring
2400
c seems only to be a wrapper for the former. In sprgen the information
2401
c is available on the values of xx(i) when the event is accepted.
2403
c There is a serious logical mistake in spring: the initialization
2404
c is performed only if ibases is different from zero. Ibases is
2405
c a flag which is set to 1 in bsinit, and it is connected to bases,
2406
c and not to spring. Thus, if one runs spring by using the bases data
2407
c files saved at a earlier stage, spring is initialized only the first
2408
c time, and this could force the code to crash, but it could also
2409
c lead to incorrect results without notice.
2411
c A new common block /sprng0/ispring has therefore been introduced,
2412
c which controls the initialization of spring.
2413
************************************************************************
2414
* ================================== *
2415
SUBROUTINE SPRING(FUNC, MXTRY, MXEVTS )
2416
* ================================== *
2417
* Main Program for the Event generation program SPRING. *
2419
* Coded by S.Kawabata September '84 *
2421
* Modified by S. Frixione in order to exclude the built-in *
2422
* histogramming package. The common block bfxlhisto has been added *
2423
* The common block SPRNG0 has been added, to correct the logical *
2424
* flaw pointed out in the comments above. A further modification *
2425
* occurred on 6/4/2005. Entry MXEVTS has been added to this routine, *
2426
* in order to modify the warning relevant to the number of *
2427
* mis-generations, now issued only after the mis-generations are *
2428
* than 1% of the total number events to be generated *
2429
************************************************************************
2431
IMPLICIT REAL*8 (A-H,O-Z)
2433
COMMON/SPRNG0/ISPRING
2434
COMMON /BASE0/ NDUM,IBASES
2435
PARAMETER (MXDIM = 50, NDMX = 50, LENG = 32768)
2436
COMMON /BASE1/ XL(MXDIM),XU(MXDIM),NDIM,NWILD,
2438
COMMON /BASE4/ XI(NDMX,MXDIM),DX(MXDIM),DXD(LENG),DXP(LENG),
2439
. ND,NG,NPG,MA(MXDIM)
2440
COMMON /BDATE/ IDATE(3),ITIME(2)
2441
COMMON /BSCNTL/ INTV, IPNT, NLOOP, MLOOP
2443
COMMON /SPRNG1/ XND, DXG, XJAC, DXMAX, NSP
2444
COMMON /SPRNG2/ MXTRYP,NEVENT, NTRIAL,MISS
2446
REAL*4 TIMEBS,TIMINT,TIMESP,TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
2447
COMMON /BTIME1/ TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
2448
COMMON /BTIME2/ TIMEBS(0:2),TIMINT,TIMESP(0:2)
2450
COMMON/BFXLHISTO/LHISTO
2452
*----------------------------- Entry point ----------------------------*
2454
*======================================================================*
2455
* Initialization of the program *
2456
*======================================================================*
2457
*----------------------------------------------------------------------*
2458
* initialize timer etc. *
2459
*----------------------------------------------------------------------*
2461
c IF( IBASES .GT. 0 ) THEN
2462
IF( ISPRING .GT. 0 ) THEN
2466
CALL BSTIME( TIME0, 0 )
2481
IF( MXTRY .LT. 10 ) MXTRY = 50
2483
IF( MXTRY .GT. 50) NBIN = 50
2491
CALL SHINIT( MXTRY1 )
2494
*----------------------------------------------------------------------*
2495
* Make the cumulative probability distribution *
2496
*----------------------------------------------------------------------*
2504
* CALL BSPRNT( 4, MCALL, IDUM2, IDUM3, IDUM4 )
2512
IF( DXD( I ) .LT. 0.0D0 ) THEN
2515
. /5X,'********** FATAL ERROR IN SPRING **********',
2516
. /5X,'* Negative probability was found *',
2517
. /5X,'* in the ',I6,'-th Hypercube. *',
2518
. /5X,'*******************************************')
2521
DXMAX = DXMAX + DXD( I )
2524
* =====================
2526
* =====================
2528
* =====================
2530
* =====================
2531
c IF( IBASES .EQ. 1 ) THEN
2532
IF( ISPRING .EQ. 1 ) THEN
2535
. 1X,'**************************************************',
2536
. /1X,'* Flag ISPRING was not equal to "0". *',
2538
. /1X,'* Process was terminated by this error. *',
2539
. /1X,'* Call S.Kawabata. *',
2540
. /1X,'**************************************************')
2544
*======================================================================*
2545
* Event generation *
2546
*======================================================================*
2547
* =====================
2548
500 CALL BSUTIM( 1, 1 )
2549
* =====================
2551
* ==================================
2552
CALL SPRGEN( FUNC, MXTRY, IRET)
2553
* ==================================
2555
* =====================
2557
* =====================
2563
IF( IRET .LE. MXTRY ) THEN
2564
NTRIAL =NTRIAL + IRET
2570
NTRIAL =NTRIAL + IRET - 1
2572
IF( MISFLG .EQ. 0 .AND. MISS .GT. INT(0.01*MXEVTS) ) THEN
2574
9600 FORMAT(1X,'****************************************',
2575
. '****************************************',
2576
. /1X,'* (((( Warning )))) ',
2580
. /1X,'* The number of mis-generations is foun',
2581
. 'd more than',I3,' times. *')
2583
9610 FORMAT(1X,'* ',
2585
. /1X,'*(( Suggestion )) ',
2587
. /1X,'* (1) Try integration again with larger ',
2588
. 'number of sample points than this job. *',
2591
. /1X,'* (2) The integral variables are not sui',
2592
. 'ted for the function. *',
2593
. /1X,'* Take another integral variables !!',
2597
. /1X,'****************************************',
2598
. '****************************************')
2603
* =====================
2604
600 CALL BSUTIM( 1, 1 )
2605
* =====================
2611
************************************************************************
2612
* =================== *
2614
* =================== *
2616
* To check user's initialization parameters. *
2618
* Coded by S.Kawabata April '94 *
2620
************************************************************************
2622
IMPLICIT REAL*8 (A-H,O-Z)
2623
PARAMETER ( MXDIM = 50)
2624
COMMON /BPARM1/ XL(MXDIM),XU(MXDIM),NDIM,NWILD,
2626
COMMON /BPARM2/ ACC1,ACC2,ITMX1,ITMX2
2628
COMMON /BASE0/ JFLAG,IBASES
2629
COMMON /BASE1/ XLT(MXDIM),XUT(MXDIM),NDIMT,NWILDT,
2631
COMMON /BASE2/ ACC1T,ACC2T,ITMX1T,ITMX2T
2632
COMMON /BSCNTL/ INTV, IPNT, NLOOP, MLOOP
2634
IF( NDIM .NE. NDIMT ) THEN
2635
WRITE(6,9100) NDIM,NDIMT
2637
. 5X,'*************************************************',
2639
. /5X,'* Given NDIM(',I6,' ) does not match *',
2640
. /5X,'* to NDIM(',I6,' ) in BASES. *',
2642
. /5X,'* Process was terminated due to this error. *',
2644
. /5X,'*************************************************')
2648
IF( NWILD .NE. NWILDT ) THEN
2649
WRITE(6,9110) NWILD,NWILDT
2651
. 5X,'*************************************************',
2653
. /5X,'* Given NWILD(',I6,' ) does not match *',
2654
. /5X,'* to NWILD(',I6,' ) in BASES. *',
2656
. /5X,'* Process was terminated due to this error. *',
2658
. /5X,'*************************************************')
2663
IF( XL(I) .NE. XLT(I) ) THEN
2664
WRITE(6,9200) I,XL(I),I,XLT(I)
2666
. 5X,'*************************************************',
2668
. /5X,'* Given XL(',I3,' ) = ',D15.8,' *',
2669
. /5X,'* does not match to *',
2670
. /5X,'* to XL(',I3,' ) = ',D15.8,' in BASES *',
2672
. /5X,'* Process was terminated due to this error. *',
2674
. /5X,'*************************************************')
2677
IF( XU(I) .NE. XUT(I) ) THEN
2678
WRITE(6,9210) I,XU(I),I,XUT(I)
2680
. 5X,'*************************************************',
2682
. /5X,'* Given XU(',I3,' ) = ',D15.8,' *',
2683
. /5X,'* does not match to *',
2684
. /5X,'* to XU(',I3,' ) = ',D15.8,' in BASES *',
2686
. /5X,'* Process was terminated due to this error. *',
2688
. /5X,'*************************************************')
2697
***********************************************************************
2698
*============================ *
2699
SUBROUTINE SPINFO( LU )
2700
*============================ *
2702
* Print the information for *
2703
* (1) BASES parameters *
2704
* (2) Computer time information *
2705
* (3) Convergency behavior of the Grid optimization step *
2706
* (4) Convergency behavior of the integration step *
2708
* LU : Logical unit number of printer *
2710
* by S.Kawabata March 1994 at KEK
2712
* Modified by S. Frixione in order to exclude the built-in *
2713
* histogramming package. The common block bfxlhisto has been added *
2714
***********************************************************************
2716
IMPLICIT REAL*8 (A-H,O-Z)
2717
COMMON /BDATE/ IDATE(3),ITIME(2)
2718
COMMON /BSCNTL/ INTV, IPNT, NLOOP, MLOOP
2720
COMMON /SPRNG2/ MXTRY,NEVENT, NTRIAL, MISS
2722
PARAMETER ( NHS = 50, NSC = 50 )
2723
INTEGER*4 XHASH,DHASH,NHIST,MAPL,IFBASE,NSCAT,MAPD
2724
* COMMON/PLOTH/ XHASH(ILH,13),DHASH(IDH,14),IFBASE(ILH),
2725
* . MAXL, NHIST, MAPL(4,ILH),
2726
* . MAXD, NSCAT, MAPD(4,IDH),
2728
COMMON/PLOTH/ XHASH(NHS+1,13),DHASH(NSC+1,14),IFBASE(NHS),
2729
. NHIST, MAPL(4,NHS),
2730
. NSCAT, MAPD(4,NSC),
2733
REAL*4 TIMEBS,TIMINT,TIMESP,TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
2734
COMMON /BTIME1/ TIME0,RTIME,TIMEB1,TIMEB2,TIMES1
2735
COMMON /BTIME2/ TIMEBS(0:2),TIMINT,TIMESP(0:2)
2738
COMMON/BFXLHISTO/LHISTO
2742
IF( IPNT .EQ. 0 ) THEN
2748
9300 FORMAT(/1H1,////1H )
2749
9350 FORMAT(A1,////1X)
2750
WRITE(LU,9360) (IDATE(I),I=1,3),(ITIME(J),J=1,2)
2751
9360 FORMAT(55X,'Date: ',I2,'/',I2,'/',I2,2X,I2.2,':',I2.2)
2754
. 8X,'**********************************************************',
2756
./8X,'* SSSSS PPPPPP RRRRRR IIIII N NN GGGGG *',
2757
./8X,'* SS SS PP PP RR RR III NN NN GG GG *',
2758
./8X,'* SS PP PP RR RR III NNN NN GG *',
2759
./8X,'* SSSSS PPPPPP RRRRR III NNNN NN GG GGGG *',
2760
./8X,'* SS PP RR RR III NN NNNN GG GG *',
2761
./8X,'* SS SS PP RR RR III NN NNN GG GG *',
2762
./8X,'* SSSSS PP RR RR IIIII NN NN GGGGG *',
2764
./8X,'* SPRING Version 5.1 *',
2765
./8X,'* coded by S.Kawabata KEK, March 1994 *',
2766
./8X,'**********************************************************')
2768
EFF = FLOAT(NEVENT)/FLOAT(NTRIAL)*100.D0
2769
CALL BSTIME( RTIME, 1 )
2770
XTIME = RTIME - TIMES1
2771
WRITE(LU,9500) NEVENT,EFF,(TIMESP(I),I=0,2),XTIME,MXTRY,MISS
2772
9500 FORMAT(/5X,'Number of generated events =',I10,
2773
. /5X,'Generation efficiency =',F10.3,' Percent',
2774
. /5X,'Computing time for generation =',F10.3,' Seconds',
2775
. /5X,' for Overhead =',F10.3,' Seconds',
2776
. /5X,' for Others =',F10.3,' Seconds',
2777
. /5X,'GO time for event generation =',F10.3,' Seconds',
2778
. /5X,'Max. number of trials MXTRY =',I10,' per event',
2779
. /5X,'Number of mis-generations =',I10,' times')
2789
C***********************************************************************
2790
C*==================================== *
2791
C* SUBROUTINE SPRGEN( F, MXTRY, NTRY ) *
2792
C*==================================== *
2794
C* Generation of events according to the probability density *
2795
C* which is stored in a disk file. *
2797
C* Coded by S.Kawabata at July,1980 *
2798
C* Update S.Kawabata September '84 *
2800
C* Modified by S. Frixione in order to exclude the built-in *
2801
C* histogramming package. The common block bfxlhisto has been added *
2802
C***********************************************************************
2804
SUBROUTINE SPRGEN(F,MXTRY,NTRY)
2806
IMPLICIT REAL*8 (A-H,O-Z)
2809
PARAMETER (MXDIM = 50, NDMX = 50, LENG = 32768)
2810
COMMON /BASE1/ XL(MXDIM),XU(MXDIM),NDIM,NWILD,
2812
COMMON /BASE4/ XI(NDMX,MXDIM),DX(MXDIM),DXD(LENG),DXP(LENG),
2813
. ND,NG,NPG,MA(MXDIM)
2815
COMMON /SPRNG1/ XND, DXG, XJAC, DXMAX, NSP
2818
COMMON/BFXLHISTO/LHISTO
2819
DIMENSION Y(MXDIM),KG(MXDIM)
2823
RX = DRN(IDUMY)*DXMAX
2825
C -------------- Binary Search --------------------------------
2830
300 IC = (IPMIN+IPMAX)/2
2831
IF(RX .LT. DXD(IC)) THEN
2836
IF(IPMAX-IPMIN .GT. 2) GO TO 300
2840
IF(DXD(IC) .LT. RX) GO TO 350
2842
C --------------------------------------------------------------------
2843
C Identify the hypecube number from sequential number IC
2844
C --------------------------------------------------------------------
2850
KG(NWILD) = IX/MA(NWILD) + 1
2851
IF( NWILD .GT. 1 ) THEN
2852
DO 400 J = 1,NWILD-1
2853
NUM = MOD(IX,MA(J+1))
2854
KG(J) = NUM/MA(J) + 1
2858
C ------------------------------------------------------------------
2859
C Sample and test a event
2860
C ------------------------------------------------------------------
2862
DO 600 NTRY = 1,MXTRY
2865
IF( J .LE. NWILD) THEN
2866
XN = (KG(J)-DRN(IDUMY))*DXG+ONE
2868
XN = ND*DRN(IDUMY) + ONE
2875
XO = XI(IAJ,J)-XI(IAJ-1,J)
2876
RC = XI(IAJ-1,J)+(XN-IAJ)*XO
2878
Y(J) = XL(J) + RC*DX(J)
2887
IF( FX .GT. 0.0D0 ) THEN
2888
* IF( DRN(IDUMY) .LE. FUNCT ) GO TO 700
2890
IF( XJ .LE. FUNCT ) GO TO 700
2891
* IF( XJ .LE. FUNCT ) THEN
2892
* WRITE(6,9999) NTRY,IC,FF,WGT,XJ,FUNCT
2893
*9999 FORMAT(1X,'NTRY,IC,FF,WGT,XJ,FUNCT = ',2I5,4E12.4)
2897
. IF( FX .LT. 0.0D0 ) THEN
2900
. /5X,'********** FATAL ERROR IN SPRING **********',
2901
. /5X,'* A negative value of function was found *',
2902
. /5X,'* in the ',I6,'-th Hypercube. *',
2903
. /5X,'*******************************************')
2905
9405 FORMAT(5X,'------',3('+---------------'),'+')
2907
9410 FORMAT(5X,' i XL(i) X ',
2911
WRITE(6,9420) I,XL(I),Y(I),XU(I)
2912
9420 FORMAT(5X,I5,1P,3(' ',E14.6))
2930
c End of SPRING routines
2935
c Begin of utilities
2938
subroutine shinit(mxtry1)
2939
write(6,*)'shinit is not here'
2945
write(6,*)'shrset is not here'
2950
subroutine shfill(iret)
2951
write(6,*)'shfill is not here'
2957
write(6,*)'shupdt is not here'
2962
subroutine sphist(lu)
2963
write(6,*)'sphist is not here'
2969
write(6,*)'shcler is not here'
2974
subroutine bhinit(lu)
2975
write(6,*)'bhinit is not here'
2981
write(6,*)'bhrset is not here'
2987
write(6,*)'bhsave is not here'
added
removed
Template/NLO/MCatNLO/srcCommon/mcatnlo_int.f
