c
c PARALLEL RECURSIVE SPECTRAL BISECTION 
c
c This file has procedures  for Parallel Recursive Spectral Bisection
c Method (RSB) for partitioning unstructured meshes.  It is based on the
c sequential version provided by Horst Simon 
c ----- Ravi  pravi@cs.umd.edu 
c For help or information regarding the SEQUENTIAL CODE contact either:
c     Stephen T. Barnard 415-604-0171 barnard@nas.nasa.gov 
c     or   Horst D. Simon 415-604-4322 simon@nas.nasa.gov
c
c

c ***********************************partition**********

      subroutine  partition (
     .                       neqns, xadj, adj, key, rwork, rnwork,
     .                       iwork, inwork,color, mask, pfac,
     .                       npart,
     .                       msglvl)
c	 implicit none
#if defined CM5
#include "../include/machine/cm5.fort.h"
#elif defined INTEL
#include "../include/machine/intel.fort.h"
#elif defined SP1
#include "../include/machine/sp1.fort.h"
#endif

c
c**************************************************************
c  
         integer    sched, noffproc
         common /parti/sched, noffproc
C        common /rsb/ maxlan
         integer    min_neqns
         integer    i     ,  neqns , rnwork, inwork 
         integer    j     , k     , m,      npart
         integer    xadj(neqns+1), adj(*)
         integer    color(neqns), mask(*),pfac, ihalf,
     1              error
    	 integer    key(*)
         integer    new_rnwork
    	 real       rwork(rnwork) 
         integer    iwork(inwork)
    	 external   qsort, genord
    	 integer    msglvl

         integer ndiagloc,nfied
         integer diagloc,fied,ncolor


c
          integer maxlan,blksiz,output,maxrqi
              real tol,rqitol
          data maxlan,tol,blksiz,output/400,0.0001,12,6/
          data rqitol,maxrqi/0.0001,100/
          integer base, node, newcol, invmap, size
          integer teqns, twork
          integer temp, temp2

          integer level
          real    vector(1000)
          real    map(1000), pmask(1000)
    
c------------------------------------------------------------------
c

      if  ( neqns .le. 0 )  return     

      ndiagloc   = neqns
      nfied      = neqns+1
      new_rnwork = rnwork-(ndiagloc+nfied)

      diagloc = 1
      fied    = diagloc  + ndiagloc
      rnwork  = fied + nfied

      do  100  i = 1, neqns  
        mask(i)   = 1
        key(i)    = i
        color(i) = 0
 100  continue               
c                               
c     ------------------------------------------
c     ... loop over log(number of partitions)
c         in each iteration of the loop
c         partition the existing domains in half
c     -----------------------------------------
c
      do 1000  k = 1, pfac
c
        if (msglvl.ge.1) then
            if(iPARTI_mynode().eq.0) print *,"k=",k," j=",j
        end if

         call set_diag(neqns,xadj,adj,rwork(diagloc),
     .          color, mask)
c
          ncolor = 2 ** (k-1)
          call find_evec(
     .      neqns, xadj, adj, rwork(diagloc),
     .      rwork(rnwork), new_rnwork,
     .      maxrqi,maxlan,tol,blksiz,output,rqitol,msglvl,
     .      rwork(fied),mask,iwork(1), inwork,level, color,ncolor)


c
c     ----------------------------------------
c     ... fudge Fiedler vector in the presence
c         of several components
c     ----------------------------------------

          call fudge(neqns, rwork(fied))


c       ----------------------------------------
c       ... for each domain already colored
c           partition in half, i.e. cut in
c           half the domains colored with color j
c       ----------------------------------------
c      

          invmap = neqns+2*ncolor

          do i = 1, ncolor
            base = iwork(neqns+i)
            size = iwork(neqns+ncolor+i)
            teqns = size
            do m = 1, size
              rwork(diagloc+m-1) = rwork(fied+iwork(base+m)-1)
              iwork(invmap+m)    = iwork(base+m)
            end  do
	     call crib_sort(rwork,size,mask)


c            call qsort (size,rwork(diagloc),mask,error)
c             call printit(mask, rwork(diagloc), size)

            call PARTI_gisum(teqns, 1, twork)
            ihalf = teqns/2
c            print *, 'size', size, 'ihalf', ihalf
c            newcol = 2*(i-1)
c            do  m  = 1, ihalf
c              node = mask(m)
c              color(iwork(base+node)) = newcol
c            end do
c
c            newcol = 2*(i-1) + 1
c            do m = ihalf + 1, size
c              node = mask(m)
c              color(iwork(base+node)) = newcol
c            end do

             do m = 1, size
               node = iwork(invmap+m)
               if (mask(m) .le. ihalf) then
c                   color(iwork(base+node)) = 2*(i-1)
                   color(node) = 2*(i-1)
               else
c                   color(iwork(base+node)) = 2*(i-1)+1
                   color(node) = 2*(i-1)+1
               end if
             end do

           
          end do
c

 1000   continue


        return
c              
      end      

c**********************************************************************

      subroutine find_evec( 
     .  neqns,xadj,adj,diag,
     .  rwork,rnwork,
     .  maxrqi,maxlan,tol,blksiz,output,rqitol,msglvl,
     .  fied,mask,iwork, inwork,level, color, ncolor)

      implicit none

      integer sched, noffproc
      common /parti/sched, noffproc
      integer  inwork, rnwork
      integer neqns,xadj(*),adj(*),diag(*), mask(*)
      integer rwork(*),color(*),msglvl
      integer iwork(inwork)
      real    fied(*)
      integer level
      integer maxrqi,maxlan,blksiz,output
      real    tol,rqitol
c      integer ninds,nfied_c,nu0,nmark,nxadj_c,nadj_c,ndiag_c,
c     .  nvmap,nvmap_inv,
c      integer fied_c,u0,mark,xadj_c,adj_c,diag_c,vmap,vmap_inv
c      integer next_level

      integer error, ncolor

c       **** find lanczos solution for coarse mesh ***

        call genord (neqns, xadj, adj, diag, mask,
     .    maxlan, tol, blksiz, output, msglvl,
     .    rwork(1), rnwork,iwork(1), inwork, 
     .    color, ncolor, error, fied)

        if (error .ne. 0 .and. msglvl.gt.2) then
          write (*,*) ' error from Lanczos and Fiedler = ', error
        endif


      end
c     
c**********************************************************************
c

      subroutine set_diag(neqns,xadj,adj,diag,color,mask)
      implicit none
      integer neqns,xadj(*),adj(*),mask(*),color(*)
      real diag(*)
      integer sched, noffproc
      common /parti/sched, noffproc

      integer i,j,n
c
      call igather(color(neqns+1), color, sched)

      do i=1,neqns
        n = 0
        do j = xadj(i), xadj(i+1)-1
         if (color(i) .eq. color(adj(j))) then
            mask(j) = 1
            n = n+1
         else
            mask(j) = 0
         endif
        end do
        diag(i) = -max(1,n)
      enddo

      end

c**********************************************************************
c**********************************************************************
      subroutine fudge (neqns, evec)
c
      integer    neqns, i
      real       evec(neqns)
c
c     -------------------------
c     ... create fudged Fiedler
c     -------------------------
c
      do 100 i = 1, neqns
         evec(i) = 2.0 + evec(i)
100   continue
      return
      end

c**********************************************************************
c

