PROGRAM par_ctmp USE kiss_random IMPLICIT NONE INCLUDE 'mpp/shmem.fh' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! REAL(8), PARAMETER :: mstop=0.0e0, K=0.630000e0, H=-0.180000e0 INTEGER, PARAMETER :: LX=32,LY=32,NX=16,NY=16, t_size=1000000 INTEGER, PARAMETER :: NL=LX*LY,NP=NX*NY, num_esc=10 INTEGER, PARAMETER :: dt_MCSS=1, tstop=1000000/dt_MCSS+1 !enter in [MCSS] REAL(8), DIMENSION(1:10) :: p REAL(8) :: r INTEGER :: nproc,mytid, count1, count2, count3, count4, count0 INTEGER :: i,j,l,m,n,ncount,ibc,rate,ind1,ind2,ind3,ind4,denum INTEGER(4) :: kiss1, kiss2, kiss3 INTEGER, DIMENSION(4) :: NNP INTEGER, DIMENSION(NL,4) :: NNSP,NNSL INTEGER, DIMENSION(NL) :: LOOC INTEGER(4), DIMENSION(NL), SAVE :: CLSB INTEGER, SAVE :: time_MCSS REAL(8), SAVE :: ltime INTEGER(4) :: lmagn INTEGER(4), DIMENSION(t_size,2), SAVE :: tmagn INTEGER, SAVE :: total INTEGER :: lstop, num_sample REAL(8) :: avg_time, std_time, err_time, avg_cpu, avg_updt, plus REAL(8) :: avg_dt,avg_dt_esc, avg_pass, avg_fail, avg_kernel REAL(8), DIMENSION(4), SAVE :: tmp_data, tmp_sum INTEGER, DIMENSION(max(3,SHMEM_REDUCE_MIN_WRKDATA_SIZE)), SAVE :: PWRK INTEGER, DIMENSION(SHMEM_REDUCE_SYNC_SIZE), SAVE :: & PSYNC=SHMEM_SYNC_VALUE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! nproc=N$PES mytid=MY_PE() !!! write(*,*) mytid,'job has started' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! p=(/min(1.0e0,exp(-8e0*K-2e0*H)),& min(1.0e0,exp(-4e0*K-2e0*H)),& min(1.0e0,exp(-2e0*H)), & min(1.0e0,exp(+4e0*K-2e0*H)),& min(1.0e0,exp(+8e0*K-2e0*H)),& min(1.0e0,exp(+8e0*K+2e0*H)),& min(1.0e0,exp(+4e0*K+2e0*H)),& min(1.0e0,exp(+2e0*H)), & min(1.0e0,exp(-4e0*K+2e0*H)),& min(1.0e0,exp(-8e0*K+2e0*H)) /) lstop=INT(mstop*NL*NP) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(nproc /= NP) then if(mytid == 0) write(*,*) 'nproc not = NP' goto 1000 end if avg_time=0.0e0 ; std_time=0.0e0 ; avg_cpu=0.0e0 avg_dt=0.0e0 ; avg_pass=0.0e0 ; avg_fail=0.0e0 ; avg_kernel=0.0e0 avg_updt=0.0e0 ; num_sample=num_esc escapes: do n=1,num_esc call SYSTEM_CLOCK(COUNT=count0) call SYSTEM_CLOCK(COUNT=count1); kiss1=MOD(count1,2147483647)+1 call SYSTEM_CLOCK(COUNT=count2); kiss2=MOD(count2,2147483647)+1 call SYSTEM_CLOCK(COUNT=count3); kiss3=MOD(count3,2147483647)+1 kiss1=MOD(kiss1,2580_4)+1_4 !!! denum=mytid+1 !!! kiss1=1991/denum; kiss2=19896321/denum; kiss3=9832678/denum call irandset(kiss1,kiss2,kiss3) call VIRTOP call INIT call SHMEM_BARRIER_ALL() !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ncount=1 ; time_MCSS=1 ind1=0 ; ind2=0 ; ind3=0 ; ind4=0 call CTMP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call SYSTEM_CLOCK(COUNT=count4,COUNT_RATE=rate) if(total<=lstop) then plus=1.0e0*dt_MCSS*(time_MCSS-1) avg_time=avg_time+plus std_time=std_time+plus**2 avg_cpu=avg_cpu+1.0e0*(count4-count0)/(1.0e0*rate) avg_updt=avg_updt+1.0e0*plus*rate/(1.0e0*(count4-count0)) avg_dt=avg_dt+avg_dt_esc avg_pass=avg_pass+1.0e0*(ind3-ind4)/(1.0e0*(ind3+ind1-ind2)) avg_fail=avg_fail+1.0e0*ind4/(1.0e0*(ind3+ind1-ind2)) avg_kernel=avg_kernel+1.0e0*(ind1-ind2)/(1.0e0*(ind3+ind1-ind2)) else num_sample=num_sample-1 end if ! if(mytid==0) then ! write(*,*) mytid,'escape#',n,'is done, cpu_sec=',& ! 1.0e0*(count4-count0)/(1.0e0*rate),& ! ' esc_time[MCSS]=',1.0e0*dt_MCSS*(time_MCSS-1) ! write(*,*) mytid,'magn=',total/(1.0*NL*NP),'avg_dt_esc=',avg_dt_esc ! write(*,*) mytid,'pass:',1.0e0*(ind3-ind4)/(1.0e0*(ind3+ind1-ind2)),& ! 'fail:',1.0e0*ind4/(1.0e0*(ind3+ind1-ind2)),& ! 'kernel:',1.0e0*(ind1-ind2)/(1.0e0*(ind3+ind1-ind2)) ! end if call SHMEM_BARRIER_ALL() end do escapes tmp_data(1)=avg_dt tmp_data(2)=avg_pass ; tmp_data(3)=avg_fail ; tmp_data(4)=avg_kernel call SHMEM_REAL8_SUM_TO_ALL(tmp_sum,tmp_data,4,0,0,nproc,PWRK,PSYNC) if(mytid==0) then avg_time=avg_time/(1.0e0*num_sample) std_time=& (std_time-(1.0e0*num_sample)*(avg_time**2))/(1.0e0*(num_sample-1)) err_time=SQRT(std_time/(1.0e0*num_sample)) avg_cpu=avg_cpu/(1.0e0*num_sample) avg_updt=avg_updt/(1.0e0*num_sample) avg_dt=tmp_sum(1)/(1.0e0*num_sample*nproc) avg_pass=tmp_sum(2)/(1.0e0*num_sample*nproc) avg_fail=tmp_sum(3)/(1.0e0*num_sample*nproc) avg_kernel=tmp_sum(4)/(1.0e0*num_sample*nproc) write(*,*) 'cont. time, asynchro. parallel metropolis simulation' write(*,*) '**********************************************************' write(*,*) 'NX=',NX,' NY=',NY write(*,*) 'LX=',LX,' LY=',LY write(*,*) 'K=',K,' H=',H write(*,*) 'm_cutoff=',mstop write(*,*) 't_cutoff=',(tstop-1)*dt_MCSS,' MCSS' write(*,*) 'dt_MCSS=',dt_MCSS write(*,*) 'number of escape trials=',num_esc write(*,*) 'number of valid samples=',num_sample write(*,*) '**********************************************************' write(*,*) '=',avg_time,'+-',err_time,' MCSS' write(*,*) '=',avg_cpu,' cpu_sec' write(*,*) '=',avg_updt*NL write(*,*) '_=',avg_time*NL/avg_cpu write(*,*) 'avg_dt=',avg_dt,' MCS max_dt=',1.00 write(*,*) '**********************************************************' write(*,*) 'avg_pass_boundary=',avg_pass write(*,*) 'avg_fail_boundary=',avg_fail write(*,*) 'avg_kernel=',avg_kernel end if 1000 continue !!!!!!!!!!!!!! internal routines !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CONTAINS SUBROUTINE VIRTOP ! upper neighbor block: NNP(1)=mytid+NX ;if(mytid >= (NY-1)*NX) NNP(1)=mytid-(NY-1)*NX ! right neighbor block: NNP(2)=mytid+1 ; if(mod(mytid,NX) == NX-1) NNP(2)=mytid-NX+1 ! lower neighbor block: NNP(3)=mytid-NX ; if(mytid < NX) NNP(3)=mytid+(NY-1)*NX ! left neighbor block: NNP(4)=mytid-1 ; if(mod(mytid,NX) == 0) NNP(4)=mytid+NX-1 do i=1,NL NNSP(i,1)=mytid ; NNSL(i,1)=i+LX if(i >= (LY-1)*LX+1) then NNSP(i,1)=NNP(1) ; NNSL(i,1)=i-(LY-1)*LX ; endif NNSP(i,2)=mytid ; NNSL(i,2)=i+1 if(mod(i,LX)==0) then NNSP(i,2)=NNP(2) ; NNSL(i,2)=i-LX+1 ; endif NNSP(i,3)=mytid ; NNSL(i,3)=i-LX if(i <= LX) then NNSP(i,3)=NNP(3) ; NNSL(i,3)=i+(LY-1)*LX ; endif NNSP(i,4)=mytid ; NNSL(i,4)=i-1 if(mod(i,LX)==1) then NNSP(i,4)=NNP(4) ; NNSL(i,4)=i+LX-1 ; endif end do return END SUBROUTINE SUBROUTINE INIT ltime=0.0e0 ; lmagn=NL ; tmagn=0 ; total=NL*NP ; CLSB=1 do i=1,NL if( (i<=LX).OR.(i>=(LY-1)*LX+1).OR.(mod(i,LX)==0).OR. & (mod(i,LX)==1) ) then LOOC(i)=0 else LOOC(i)=1 end if end do !!!!!!! determine initial local update time !!!!!!!!!!!!!!!!!!!!!!!! r=urandxx() ltime=ltime-LOG(r) return END SUBROUTINE SUBROUTINE CTMP INTEGER :: i,ks,kps, nn,nc, nnumber REAL(8) :: dtime, next_time, meas_time REAL(8), DIMENSION(2) :: nntime INTEGER(4) :: spin,updt !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! nc=0 ; avg_dt_esc=0.0e0 meas_time=1.0e0*ncount*NL*dt_MCSS main: do !!! begin main loop !DIR$ SUPPRESS CLSB,tmagn,total,time_MCSS nc=nc+1 !!!!!!! pick site i !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! r=urandxx() ; i=INT(r*NL)+1 ks=LOOC(i) !!! write(*,*) mytid,'picked site:',i !!! test !!!!!!! update site i and its nearest neighbors ind1=ind1+1 SELECT CASE(ks) CASE(0) !!! COMMUNICATION REQUIRED IN SEGMENT BELOW !!! Check if local time <= nearest neighbor time(s)!!! ind2=ind2+1 10 nnumber=0 do nn=1,4 if(NNSP(i,nn)/=mytid) then nnumber=nnumber+1 call SHMEM_GET(nntime(nnumber),ltime,1,NNSP(i,nn)) !!!write(*,*) mytid,'time received from:',NNSP(i,nn),& !!! test !!! 'nntime=',nntime(nnumber) !!! test end if end do ind3=ind3+1 if( ltime > MINVAL(nntime(1:nnumber)) ) then ind4=ind4+1 !!! call SHMEM_BARRIER_ALL() if( (total<=lstop).OR.(time_MCSS>=tstop) ) then ! CHECK STOP exit main end if !!! if(mytid==0) PAUSE 'end od main loop' !!! test !!! call SHMEM_BARRIER_ALL() goto 10 end if !!!!!!! determine state of site !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(CLSB(i)<=5) then spin=1 else spin=-1 end if updt=5*spin kps=CLSB(i) !!! call SHMEM_BARRIER_ALL() r=urandxx() if(r <= p(kps)) then call SHMEM_INT4_ADD(CLSB(i),updt,mytid) do nn=1,4 if(LOOC(NNSL(i,nn))==1) then CLSB(NNSL(i,nn))=CLSB(NNSL(i,nn))+spin else call SHMEM_INT4_ADD(CLSB(NNSL(i,nn)),spin,NNSP(i,nn)) end if end do lmagn=lmagn-2*spin end if CASE(1) !!! !!!!!!! determine state of site !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(CLSB(i)<=5) then spin=1 else spin=-1 end if updt=5*spin kps=CLSB(i) !!! call SHMEM_BARRIER_ALL() r=urandxx() if(r <= p(kps)) then CLSB(i)=CLSB(i)+updt do nn=1,4 if(LOOC(NNSL(i,nn))==1) then CLSB(NNSL(i,nn))=CLSB(NNSL(i,nn))+spin else call SHMEM_INT4_ADD(CLSB(NNSL(i,nn)),spin,mytid) end if end do lmagn=lmagn-2*spin end if END SELECT !!!!!!! determine local time increment, dtime !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! r=urandxx() dtime=-LOG(r) next_time=ltime+dtime !!!!!!!! find next_update_time avg_dt_esc=avg_dt_esc+1.0e0*dtime !!!!!!!!!! trace avg_dt if( (total<=lstop).OR.(time_MCSS>=tstop) ) then !!! CHECK STOP exit main end if !!! write(*,*) mytid,'time=',ltime,'spin(',i,')=',spin(i),& !!! test !!! 'lmagm=',lmagn !!! test if((ltime<=meas_time).and.(meas_time