program parallel
parameter(nmu=20,nsig=20,nd=11,ng=10,nc=16)
parameter (nsim=2000)
integer k1,k2,k3,k4,k5,i,j,m,n
real xmu,sigma,d,g,c
data xmu1/.02/,xmu2/.4/
data sig1/.2/,sig2/4./
data d1/1./,d2/2./
data g1/.1/,g2/1./
data c1/1./,c2/4./
data iseed/-9/
open(unit=10,file='p10')
open(unit=20,file='S1T0p.txt')
r=ran2(iseed)
read(20,*) n
read(20,*) m
do k1=1,nmu
xmu=xmu1+(xmu2-xmu1)*(k1-1)/float(nmu-1)
do k2=1,nsig
sigma=sig1+(sig2-sig1)*(k2-1)/float(nsig-1)
do k3=1,nd
d=d1+(d2-d1)*(k3-1)/float(nd-1)
do k4=1,ng
g=g1+(g2-g1)*(k4-1)/float(ng-1)
do k5=1,nc
c=c1+(c2-c1)*(k5-1)/float(nc-1)
call varpar(nsim,n,m,d,g,c,sigma,xmu)
enddo
enddo
enddo
enddo
enddo
end
subroutine varpar(nsim,n,m,d,g,c,sigma,xmu)
integer nsim,n,m,it,icorr,ntime
real d,g,c,sigma,xmu,totcor,totstp,stepsq,pcorr
real perr,avst,errvar,stvar
totcor=0
totstp=0
totsq=0
stepsq=0
do it=1,nsim
call trpar(g,c,sigma,xmu,d,n,m,icorr,ntime)
totcor=totcor+icorr
totstp=totstp+ntime
stepsq=stepsq+ntime**2
enddo
pcorr=totcor/nsim
perr=1.-pcorr
avst=totstp/nsim
errvar=pcorr*perr/nsim
stvar=(stepsq-avst**2/nsim)/float(nsim-1)
write(10,*) perr,avst,stvar
return
end
subroutine trpar(g,c,sigma,xmu,d,n,m,icorr,ntime)
parameter(tb=20)
integer typetr
integer n,m,icorr,ntime,ndist
real g,c,sigma,xmu,d,distb0,db,tbr
typetr=0
ntime=0
ndist=n-m
distb0=-d*tb
db=rlxb(n,c,distb0)
tbr=rlxt(n,c,tb)
call walker(n,m,xmu,sigma,db,tb,tbr,g,ntime,typetr)
if(m .gt. 0) then
if(typetr .eq. 1) then
icorr=1
else
icorr=0
endif
else
if(typetr .eq. 0) then
icorr=1
else
icorr=0
endif
endif
return
end
function rlxb(nel,c,distb0)
integer nel
real c,distb0
if(nel .eq. 4) rlxb=distb0/(c*c)
if(nel .eq. 2) rlxb=distb0/c
if(nel .eq. 1) rlxb=distb0
return
end
function rlxt(nel,c,targb)
integer nel
real c,targb
if(nel .eq. 4) rlxt=targb/(c*c)
if(nel .eq. 2) rlxt=targb/c
if(nel .eq. 1) rlxt=targb
return
end
subroutine walker (nel,ntarg,xmu,sig,dbrlx,tb,tbr,g,n,idec)
parameter(nelmax=4)
integer dstvec(nelmax),tv(nelmax),markd,markt
integer nel,ntarg,n,idec,ndist,i,k
real wvec(nelmax)
real xmu,sig,dbrlx,tb,tbr,g,typeel
dimension xmean(nelmax)
n=0
ndist=nel-ntarg
do i=1,nel
dstvec(i)=3
tv(i)=3
wvec(i)=0.
enddo
do i=1,ntarg
xmean(i)=xmu
enddo
do i=1,ndist
xmean(i+ntarg)=-xmu
enddo
10 continue
do i=1,nel
dstvec(i)=3
tv(i)=3
enddo
n=n+1
do i=1,nel
call subwlk(xmean(i),sig,dbrlx,tb,tbr,g,wvec(i),nel,typeel)
if (typeel .eq. 1) then
idec=1
return
endif
if (typeel .eq. 2) then
tv(i)=0
markt=0
do k=1,nel
markt=markt+tv(k)
enddo
if (markt .eq. 0) then
idec=1
return
endif
endif
if (typeel .eq. 0) then
dstvec(i)=0
markd=0
do k=1,nel
markd=markd+dstvec(k)
enddo
if (markd .eq. 0) then
idec=0
return
endif
endif
enddo
go to 10
end
subroutine subwlk(xm,sig,dbrlx,tb,tbr,g,w,nel,typeel)
integer nel
real xgam
real xm,sig,dbrlx,tb,tbr,g,w,typeel,step
iran=1
xgam=0
step=gasdev(iran)*sig+xm
xgam=(-1.)*g
step=step*(float(nel) ** xgam)
w=w+step
if(w .ge. tb) then
typeel=1
return
endif
if(w .ge. tbr) then
typeel=2
return
endif
if(w .le. dbrlx) then
typeel=0
return
endif
typeel=3
return
end
FUNCTION gasdev(idum)
INTEGER idum
REAL gasdev
CU USES ran1
INTEGER iset
REAL fac,gset,rsq,v1,v2,ran1
SAVE iset,gset
DATA iset/0/
if (iset.eq.0) then
1 v1=2.*ran2(idum)-1.
v2=2.*ran2(idum)-1.
rsq=v1**2+v2**2
if(rsq.ge.1..or.rsq.eq.0.)goto 1
fac=sqrt(-2.*alog(rsq)/rsq)
gset=v1*fac
gasdev=v2*fac
iset=1
else
gasdev=gset
iset=0
endif
return
END
FUNCTION ran2(idum)
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
REAL ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1)
PARAMETER (IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211)
PARAMETER (IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7)
PARAMETER (RNMX=1.-EPS)
INTEGER idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if (idum.le.0) then
idum=max(-idum,1)
idum2=idum
do 11 j=NTAB+8,1,-1
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.lt.0) idum=idum+IM1
if (j.le.NTAB) iv(j)=idum
11 continue
iy=iv(1)
endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.lt.0) idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if (idum2.lt.0) idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1)iy=iy+IMM1
ran2=min(AM*iy,RNMX)
return
END
