| 1 | c second version of vegas with double precision
|
|---|
| 2 | c...............
|
|---|
| 3 |
|
|---|
| 4 | SUBROUTINE vegas(region,ndim,fxn,init,ncall,itmx,nprn,
|
|---|
| 5 | > tgral,sd,chi2a)
|
|---|
| 6 | INTEGER init,itmx,ncall,ndim,nprn,NDMX,MXDIM
|
|---|
| 7 | REAL*8 tgral,chi2a,sd,region(2*ndim),fxn,ALPH,TINY
|
|---|
| 8 | PARAMETER (ALPH=1.5,NDMX=50,MXDIM=10,TINY=1.d-30)
|
|---|
| 9 | EXTERNAL fxn
|
|---|
| 10 | C USES fxn,ran2,rebin
|
|---|
| 11 | C
|
|---|
| 12 | INTEGER i,idum,it,j,k,mds,nd,ndo,ng,npg,ia(MXDIM),kg(MXDIM)
|
|---|
| 13 | REAL*8 calls,dv2g,dxg,f,f2,f2b,fb,rc,ti,tsi,wgt,xjac,xn,xnd,xo,
|
|---|
| 14 | > d(NDMX,MXDIM),di(NDMX,MXDIM),dt(MXDIM),dx(MXDIM),
|
|---|
| 15 | > r(NDMX),x(MXDIM),xi(NDMX,MXDIM),xin(NDMX),ran2
|
|---|
| 16 | DOUBLE PRECISION schi,si,swgt
|
|---|
| 17 | COMMON /ranno/ idum
|
|---|
| 18 | SAVE
|
|---|
| 19 | if(init.le.0)then
|
|---|
| 20 | mds=1
|
|---|
| 21 | ndo=1
|
|---|
| 22 | do 11 j=1,ndim
|
|---|
| 23 | xi(1,j)=1
|
|---|
| 24 | 11 enddo
|
|---|
| 25 | endif
|
|---|
| 26 | if (init.le.1) then
|
|---|
| 27 | si=0
|
|---|
| 28 | swgt=0
|
|---|
| 29 | schi=0
|
|---|
| 30 | endif
|
|---|
| 31 | if (init.le.2)then
|
|---|
| 32 | nd=NDMX
|
|---|
| 33 | ng=1
|
|---|
| 34 | if(mds.ne.0)then
|
|---|
| 35 | ng=(ncall/2.d0+0.25d0)**(1.d0/ndim)
|
|---|
| 36 | mds=1
|
|---|
| 37 | if((2*ng-NDMX).ge.0)then
|
|---|
| 38 | mds=-1
|
|---|
| 39 | npg=ng/NDMX+1
|
|---|
| 40 | nd=ng/npg
|
|---|
| 41 | ng=npg*nd
|
|---|
| 42 | endif
|
|---|
| 43 | endif
|
|---|
| 44 | k=ng**ndim
|
|---|
| 45 | npg=max(ncall/k,2)
|
|---|
| 46 | calls=npg*k
|
|---|
| 47 | dxg=1.d0/ng
|
|---|
| 48 | dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-1.d0)
|
|---|
| 49 | xnd=nd
|
|---|
| 50 | dxg=dxg*xnd
|
|---|
| 51 | xjac=1.d0/calls
|
|---|
| 52 | do 12 j=1, ndim
|
|---|
| 53 | dx(j)=region(j+ndim)-region(j)
|
|---|
| 54 | xjac=xjac*dx(j)
|
|---|
| 55 | 12 enddo
|
|---|
| 56 | if(nd.ne.ndo)then
|
|---|
| 57 | do 13 i=1,nd
|
|---|
| 58 | r(i)=1.d0
|
|---|
| 59 | 13 enddo
|
|---|
| 60 | do 14 j=1,ndim
|
|---|
| 61 | call rebin(ndo/xnd,nd,r,xin,xi(1,j))
|
|---|
| 62 | 14 enddo
|
|---|
| 63 | ndo=nd
|
|---|
| 64 | endif
|
|---|
| 65 | if(nprn.ge.0) write(*,200) ndim,calls,it,itmx,nprn,
|
|---|
| 66 | * ALPH,mds,nd,(j,region(j),j,region(j+ndim),j=1,ndim)
|
|---|
| 67 | endif
|
|---|
| 68 | do 28 it=1,itmx
|
|---|
| 69 | ti=0.d0
|
|---|
| 70 | tsi=0.d0
|
|---|
| 71 | do 16 j=1,ndim
|
|---|
| 72 | kg(j)=1
|
|---|
| 73 | do 15 i=1,nd
|
|---|
| 74 | d(i,j)=0.d0
|
|---|
| 75 | di(i,j)=0.d0
|
|---|
| 76 | 15 enddo
|
|---|
| 77 | 16 enddo
|
|---|
| 78 | 10 continue
|
|---|
| 79 | fb=0.d0
|
|---|
| 80 | f2b=0.d0
|
|---|
| 81 | do 19 k=1,npg
|
|---|
| 82 | wgt=xjac
|
|---|
| 83 | do 17 j=1,ndim
|
|---|
| 84 | xn=(kg(j)-ran2(idum))*dxg+1.d0
|
|---|
| 85 | ia(j)=max(min(int(xn),NDMX),1)
|
|---|
| 86 | if(ia(j).gt.1)then
|
|---|
| 87 | xo=xi(ia(j),j)-xi(ia(j)-1,j)
|
|---|
| 88 | rc=xi(ia(j)-1,j)+(xn-ia(j))*xo
|
|---|
| 89 | else
|
|---|
| 90 | xo=xi(ia(j),j)
|
|---|
| 91 | rc=(xn-ia(j))*xo
|
|---|
| 92 | endif
|
|---|
| 93 | x(j)=region(j)+rc*dx(j)
|
|---|
| 94 | wgt=wgt*xo*xnd
|
|---|
| 95 | 17 enddo
|
|---|
| 96 | f=wgt*fxn(x,wgt)
|
|---|
| 97 | f2=f*f
|
|---|
| 98 | fb=fb+f
|
|---|
| 99 | f2b=f2b+f2
|
|---|
| 100 | do 18 j=1,ndim
|
|---|
| 101 | di(ia(j),j)=di(ia(j),j)+f
|
|---|
| 102 | if(mds.ge.0) d(ia(j),j)=d(ia(j),j)+f2
|
|---|
| 103 | 18 enddo
|
|---|
| 104 | 19 enddo
|
|---|
| 105 | f2b=sqrt(f2b*npg)
|
|---|
| 106 | f2b=(f2b-fb)*(f2b+fb)
|
|---|
| 107 | if (f2b.le.0) f2b=TINY
|
|---|
| 108 | ti=ti+fb
|
|---|
| 109 | tsi=tsi+f2b
|
|---|
| 110 | if(mds.lt.0)then
|
|---|
| 111 | do 21 j=1,ndim
|
|---|
| 112 | d(ia(j),j)=d(ia(j),j)+f2b
|
|---|
| 113 | 21 enddo
|
|---|
| 114 | endif
|
|---|
| 115 | do 22 k=ndim,1,-1
|
|---|
| 116 | kg(k)=mod(kg(k),ng)+1
|
|---|
| 117 | if(kg(k).ne.1) goto 10
|
|---|
| 118 | 22 enddo
|
|---|
| 119 | tsi=tsi*dv2g
|
|---|
| 120 | wgt=1.d0/tsi
|
|---|
| 121 | si=si+dble(wgt)*dble(ti)
|
|---|
| 122 | schi=schi+dble(wgt)*dble(ti)**2
|
|---|
| 123 | swgt=swgt+dble(wgt)
|
|---|
| 124 | tgral=si/swgt
|
|---|
| 125 | chi2a=max((schi-si*tgral)/(it-.99d0),0.d0)
|
|---|
| 126 | sd=sqrt(1./swgt)
|
|---|
| 127 | tsi=sqrt(tsi)
|
|---|
| 128 | if(nprn.ge.0)then
|
|---|
| 129 | write(*,201) it,ti,tsi,tgral,sd,chi2a
|
|---|
| 130 | if(nprn.ne.0)then
|
|---|
| 131 | do 23 j=1,ndim
|
|---|
| 132 | write(*,202) j,(xi(i,j),di(i,j),
|
|---|
| 133 | * i=1+nprn/2,nd,nprn)
|
|---|
| 134 | 23 enddo
|
|---|
| 135 | endif
|
|---|
| 136 | endif
|
|---|
| 137 | do 25 j=1,ndim
|
|---|
| 138 | xo=d(1,j)
|
|---|
| 139 | xn=d(2,j)
|
|---|
| 140 | d(1,j)=(xo+xn)/2.d0
|
|---|
| 141 | dt(j)=d(1,j)
|
|---|
| 142 | do 24 i=2,nd-1
|
|---|
| 143 | rc=xo+xn
|
|---|
| 144 | xo=xn
|
|---|
| 145 | xn=d(i+1,j)
|
|---|
| 146 | d(i,j)=(rc+xn)/3.d0
|
|---|
| 147 | dt(j)=dt(j)+d(i,j)
|
|---|
| 148 | 24 enddo
|
|---|
| 149 | d(nd,j)=(xo+xn)/2.d0
|
|---|
| 150 | dt(j)=dt(j)+d(nd,j)
|
|---|
| 151 | 25 enddo
|
|---|
| 152 | do 27 j=1,ndim
|
|---|
| 153 | rc=0.d0
|
|---|
| 154 | do 26 i=1,nd
|
|---|
| 155 | if(d(i,j).lt.TINY) d(i,j)=TINY
|
|---|
| 156 | r(i)=((1.d0-d(i,j)/dt(j))/(log(dt(j))-log(d(i,j))))**ALPH
|
|---|
| 157 | rc=rc+r(i)
|
|---|
| 158 | 26 enddo
|
|---|
| 159 | call rebin(rc/xnd,nd,r,xin,xi(1,j))
|
|---|
| 160 | 27 enddo
|
|---|
| 161 | 28 enddo
|
|---|
| 162 | return
|
|---|
| 163 | 200 FORMAT(/' input parameters for vegas: ndim=',i3,' ncall=',f8.0
|
|---|
| 164 | * /28x,' it=',i5,' itmx=',i5
|
|---|
| 165 | * /28x,' nprn=',i3,' alph=',f5.2/28x,' mds=',i3,' nd=',i4
|
|---|
| 166 | * /(30x,'xl(',i2,')= ',g11.4,' xu(',i2,')= ',g11.4))
|
|---|
| 167 | 201 FORMAT(/' iteration no.',I3,': ','integral =',g14.7,'+/- ',g9.2
|
|---|
| 168 | * /' all iterations: integral =',g14.7,'+/- ',g9.2,
|
|---|
| 169 | * ' chi**2/it''n =',g9.2)
|
|---|
| 170 | 202 FORMAT(/' data for axis ',I2/' X delta i',
|
|---|
| 171 | * ' x delta i ',' xdelta i ',
|
|---|
| 172 | * /(1x,f7.5,1x,g11.4,5x,f7.5,1x,g11.4,5x,f7.5,1x,g11.4))
|
|---|
| 173 | END
|
|---|
| 174 |
|
|---|
| 175 | c...............................................................
|
|---|
| 176 | SUBROUTINE rebin(rc,nd,r,xin,xi)
|
|---|
| 177 | INTEGER nd
|
|---|
| 178 | REAL*8 rc,r(*),xi(*),xin(*)
|
|---|
| 179 | INTEGER i,k
|
|---|
| 180 | REAL*8 dr,xn,xo
|
|---|
| 181 | k=0
|
|---|
| 182 | xn=0.d0
|
|---|
| 183 | dr=0.d0
|
|---|
| 184 | do 11 i=1,nd-1
|
|---|
| 185 | 1 if(rc.gt.dr)then
|
|---|
| 186 | k=k+1
|
|---|
| 187 | dr=dr+r(k)
|
|---|
| 188 | xo=xn
|
|---|
| 189 | xn=xi(k)
|
|---|
| 190 | goto 1
|
|---|
| 191 | endif
|
|---|
| 192 | dr=dr-rc
|
|---|
| 193 | xin(i)=xn-(xn-xo)*dr/r(k)
|
|---|
| 194 | 11 enddo
|
|---|
| 195 | do 12 i=1,nd-1
|
|---|
| 196 | xi(i)=xin(i)
|
|---|
| 197 | 12 enddo
|
|---|
| 198 | xi(nd)=1.d0
|
|---|
| 199 | return
|
|---|
| 200 | END
|
|---|
| 201 | c...............................................................
|
|---|
| 202 | FUNCTION ran2(idum)
|
|---|
| 203 | INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
|
|---|
| 204 | REAL*8 ran2,AM,EPS,RNMX
|
|---|
| 205 | PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
|
|---|
| 206 | * IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,
|
|---|
| 207 | * IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2d-7,RNMX=1.d0-EPS)
|
|---|
| 208 | INTEGER idum2,j,k,iv(NTAB),iy
|
|---|
| 209 | SAVE iv,iy,idum2
|
|---|
| 210 | DATA idum2/123456789/, iv/NTAB*0/, iy/0/
|
|---|
| 211 | if (idum.le.0) then
|
|---|
| 212 | idum=max(-idum,1)
|
|---|
| 213 | idum2=idum
|
|---|
| 214 | do 11 j=NTAB+8,1,-1
|
|---|
| 215 | k=idum/IQ1
|
|---|
| 216 | idum=IA1*(idum-k*IQ1)-k*IR1
|
|---|
| 217 | if (idum.lt.0) idum=idum+IM1
|
|---|
| 218 | if (j.le.NTAB) iv(j)=idum
|
|---|
| 219 | 11 enddo
|
|---|
| 220 | iy=iv(1)
|
|---|
| 221 | endif
|
|---|
| 222 | k=idum/IQ1
|
|---|
| 223 | idum=IA1*(idum-k*IQ1)-k*IR1
|
|---|
| 224 | if (idum.lt.0) idum=idum+IM1
|
|---|
| 225 | k=idum2/IQ2
|
|---|
| 226 | idum2=IA2*(idum2-k*IQ2)-k*IR2
|
|---|
| 227 | if (idum2.lt.0) idum2=idum2+IM2
|
|---|
| 228 | j=1+iy/NDIV
|
|---|
| 229 | iy=iv(j)-idum2
|
|---|
| 230 | iv(j)=idum
|
|---|
| 231 | if(iy.lt.1) iy=iy+IMM1
|
|---|
| 232 | ran2=min(AM*iy,RNMX)
|
|---|
| 233 | return
|
|---|
| 234 | END
|
|---|
| 235 | c....................................................................
|
|---|
| 236 |
|
|---|
