PROGRAM par_dpt_ctmp_J USE kiss_random IMPLICIT NONE INCLUDE 'mpp/shmem.fh' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! REAL(8), PARAMETER :: T_J=0.225e0, H_J=-3.0e0 REAL(8), PARAMETER :: K=1.0e0/T_J, H=H_J/T_J INTEGER, PARAMETER :: LX=64,LY=64,NX=16,NY=16 INTEGER, PARAMETER :: NL=LX*LY,NP=NX*NY, NUM_PER=10000 REAL(8), PARAMETER :: HALF_PER=20.00e0 !in units of [MCSS] !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! REAL(8), DIMENSION(10,2) :: p INTEGER, DIMENSION(10) :: DE REAL(8) :: r INTEGER :: nproc,mytid INTEGER :: i,j,l,m,n,ncount,denum, count1, count2, count3 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 REAL(8), SAVE :: ltime INTEGER :: lmagn, lenerg INTEGER :: sign_H, index_H REAL(8), DIMENSION(4), SAVE :: tmp_data, tmp_sum REAL(8), 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(1,1)=1.0/(1.0+exp(8e0*K+2e0*H)); p(1,2)=1.0/(1.0+exp(8e0*K-2e0*H)) p(2,1)=1.0/(1.0+exp(4e0*K+2e0*H)); p(2,2)=1.0/(1.0+exp(4e0*K-2e0*H)) p(3,1)=1.0/(1.0+exp(2e0*H)); p(3,2)=1.0/(1.0+exp(-2e0*H)) p(4,1)=1.0/(1.0+exp(-4e0*K+2e0*H)); p(4,2)=1.0/(1.0+exp(-4e0*K-2e0*H)) p(5,1)=1.0/(1.0+exp(-8e0*K+2e0*H)); p(5,2)=1.0/(1.0+exp(-8e0*K-2e0*H)) p(6,1)=1.0/(1.0+exp(-8e0*K-2e0*H)); p(6,2)=1.0/(1.0+exp(-8e0*K+2e0*H)) p(7,1)=1.0/(1.0+exp(-4e0*K-2e0*H)); p(7,2)=1.0/(1.0+exp(-4e0*K+2e0*H)) p(8,1)=1.0/(1.0+exp(-2e0*H)); p(8,2)=1.0/(1.0+exp(2e0*H)) p(9,1)=1.0/(1.0+exp(4e0*K-2e0*H)); p(9,2)=1.0/(1.0+exp(4e0*K+2e0*H)) p(10,1)=1.0/(1.0+exp(8e0*K-2e0*H)); p(10,2)=1.0/(1.0+exp(8e0*K+2e0*H)) DE(1)=8 DE(2)=4 DE(3)=0 DE(4)=-4 DE(5)=-8 DE(6)=-8 DE(7)=-4 DE(8)=0 DE(9)=4 DE(10)=8 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(nproc /= NP) then if(mytid == 0) write(*,*) 'nproc not = NP' goto 1000 end if if(mytid == 0) then write(*,*) '# cont. time, async. par. mp. for hyteresis' write(*,*) '# ********************************************************' write(*,*) '# NX=',NX,' NY=',NY write(*,*) '# LX=',LX,' LY=',LY write(*,*) '# K=',K,' H=',H write(*,*) '# HALF_PER=',HALF_PER,' NUM_PER=',NUM_PER write(*,*) '# ********************************************************' end if 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=1511/denum; kiss2=19896321/denum; kiss3=9832678/denum call irandset(kiss1,kiss2,kiss3) call VIRTOP call INIT call SHMEM_BARRIER_ALL() !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call CTMP_HYST !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call SHMEM_BARRIER_ALL() 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; lenerg=-2*NL ; CLSB=1 sign_H=-1; index_H=(3+sign_H)/2 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_HYST INTEGER :: i,ks,kps, nn, nnumber REAL(8) :: dtime, next_time, switch_time,meas_time REAL(8), DIMENSION(2) :: nntime INTEGER(4) :: spin,updt REAL(8) :: local_Q, local_E, local_m, norm_Q, norm_m REAL(8) :: old_time, delta_time, next_half_time INTEGER :: COUNT_HALF_PER, STOP_HALF_PER !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! next_half_time=HALF_PER*NL STOP_HALF_PER=2*NUM_PER COUNT_HALF_PER=0 norm_Q=1.0*HALF_PER*NL*NL*nproc norm_m=1.0*NL*nproc local_Q=0.0; local_E=0.0; old_time=0.0 tmp_data=0.0; tmp_sum=0.0 main: do !!! begin main loop !DIR$ SUPPRESS CLSB !!!!!!! pick site i !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! r=urandxx() ; i=INT(r*NL)+1 ks=LOOC(i) !!!!!!! update site i and its nearest neighbors SELECT CASE(ks) CASE(0) !!! COMMUNICATION REQUIRED IN SEGMENT BELOW !!! Check if local time <= nearest neighbor time(s)!!! 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)) end if end do if( ltime > MINVAL(nntime(1:nnumber)) ) then !!! call SHMEM_BARRIER_ALL() !!! 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,index_H)) 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! delta_time=ltime-old_time local_Q=local_Q+lmagn*delta_time local_E=local_E+lenerg*delta_time old_time=ltime !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! lmagn=lmagn-2*spin lenerg=lenerg+DE(kps) end if CASE(1) !!! NO COMMUNICATION WILL BE REQUIRED !!!!!!! 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,index_H)) 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! delta_time=ltime-old_time local_Q=local_Q+lmagn*delta_time local_E=local_E+lenerg*delta_time old_time=ltime !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! lmagn=lmagn-2*spin lenerg=lenerg+DE(kps) end if END SELECT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! r=urandxx() dtime=-LOG(r) next_time=ltime+dtime !!!!!!!! determine next_update_time call SHMEM_QUIET ltime=next_time if(next_time>next_half_time) then delta_time=next_half_time-old_time local_Q=local_Q+lmagn*delta_time local_E=local_E+lenerg*delta_time tmp_data(1)=local_Q tmp_data(2)=local_E tmp_data(3)=1.0*lmagn call SHMEM_BARRIER_ALL() !!!!!!!! half-period synchronization call SHMEM_REAL8_SUM_TO_ALL(tmp_sum,tmp_data,4,0,0,nproc,PWRK,PSYNC) COUNT_HALF_PER=COUNT_HALF_PER+1 if(mytid==0) then tmp_sum(1)=tmp_sum(1)/(norm_Q) tmp_sum(2)=tmp_sum(2)/(norm_Q) tmp_sum(3)=tmp_sum(3)/(norm_m) 20 format(I6,3E23.15) write(*,20) COUNT_HALF_PER, tmp_sum(1),tmp_sum(2),tmp_sum(3) end if if(COUNT_HALF_PER >= STOP_HALF_PER) exit main call SHMEM_BARRIER_ALL() !!!!!!!! half-period synchronization local_Q=0.0; local_E=0.0 sign_H=-sign_H; index_H=(3+sign_H)/2 old_time=next_half_time next_half_time=HALF_PER*NL*(COUNT_HALF_PER+1) end if !!! call SHMEM_BARRIER_ALL() end do main !!! end of main loop return END SUBROUTINE END PROGRAM !!!!!!!! end of program !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!