#include <stdio.h>

/* ---------------------------------------------------------------------- */
/*    THE POINTS ARE DIVIDED INTO CLUSTERS ACCORDING TO FROZEN BONDS      */
/*    IN A 2D SQUARE LATTICE WITH NEAREST NEIGHBOR CONNECTTIVITY          */
/*                                                                        */
/*    N spins, each spin has up to 4 forzen bonds                         */
/*                                                                        */
/* input                                                                  */
/*     - int N: number of spins                                           */
/*     - int Bond[i][k]:  i=0..N-1, k=0..3                                */
/*           index of spins connected by frozen bonds to spin i           */
/*           negative index means no bond                                 */            
/*                                                                        */
/*       example:                                                         */
/*           1 9 10 -1  - spin 0 is connected to spins 1,9,10             */
/*           0 2 -1 -1  - spin 1 is connected to spins 0,2                */
/*       note that Bonds matrix can state each bond either once or twice  */
/*                                                                        */
/* returns                                                                */
/*     - the number of clusters                                           */
/*     - int Block[i] is the index of the cluster to which spin i belongs */
/* ---------------------------------------------------------------------- */
int  Coarsening(int N, int **Bond, int *Block)
{
int  *Stack;
int  ns,i,k,i0;
int  nblock = -1;    /* number of clusters (blocks) generated */

   /* allocate memory for stack */
   if ( ( Stack = (int*)malloc(N*sizeof(int)) ) == NULL ) {
     fprintf(stderr,"\nCan not allocate memory\n");
     return(0);
   }

   for(i = 0; i < N; i++) Block[i] = -1;   /* init blocks to -1 */

   for(i = 0; i < N; i++)
      if (Block[i] == -1){        /* point i was not labeled yet */
 	ns = 0;                   /* set stack index to 0 */
	nblock ++;                
	Block[i] = nblock;
	Stack[ns] = i;            /* insert point i to stack */
	while(ns >= 0) {          /* while stack is not empty */
	  i0 = Stack[ns];
	  ns--;
	  for(k = 0; k < 4; k++){
	    if ( ( Bond[i0][k] >= 0 ) && ( Bond[i0][k] < N ) && 
		 ( Block[Bond[i0][k]] == -1 ) ) {  /* found new neighbor */
	      Block[ Bond[i0][k] ] = nblock;
              ns++;
	      Stack[ns] = Bond[i0][k];
            }
          }
	}
      }

   nblock++;

   free(Stack);
   return(nblock);
}