!C----------------------------------------------------------------------C
!C                                                                      C
!C     Parallel Swendsen-Wang cluster Monte Carlo algorithm for the     C
!C     2-dimensional Ising model.                                       C
!C                                                                      C
!C     Uses the Get/Send connected component labeling algorithm         C
!C     described in the paper "New SIMD Algorithms for Cluster          C
!C     Labeling on Parallel Computers", by John Apostolakis, Paul       C
!C     Coddington and Enzo Marinari, NPAC technical report SCCS-279.    C
!C     Int. J. Mod. Phys. C 4, 749 (1993).                              C
!C                                                                      C
!C     John Apostolakis and Paul Coddington,                            C
!C     May 1992, August 1992, June 1993.                                C
!C                                                                      C
!C     Fortran 90 version by Paul Coddington, paulc@npac.syr.edu        C
!C     Version 1.0, June 1997.                                          C
!C                                                                      C
!C     This is just the HPF version with FORALL statements replaced     C
!C     with DO loops, and HPF compiler directives and a call to the     C
!C     HPF intrinsic NUMBER_OF_PROCESSORS removed. FORALL is not part   C
!C     of standard Fortran 90 but is in the Fortran 95 specification.   C
!C                                                                      C
!C----------------------------------------------------------------------C

      PROGRAM SWcluster
      USE global_variables
      INTEGER m, seed, start
      REAL time1(2), time2(2) 

!C
!C  Initialize.
!C
!C --- Read in parameters.
      CALL input(seed,start)
 
!C --- Initialize the random number generator.
      CALL HPF_RANDOM_SEED(seed)
 
!C --- Initialize the spins.
      CALL init_spins(start)
 
!C --- Initialize the measurement data.
      CALL init_measure
 
!C --- Set up the initial labels.
      CALL init_labels
 
!C --- Set up Swendsen-Wang bond probability.
      p = 1.0 - EXP(-2.0*beta)
      PRINT *, ' Swendsen-Wang probability =',p

!C --- Initialize timer. 
!C      CALL dtime(time1)

!C
!C  Thermalize the configuration.
!C
      CALL SWupdate(therm_sweeps)
!C
!C  Swendsen-Wang update, taking measurements every 'sweeps' iterations.
!C
      DO m = 1, measurements
        CALL SWupdate(sweeps)
        CALL measure
      END DO
 
!C --- Get time taken. 
!C      CALL dtime(time2)

!C
!C  Print out the results.
!C
      CALL output
!C     CALL picture

      PRINT *
!C      PRINT *, 'Time taken = ', time2(1) - time1(1)

      STOP
      END




!C----------------------------------------------------------------------C
!C                                                                      C
!C     Swendsen-Wang cluster Monte Carlo update.                        C
!C                                                                      C
!C----------------------------------------------------------------------C
 
      SUBROUTINE SWupdate(iterations)
      USE global_variables      
      INTEGER i, j, iter, iterations

      DO iter=1,iterations
!C
!C  Set up bonds with probability p between neighboring lattice sites
!C  (in both directions) with the same spin value.
!C
        CALL HPF_RANDOM_NUMBER(ranum)
        bond(0,:,:) = (CSHIFT(spin,1,1).EQ.spin).AND.(ranum.LT.p)
        CALL HPF_RANDOM_NUMBER(ranum)
        bond(1,:,:) = (CSHIFT(spin,1,2).EQ.spin).AND.(ranum.LT.p)
        
!C
!C  Set up initial labels, which are just two times the processor number, 
!C  and random spins for each label. For the Ising model there are only
!C  two possible states, so we set the last bit of the original cluster 
!C  label to be the spin state, i.e. even labels correspond to spin +1, 
!C  and odd labels to spin -1.
!C
        label = init_label
        CALL HPF_RANDOM_NUMBER(ranum)
        WHERE (ranum.LT.0.5) label = label + 1
!C
!C  Label the clusters.
!C       
        CALL cluster_label
!C
!C  Update the spins.
!C
        WHERE (MOD(label,2).EQ.0)
          spin = 1
        ELSEWHERE
          spin = -1
        END WHERE

      END DO
      RETURN
      END 



