!  Connected components in the grid.  -- for XMX by YMX array
!  File:  grcc.f90                              Ernest Sibert - September 1996
!
!  Axis convention:
!         0     1       --- y --->              YMX-1
!      +-----+-----+-----+-----+-----+- . . . -+-----+
!    0 |     |     |     |     |     |         |     |
!      +-----+-----+-----+-----+-----+- . . . -+-----+
!    1 |     |     |     |     |     |         |     |
!  |   +-----+-----+-----+-----+-----+- . . . -+-----+
!  x   |     |     |     |     |     |         |     |
!  |   +-----+-----+-----+-----+-----+- . . . -+-----+
!  v   :     :     :     :     :     :         :     :
!      +-----+-----+-----+-----+-----+- . . . -+-----+
! XMX-1|     |     |     |     |     |         |     |
!      +-----+-----+-----+-----+-----+- . . . -+-----+

!  Parallel vars for cell data.  These are in common block /GRCCPE/
!  defined in grcc.h.  Parameters are defined by grcc.hp, which is
!  included automatically with grcc.h.
!  All arrays are (0:XMX-1, 0:YMX-1)
!
!  integer:  cc    Component number
!            cw    Working component
!             x    x coord (axis 1)
!             y    y coord (axis 2)
!	     id    Unique identifying number
!
!  logical:   nb,  Northern barrier
!             eb,  Eastern barrier
!             wb,  Western barrier
!             sb;  Southern barrier


!  Construct array of random integers in the range [0..]up

      function Ranarray(up)
      INCLUDE 'grcc.h'
      integer up
      integer Ranarray(0:XMX-1, 0:YMX-1)

      real W(XMX*YMX)

      call Random_Number(W)               ! Get random numbers in [0..]1.0
    ! Reshape to XMX by YMX, scale, convert to integer
      Ranarray = int(up*RESHAPE(W, (/ XMX, YMX /) ) )
      end function Ranarray


!  Setup random test with probability of barrier = p/100.
!  Returns total number of barriers installed [nb + wb]
!  Expects 0 <= p <= 100.

      integer function setup(p)
      integer p

      INCLUDE 'grcc.h'

      INTERFACE
      function Ranarray(up)
      INCLUDE 'grcc.h'
      integer up
      integer Ranarray(0:XMX-1, 0:YMX-1)
      END INTERFACE

      integer nc, wc                         ! Barrier counts (north, west)
      integer ir(0:XMX-1, 0:YMX-1)
      integer ii

      x = SPREAD( (/ (ii, ii=0,XMX-1) /), 2, YMX)      ! x coordinate
      y = SPREAD( (/ (ii, ii=0,YMX-1) /), 1, XMX)      ! y coordinate
      id = y + YMX*x			               ! Unique identification
      nb = (x == 0)		! Northern edge
      wb = (y == 0)		! Western edge
!  Add random barriers on North and West
      ir = ranarray(100)
      where ( ir < p ) nb = .TRUE.
      ir = ranarray(100)
      where ( ir < p ) wb = .TRUE.
!  Make Eastern and Southern barriers consistent
      sb = CSHIFT(nb, DIM=1, SHIFT=1)
      eb = CSHIFT(wb, DIM=2, SHIFT=1)
!  Count barriers
      ir = 0
      where (nb) ir = 1
      where (wb) ir = ir + 1
      setup = SUM(ir)
      end function setup


!  Calculate connected components, leaving labels in cc.
!  Returns no. of cycles, not counting final test.

      integer function gccs()
      INCLUDE 'grcc.h'

      INTERFACE
	logical function gccx()
      END INTERFACE
	
      integer nc		     ! Counts number of cycles

      nc = 0
      cc = id			     ! Unique numbers
      do while ( gccx() )            ! Repeat until solved
	nc = nc + 1
      end do
      gccs = nc		             ! Result
      end function gccs


!  Perform one cycle of sweeps, test for changes.
!  Updates component numbers in cc.
!  Returns test result, .TRUE. indicates one or more changes.

      logical function gccx()
      INCLUDE 'grcc.h'
      INCLUDE 'scans.hgI'

      INTERFACE
	subroutine dcc(cp, i, j, r, c)
        INCLUDE 'grcc.hp'
        integer, dimension(0:XMX-1, 0:YMX-1) :: cp
	integer i, j, r, c
      END INTERFACE

      integer, dimension(0:XMX-1, 0:YMX-1) :: cw

      cw = cc
      cw = SegMinScan(cw, 2, wb, Inclusive, Upward)
      cw = SegMinScan(cw, 2, eb, Inclusive, Downward)
      cw = SegMinScan(cw, 1, nb, Inclusive, Upward)
      cw = SegMinScan(cw, 1, sb, Inclusive, Downward)
      gccx = ANY(cw <> cc)	 ! Any changes?
      cc = cw
      end function gccx


