C
C  Swendsen-Wang Monte Carlo update of the Ising model
C  (Swendsen and Wang, Phys. Rev. Lett. 58 (1987) 86).
C

	SUBROUTINE SWupdate(iterations, p)
	INCLUDE "common.h"
	INTEGER iterations
	REAL p
	INTEGER iter,i,j,n,num_clusters
	INTEGER label(MAX_SIDE,MAX_SIDE)
	INTEGER new_spin(MAX_SITES)
	LOGICAL ibond(MAX_SIDE,MAX_SIDE), jbond(MAX_SIDE,MAX_SIDE)
	INTEGER Npoint(MAX_SITES)
	REAL ranmar

	DO iter=1,iterations

C ----  Set up the bonds between sites.
	DO i=1,side
	DO j=1,side
	  IF (( spin(i,j).EQ.spin(pos(i),j) ).AND.( p.GT.ranmar() )) THEN
	      ibond(i,j) = .TRUE.
	  ELSE
	    ibond(i,j) = .FALSE.
	  ENDIF
	  IF (( spin(i,j).EQ.spin(i,pos(j)) ).AND.( p.GT.ranmar() )) THEN
	      jbond(i,j) = .TRUE.
	  ELSE
	    jbond(i,j) = .FALSE.
	  ENDIF
	ENDDO
	ENDDO

C ----  Label the connected clusters of sites.
	CALL cluster_label(ibond,jbond,label,Npoint)

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
	    IF ( ranmar().LT.0.5 ) THEN
	      new_spin(n) = 1
	    ELSE
	      new_spin(n) = -1
	    ENDIF
	  ENDIF
	ENDDO

C ----  Update the spins at every site according to their cluster label.
	DO i=1,side
	DO j=1,side
	  spin(i,j) = new_spin( label(i,j) )
	ENDDO
	ENDDO

	accept = num_clusters
        
	ENDDO
	RETURN 
	END


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

	SUBROUTINE cluster_label(ibond,jbond,label,Npoint)
	INCLUDE "common.h"
	LOGICAL ibond(MAX_SIDE,MAX_SIDE), jbond(MAX_SIDE,MAX_SIDE)
	INTEGER label(MAX_SIDE,MAX_SIDE)
	INTEGER cluster_number
	INTEGER i,j,l,n,m(4),current_label,neighbors,min,sum
	LOGICAL different
	INTEGER Npoint(MAX_SITES)
	INTEGER classify
	REAL ranmar

        cluster_number = 0
	DO i=1,side
	DO j=1,side

	  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(neg(i),j) .AND. (i.NE.1) ) THEN
	    current_label = label(neg(i),j)
	    neighbors = neighbors + 1
	    m(neighbors) = classify(Npoint,current_label)
	  ENDIF
	  IF ( jbond(i,neg(j)) .AND. (j.NE.1) ) THEN
	    current_label = label(i,neg(j))
	    neighbors = neighbors + 1
	    m(neighbors) = classify(Npoint,current_label)
	  ENDIF

C ------  Periodic boundary conditions.
	  IF ( ibond(i,j) .AND. (i.EQ.side) ) THEN
	    current_label = label(1,j)
	    neighbors = neighbors + 1
	    m(neighbors) = classify(Npoint,current_label)
	  ENDIF
	  IF ( jbond(i,j) .AND. (j.EQ.side) ) THEN
	    current_label = label(i,1)
	    neighbors = neighbors + 1
	    m(neighbors) = classify(Npoint,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)) different = .FALSE.
              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,side
	DO j=1,side

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

	ENDDO
	ENDDO

C	 CALL check_labeling(ibond,jbond,label,Npoint, cluster_number)

	RETURN 
	END


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

	FUNCTION classify(Npoint,current_label)
	INCLUDE "common.h"
	INTEGER Npoint(MAX_SITES)
	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


C---------------------------------------------------------------------------

C
C   Check that the labeling is consistent.
C

	SUBROUTINE check_labeling(ibond,jbond,label,Npoint, cluster_number)
	INCLUDE "common.h"
	LOGICAL ibond(MAX_SIDE,MAX_SIDE), jbond(MAX_SIDE,MAX_SIDE)
	INTEGER label(MAX_SIDE,MAX_SIDE)
	INTEGER Npoint(MAX_SITES)
	INTEGER cluster_number
	INTEGER i,j,n

	DO i=1,side
	DO j=1,side
	  IF ((label(i,j).LT.1).OR.(label(i,j).GT.cluster_number)) THEN
	    WRITE(6,*) 'label',i,j,' = ',label(i,j)
	  ENDIF
	  IF (( ibond(i,j) ).AND.( label(i,j).NE.label(pos(i),j) )) THEN
	    WRITE(6,*) 'OOPS! labeling is incorrect!'
	  ENDIF
	  IF (( jbond(i,j) ).AND.( label(i,j).NE.label(i,pos(j)) )) THEN
	    WRITE(6,*) 'OOPS! labeling is incorrect!'
	  ENDIF
	ENDDO
	ENDDO
	RETURN 
	END