!C----------------------------------------------------------------------C
!C                                                                      C
!C  Initialize cluster labels.                                          C
!C  We set initial labels to be twice the abstract processor number,    !C
!C  so we can use the last bit of the label to represent the spin.      C
!C                                                                      C
!C----------------------------------------------------------------------C
      
      SUBROUTINE init_labels()
      USE global_variables      
      INTEGER i,j
      DO i=0,dimm1
        DO j=0,dimm1
          init_label(i,j) = 2 * (i + (j * dim))
        END DO 
      END DO 

!!F95      FORALL (i=0:dimm1, j=0:dimm1) &
!!F95     &  init_label(i,j) = 2 * (i + (j * dim))

      RETURN
      END 



!C----------------------------------------------------------------------C
!C                                                                      C
!C     Routines for doing local/get/send steps and connection           C
!C     improvement for the connected component labeling algorithm       C
!C     described in the paper "New SIMD Algorithms for Cluster          C
!C     Labeling on Parallel Computers", by John Apostolakis, Paul       C
!C     Coddington and Enzo Marinari, NPAC technical report SCCS-279.    C
!C     Int. J. Mod. Phys. C 4, 749 (1993).                              C
!C                                                                      C
!C     John Apostolakis and Paul Coddington, May 1992, August 1992.     C
!C                                                                      C
!C----------------------------------------------------------------------C

      SUBROUTINE cluster_label
      USE global_variables      
      LOGICAL finished,done
      INTEGER MAXCOUNT,count
      PARAMETER (MAXCOUNT = 200)
!C
!C  Set up the connections, which can be changed (improved),
!C  but keep the bonds fixed, to enable checking of labels.
!C
      connection = bond
!C
!C  Iterate the algorithm until the labels don't change.
!C
      count = 0
      done = .FALSE.
      DO WHILE (.NOT.done) 
        count = count + 1
        old_label = label
        CALL local_step
        CALL send_step
        CALL get_step
        finished = ALL(old_label.EQ.label)
        done = finished.OR.(count.EQ.MAXCOUNT)
      END DO
      
      acceptance = FLOAT(count)

!C
!C  This is only for debugging purposes
!C
!C     WRITE(6,*) 'labeling iterations =',count
      CALL check_labels

      RETURN
      END 




!C----------------------------------------------------------------------!C
!C
!C   Do a local iteration:
!C   Check the labels of the site's connected neighbors in all 4 
!C   directions, and set the new label to be the minimum of these 
!C   labels and the current label.
!C   Also does connection improvement, i.e. adds a connection where 
!C   the labels are the same.
!C
!C----------------------------------------------------------------------C

      SUBROUTINE local_step
        USE global_variables          
        INTEGER new_label(0:dimm1, 0:dimm1)
        INTEGER i,j, dir, dr

!HPF$ ALIGN new_label(:,:) WITH LATTICE(:,:)

!C
!C   Loop over 2 directions (X and Y).
!C
        DO dir=0,1
!C
!C   Negative direction.
!C
!C   Instead of passing the label to the neighbor as well as the
!C   connection, if there is no connection, pass a very large label 
!C   so as not to change the neighbor's label.
!C 
        new_label = label
        WHERE (.NOT.connection(dir,:,:)) 
          new_label = MAXINT 
        END WHERE
        IF(dir .EQ. 0) THEN
           new_label = CSHIFT(new_label,-1,1)
        ELSE
           new_label = CSHIFT(new_label,-1,2)
        END IF


        WHERE (new_label.LT.label) 
          label = new_label 
        END WHERE
!C 
!C   Positive direction.
!C 
        IF(dir .EQ. 0) THEN
           new_label = CSHIFT(label,1,1)
        ELSE
           new_label = CSHIFT(label,1,2)
        END IF
!C 
!C   Connection improvement. 
!C 
        WHERE (new_label.EQ.label) 
          connection(dir,:,:) = .TRUE.
        END WHERE
        WHERE ( connection(dir,:,:).AND.(new_label.LT.label) ) 
          label = new_label 
        END WHERE

      END DO    


      RETURN 
      END 