!  Calculate connected components, displaying labels at each iteration

      integer function gcc1()
      INCLUDE 'grcc.h'

      INTERFACE
	logical function gccx()
      END INTERFACE
	
      integer nc		     ! Counts number of cycles
      character*80 pbuff             ! To pause for user signal

      nc = 0
      cc = id			     ! Unique numbers
      print *, "Initial labels:";  call dcc(cc, 0, 0, XMX, YMX)
      print '("Type p and press return > "$)';  read *, pbuff
      do while ( gccx() )            ! Repeat until solved
	nc = nc + 1
        print *, "After iteration: ", nc;  call dcc(cc, 0, 0, XMX, YMX)
        print '("Type p and press return > "$)';  read *, pbuff
      end do
      gcc1 = nc		             ! Result
      end function gcc1


!  Draw horizontal line in grid according to b in row i,
!  starting in column j for c columns.

      subroutine bline(b, i, j, c)
      INCLUDE 'grcc.hp'
      logical b(0:XMX-1, 0:YMX-1)
      integer i, j, c

      integer k

      print '("+"$)'
      k = j
      do (c) times
	if ( b(i, k) ) then
	  print '("---+"$)'
	else
	  print '("   +"$)'
	end if
	k = k + 1
      end do
      print '(1X)'	! Finishes the line
      end subroutine bline


!  Print a portion of the current grid showing barriers.
!  Uses cp for component numbers, nb, sb, wb, eb for barriers.

      subroutine dcc(cp, i, j, r, c)
      INCLUDE 'grcc.h'
      integer i, j, r, c
      integer, dimension(0:XMX-1, 0:YMX-1) :: cp

      INTERFACE
	subroutine bline(b, i, j, c)
	INCLUDE 'grcc.hp'
	logical b(0:XMX-1, 0:YMX-1)
	integer i, j, c
      END INTERFACE

      integer mm, nn, k, ix, ci

      mm = i + r - 1	! Row limit
      nn = j + c - 1	! Column limit
      ix = i
      call bline(nb, i, j, c)
      do while ( ix <= mm )
	if ( wb(ix, j) ) then
	  print '("|"$)'
	else
	  print '(" "$)'
	end if
	do, k = j, nn
	  ci = cp(ix, k)	! Suspect compiler problem with
	  print '(I3$)', ci	! cc(ix, k) in print statement
	  if ( eb(ix, k) ) then
	    print '("|"$)'
	  else
	    print '(" "$)'
	  end if
	end do	! k
	print '(1X)'	! End line
	call bline(sb, ix, j, c)
	ix = ix + 1
      end do	! ix
      end subroutine dcc


!  Main program -- manages tests and timing interactively.
!
!  Responses (prompts for additional input as needed):
!
!   s			- setup test with probability dd [0..100]
!      <dd>
!   r                   - run current example
!   d                   - run current example, display at each iteration
!   w			- print [i..]i+r by [j..]j+c  (r by c)
!      <i> <j> <r> <c>
!   q                   - quit

      program gridcc
      INCLUDE 'grcc.h'

      INTERFACE
	integer function setup(p)
	integer p
      END INTERFACE
      INTERFACE
	integer function gccs()
      END INTERFACE
      INTERFACE
	integer function gcc1()
      END INTERFACE
      INTERFACE
	subroutine dcc(cp, i, j, r, c)
        INCLUDE 'grcc.hp'
        integer, dimension(0:XMX-1, 0:YMX-1) :: cp
	integer i, j, r, c
      END INTERFACE

      character c
      integer pr, bs, nc, rg, ii, jj, rw, cl
      integer, dimension(0:XMX-1, 0:YMX-1) :: bb

      print 21, XMX, YMX
 21   FORMAT (" Grid is ", I4, " by ", I4, ".")
      call scaninit()
      bs = setup(0)
!   Read - execute - print loop
 11   print '(" grid cc> "$)'
      read *, c
      if      ( c == 's' ) then
	print '("  probability> "$)'
	read *, pr
	bs = setup(pr)
        cc = 0                  ! So display looks right
      else if ( c == 'r' ) then
	nc = gccs()
       !  Count number of regions
	bb = 0;  WHERE ( id == cc ) bb = 1;  rg = SUM(bb)
	print 22, bs, nc, rg
 22	FORMAT (1X, I7, " barriers, ", I6, " cycles,", I6, " regions.")
      else if ( c == 'd' ) then
        nc = gcc1()
	bb = 0;  WHERE ( id == cc ) bb = 1;  rg = SUM(bb)
	print 22, bs, nc, rg
      else if ( c == 'w' ) then
	print '("  i j r c> "$)'
	read *, ii, jj, rw, cl
	call dcc(cc, ii, jj, rw, cl)
      else if ( c == 'q' ) then
	go to 99		! Exit
      else
	print '(" Unrecognized input -- ignored.")'
      end if
      go to 11		! Next command
 99   continue		! For exit
      end program gridcc
