
!C----------------------------------------------------------------------C
!C                                                                      C
!C       Swendsen-Wang Monte Carlo update of the Ising model            C
!C       (Swendsen and Wang, Phys. Rev. Lett. 58 (1987) 86).            C
!C                                                                      C
!C       Self-labeling (sequential local algorithm plus local           C
!C       label propagation of edges between processors) is used         C
!C       to do the parallel component labeling. See Clive Baillie       C
!C       and Paul Coddington, Concurrency: Practice and Experience      C
!C       3, 129 (1991).                                                 C
!C                                                                      C
!C       Paul Coddington, August 1992, June 1993.                       C
!C                                                                      C
!C----------------------------------------------------------------------C

	SUBROUTINE SWupdate(iterations, p)
        INCLUDE "common.h"
	INCLUDE "cluster.h"
        INCLUDE "/usr/lpp/ppe.poe/include/mpif.h"
	INTEGER iterations
	INTEGER iter,i,j,n,num_clusters
	REAL    p, rr1,rr2,random


!C ----  Number of iterations before checking for termination
!C ----  (N.B. the value is purely empirical).
	minits = 1.8 * SQRT(FLOAT(nprocs))

!C ----  Offset to make local cluster labels global.
	offset = sites * procnum

	DO iter=1,iterations

!C ----  Set up the bonds between sites.
	DO i=1,isiz
	DO j=1,jsiz
	   rr1 = random()
	  IF (( spin(i,j).EQ.spin(i+1,j) ).AND.( p.GT.rr1 )) THEN
	    ibond(i,j) = .TRUE.
	  ELSE
	    ibond(i,j) = .FALSE.
	  ENDIF
	   rr2 = random()
	  IF (( spin(i,j).EQ.spin(i,j+1) ).AND.( p.GT.rr2 )) THEN
	    jbond(i,j) = .TRUE.
	  ELSE
	    jbond(i,j) = .FALSE.
	  END IF
	END DO
	END DO


	DO n=1,sites
	   Npoint(n) = 0
	ENDDO


!C ----  Pass bond edges.
        CALL bond_edges

!C ----  Label the connected sub-clusters of sites on each processor.
	CALL cluster_label

!C ----  Choose a random spin value for each cluster.
	num_clusters = 0
	DO n=1,sites
	  IF (Npoint(n).GT.0) THEN
	    num_clusters = num_clusters + 1
	    rr1 = random()
	    IF ( rr1.LT.0.5 ) THEN
	      new_spin(n) = 1
	    ELSE
	      new_spin(n) = -1
	    END IF
	  END IF
	END DO

!C ----  Initialize the global cluster labels to be the site number.
        DO n=0, sites+1
	  Npoint(n) = n + offset
	END DO

!C ----  Get the global cluster labels by matching up the edges.
        CALL global_label

!C ----  Update the spins at every site according to their cluster label.
	DO i=1,isiz
	DO j=1,jsiz
	  spin(i,j) = new_spin(label(i,j))
	  label(i,j) = Npoint(label(i,j))
	END DO
	END DO

!C ----  Copy the edges. 
        CALL spin_edges


!C ----  Debug.
!C       CALL check_labeling
        

	END DO
	RETURN 
	END



!!C----------------------------------------------------------------------C
!C                                                                      
!C
!C   Hoshen-Kopelmann algorithm for labeling connected clusters
!C   (J. Hoshen and R. Kopelman, Phys. Rev. B 14 (1976) 3438).
!C
!C   Every site has a pointer N. If this pointer is negative, then 
!C   it points to the label of the root site of the cluster. If it
!C   is positive, then the site is the root, and N refers to the 
!C   number of sites in the cluster, and the label is the given by
!C   the original label of that site. 
!C
!C   This differs from the sequential algorithm in that periodic
!C   boundary conditions are not imposed.
!C
!C----------------------------------------------------------------------C

	SUBROUTINE cluster_label
        INCLUDE "common.h"
	INCLUDE "cluster.h"
	INTEGER cluster_number
	INTEGER i,j,l,n,m(4),current_label,neighbors,min,sum
	LOGICAL different
	INTEGER classify,ii,jj
	REAL    rand


        cluster_number = 0

	DO i=1,isiz
	DO j=1,jsiz

	  neighbors = 0

!C ------  Check if the site is bonded to sites which have already been
!C ------  done (i.e. those in the negative directions). If so, use the
!C ------  subroutine classify() to determine the proper cluster labels 
!C ------  of the neighbors. 


	  IF ( ibond(i-1,j) .AND. (i.NE.1) ) THEN
	    current_label = label(i-1,j)
	    neighbors = neighbors + 1
	    m(neighbors) = classify(current_label)
	  ENDIF

	  IF ( jbond(i,j-1) .AND. (j.NE.1) ) THEN
	    current_label = label(i,j-1)
	    neighbors = neighbors + 1
	    m(neighbors) = classify(current_label)
	  ENDIF


	  IF (neighbors.GT.0) THEN
!C --------  Work out the proper cluster label.
	    min = m(1)
	    sum = Npoint(m(1))
!C --------  Find the smallest (i.e. proper) cluster label.
	    DO n=2,neighbors
	      IF (m(n).LT.min) min = m(n)
!C ----------  If clusters merge, add the cluster sizes.
	      different = .TRUE.
		DO l=1,n-1
		  IF (m(n).EQ.m(l))  then 
		     different = .FALSE.
		  endif   
		ENDDO
		IF (different) sum = sum + Npoint(m(n))
	    ENDDO
	    label(i,j) = min
	    Npoint(min) = sum + 1
!C --------  Set pointers to the proper cluster label.
	    DO n=1,neighbors
	       if (m(n).NE.min) Npoint(m(n)) = -min
	    ENDDO
	  ELSE 
!C --------  Give the site a new cluster label.
	    cluster_number = cluster_number + 1
	    label(i,j) = cluster_number
	    Npoint(cluster_number) = 1
          ENDIF

	ENDDO
	ENDDO


!C ----  Loop over sites again to find the size of the clusters 
!C ----  and set the labels to be the proper cluster label.
	
	DO i=1,isiz
	DO j=1,jsiz

!C ------  Work out the proper cluster label.
	  label(i,j) = classify(label(i,j))

	ENDDO
	ENDDO


	RETURN 
	END



!C----------------------------------------------------------------------C
!C 
!C  Go down the tree of pointers until you hit the root, which 
!C  is returned as the proper label.
!C
!C----------------------------------------------------------------------C

	FUNCTION classify(current_label)
        INCLUDE "common.h"
	INCLUDE "cluster.h"
	INTEGER current_label
	INTEGER classify
	INTEGER r, t 

  	r = current_label
  	t = -Npoint(r)


  	IF (t.LT.0) THEN
    	  classify = r
	  GOTO 500
  	ELSE 
    	  r = t
    	  t = -Npoint(t)
    	  IF (t.LT.0) THEN
    	    classify = r
	    GOTO 500
    	  ELSE 
 600   	    r = t
            t = -Npoint(t)
            IF (t.GE.0) GOTO 600
      	  ENDIF
    	  Npoint(current_label) = -r
    	  classify = r
      	ENDIF

 500    CONTINUE
  	RETURN 
  	END