!C----------------------------------------------------------------------C
!C  
!C   Do a send step:
!C   if the label has changed, send the new label to the root
!C   of the subcluster with the old label, i.e. the site whose
!C   original label was the old label. This site accepts the
!C   minimum of these new labels (and its own label, which may
!C   also have changed) as its new label.
!C 
!C----------------------------------------------------------------------C

      SUBROUTINE send_step
      USE global_variables        
      INTEGER new_label(0:dimm1, 0:dimm1)
      INTEGER coord_0(0:dimm1, 0:dimm1)
      INTEGER coord_1(0:dimm1, 0:dimm1)
      LOGICAL changed(0:dimm1, 0:dimm1)
      LOGICAL mask(0:dimm1, 0:dimm1)
      INTEGER i,j

!HPF$ ALIGN WITH LATTICE :: coord_0, coord_1, new_label 
!HPF$ ALIGN WITH LATTICE :: changed, mask      

!C 
!C  Find the lattice co-ordinates of the site 
!C  which originally had this label.  
!C 
      coord_0 = MOD(old_label/2,dim)
      coord_1 = old_label/(2*dim)
!C
!C   Send the new label to the root, which finds the smallest of all
!C   the new labels it receives (i.e. send_with_minimum).
!C
      changed = label.NE.old_label 
      new_label = MAXINT 
      mask = .FALSE.

      DO i=0,dimm1
        DO j=0,dimm1
          mask(i,j) = (new_label(coord_0(i,j),coord_1(i,j)).GT.label(i,j))
        END DO 
      END DO 
      mask = mask.AND.changed
      DO i=0,dimm1
        DO j=0,dimm1
          IF (mask(i,j)) THEN
            new_label(coord_0(i,j),coord_1(i,j)) = label(i,j)
          END IF 
        END DO 
      END DO 

!!F95      FORALL (i=0:dimm1, j=0:dimm1) &
!!F95     &  mask(i,j) = (new_label(coord_0(i,j),coord_1(i,j)).GT.label(i,j))
!!F95      mask = mask.AND.changed
!!F95      FORALL (i=0:dimm1, j=0:dimm1, mask(i,j)) &
!!F95     &  new_label(coord_0(i,j),coord_1(i,j)) = label(i,j)

!C
!C   Update the label if it has received a smaller label.
!C
      WHERE (label.GT.new_label) 
        label = new_label 
      END WHERE

      RETURN
      END


!C----------------------------------------------------------------------C
!C
!C  Do a get step:
!C  Get the current label from the root of the subcluster,
!C  i.e. the site whose original label is the current label 
!C  of the site doing the getting.
!C 
!C----------------------------------------------------------------------C

      SUBROUTINE get_step
      USE global_variables        
      INTEGER coord_0(0:dimm1, 0:dimm1)
      INTEGER coord_1(0:dimm1, 0:dimm1)
      INTEGER i,j

!HPF$ ALIGN WITH LATTICE :: coord_0, coord_1

!C
!C  Find the lattice co-ordinates of the site 
!C  which originally had this label.
!C
      coord_0 = MOD(label/2,dim)
      coord_1 = label/(2*dim)
!C
!C  Get the current label from this site.
!C
      DO i=0,dimm1
        DO j=0,dimm1
          label(i,j) = label(coord_0(i,j),coord_1(i,j)) 
        END DO 
      END DO 

!!F95      FORALL (i=0:dimm1, j=0:dimm1) &
!!F95     &  label(i,j) = label(coord_0(i,j),coord_1(i,j)) 

      RETURN
      END





!C----------------------------------------------------------------------C
!C                                                                      C
!C  Check that the cluster labeling is consistent.                      C
!C                                                                      C
!C----------------------------------------------------------------------C

      SUBROUTINE check_labels
      USE global_variables       
      LOGICAL check(0:dimm1, 0:dimm1)
      INTEGER numbad,numbonds
      INTEGER i, j
      LOGICAL ok

!HPF$ ALIGN check(:,:) WITH LATTICE(:,:)

!C
!C  The labels of neighboring sites should be the same, or else there 
!C  should be no bond between them.
!C
      check = ((label.EQ.CSHIFT(label,1,1)).OR.(.NOT.bond(0,:,:)) &
     &    .AND.(label.EQ.CSHIFT(label,1,2)).OR.(.NOT.bond(1,:,:)))

      ok = ALL(check)
      IF (.NOT.ok) THEN
        numbad = SUM(MERGE(0,1,check))
	WRITE(6,*) 'Check of labeling fails!! Sum =', numbad
      ENDIF

      RETURN
      END 



