c
c PARALLEL RECURSIVE SPECTRAL BISECTION
c
c This file has programs 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 lanczos.f
c
      subroutine genord (neqns, xadj, adj, diag, mask,
     .    maxlan, tol, blksiz, output, 
     .    msglvl, work, rnwork, iwork, inwork, color, 
     .    ncolor, error, fied)
c
c*********************************************************************
c
c     input
c       neqns      -  number of equations
c       xadj, adj  -  adjacency structure in Boeing Harwell format
c       diag       -  diagonal of Laplacian matrix
c       maxlan     -  maximum number of Lanczos steps
c       tol        -  tolerance: a small tol results in larger
c                     execution time, but better ordering
c       blksiz     -  block size for Lanczos algorithm
c                     a smaller block size results in higher
c                     execution time
c       output     -  output unit number
c       msglvl     -  message level
c       
c     workspace
c       work       - workspace
c       rnwork      - dimension of work, minimum required
c                    is 4* neqns + 4* maxlan + maxlan**2,
c                    where maxlan = 300, currently 
c         integer workspace
c       iwork       - workspace
c       inwork      - dimension of work, minimum required
c
c     output    
c       lanstp     - number of Lanczos steps taken
c       error      - error flag
c                  = 0    normal return
c                  = 1    no convergence in Lanczos after maxlan steps
c                  = -1 ... -neqns  EISPACK error
c                  = -1-neqns   insufficient storage 
c       fied       = Fiedler vector (fudged for several components if
c                        needed)
c
c***********************************************************************
c
c      implicit     none
c
c     --------------
c     ... parameters
c     --------------
c

#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

      integer sched, noffproc
      common /parti/sched, noffproc

      integer      neqns, rnwork, inwork, blksiz, output, msglvl, 
     1             maxlan, lanstp, error
      integer      xadj (*), adj  (*), mask(*)
      real         diag (*)
      real         tol
      real         work (rnwork)
      integer      iwork(inwork), color(*), ncolor
      real         fied (neqns)
c
c     -------------------
c     ... local variables
c     -------------------
c
      integer      i,r, u, alf, bet, betc, wkneed, z  ,j,junk ,
     .             icomp, comptr, ncomp
      integer      size, gmask, indarray, indptr, count(512),eval
      integer      mynumber, nproc, mynode, numnodes, lsize, csize
      external     opk, lanczo, icopy, fncomp, fudge
c
c     -------------------------------------------------------------
c
      error  = 0
      ncomp  = 1
c
c     ----------------------
c     ... storage allocation
c     ----------------------
c
       icomp  = 1
       comptr = icomp  + neqns
c      mask   = comptr + neqns
      u      = comptr   + neqns
      r      = u + neqns 
      alf    = r + neqns + noffproc
      bet    = alf + ncolor*maxlan
      betc   = bet + ncolor*maxlan
      z      = betc + ncolor*maxlan
      eval   = z + ncolor*maxlan
      wkneed = eval + ncolor
c
      if (wkneed .gt. rnwork) then
	  if (msglvl .gt. 0) write (output, 8901) wkneed, rnwork
	  error = -1 - neqns
	  go to 9999
      endif

      nproc     = iPARTI_numnodes()
      mynumber  = iPARTI_mynode()

      indarray  = 1
      indptr    = indarray + neqns
      size      = indptr + ncolor 
      gmask     = size + ncolor
      csize    = gmask + xadj(neqns+1)
      lsize    = csize + ncolor
      wkneed   = lsize + nproc

      if (wkneed .gt. inwork) then
	  if (msglvl .gt. 0) write (output, 8901) wkneed, inwork
	  error = -1 - neqns
	  go to 9999
      endif
      
c
c     ---------------------------
c     ... call Lanczos subroutine
c     ---------------------------
c

      do i = 1, ncolor
        iwork(size+i-1) = 0
        count(i) = 0
      end do

      do i = 1, neqns
        iwork(size+(color(i))) = iwork(size+(color(i))) + 1
      end do

      do i = 1, ncolor
       call PARTI_gcol(iwork(size+i-1), 4, iwork(lsize),4*nproc,junk)
       iwork(csize+i-1) = 0
       do j = 0, mynumber-1
         iwork(csize+i-1) = iwork(csize+i-1) + iwork(lsize+j)
       end do
      end do

      work(indptr) = 0
      do i = 1, ncolor-1
        iwork(indptr+i) = iwork(indptr+i-1) + iwork(size+i-1)
      end do

      do i = 1, neqns
        iwork(indarray+iwork(indptr+color(i))+count(color(i)+1)) = i
        count(color(i)+1) = count(color(i)+1) + 1
      end do

      call lanczo (neqns, xadj, adj, 
     1             ncomp, work(icomp), 
     2             work(comptr), diag, 
     3             work(eval), fied, work(u), work(r),
     4             work(alf), work(bet), work(betc),
     5             work(z), maxlan, tol, blksiz, lanstp, error, 
     6             msglvl, output,
     7             iwork(indptr), iwork(indarray), iwork(size),
     8             mask, iwork(csize), ncolor)

c
      if (error .ne. 0)  then
	  if (msglvl .gt. 2) write (output, 8904) error
	  if (error  .lt. 0) go to 9999
      endif
c
c     --------------------
c     ... print eigenvalue
c     --------------------
c
      if (iPARTI_mynode() .eq. 0) then
       do i = 1, ncolor
       if (msglvl .gt. 1) write(output, 8905)  work(eval+i-1)
       end do
      end if
c
 9999 continue
      return
c    
c     -----------
c     ... formats
c     -----------
c
8901  format ( /1x, ' ***** insufficient workspace ******** ',
     1         /1x, ' workspace needed   =  ', i10,
     2         /1x, ' workspace provided =  ', i10  )
c
8902  format ( /1x, ' maximum degree = ', f15.1 )
c
8904  format ( /1x, ' error in Lanczos, error = ', i10 )
c
8905  format ( /1x, ' Fiedler  Value = ', 3x, 1pe22.14, / ) 
c
8906  format ( 5e16.8 )
c
      end
c
c     --------------------------
c
c     -------------------------
c
      subroutine givens (a, b, n, eig, evec,ncolor,maxlan)
      implicit none
c
c     -----------------------------------------------
c     Givens recurrence for computing the eigenvector
c     of a tridiagonal matrix, when the eigenvalue
c     is known
c     ------------------------------------------------
c
      integer   n, ncolor, maxlan
      real      eig
      real      a(*), b(*),evec(*)
c
      integer   i
      real      sum, bd
c
c
      evec(n)   = 1.0
      evec(n-1) = (eig - a(n)) * evec(n) / b(n)
      sum       = 1.0 + evec(n-1)*evec(n-1)
      do 10 i = n-2, 1, -1
         evec(i) = ((eig - a(i+1)) * evec(i+1) - b(i+2) * evec(i+2)) /
     1               b(i+1)
         sum     = sum + evec(i)*evec(i)
10    continue
      bd = 1.0 / sqrt(sum)
      do 11 i = 1, n
         evec(i) = evec(i) * bd
   11 continue

      return
      end
c
c     ----------------------
c
c     ---------------------
c
      subroutine icopy (n, iarray, inx, jarray, iny)
      integer    n, inx, iny
      integer    iarray (n), jarray(n)
      if (inx .eq. 1 .and. iny .eq. 1) then
          do 100 i = 1, n
             jarray(i) = iarray(i)
  100     continue
      else
          do 200 i = 1, n
             jarray ( 1 + (i-1)*iny ) = iarray ( 1 + (i-1)*inx )
  200     continue
      endif
      return
      end
c
c     ---------------------
c
c     ---------------------
c
      subroutine lanczo ( neqns, xadj, 
     1                   adj, ncomp, icomp, 
     2                   comptr, diag,
     3                   eval, evec, u, r, 
     4                    alf, bet, betc,
     5                   z, maxlan, eps, blksiz, lanstp,ierr, 
     6                   msglvl, output,
     7                  indptr, indarray, size, 
     8                  mask, csize, ncolor)
c
c     +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c     special version of the Lanczos method for fiedler vector
c     modification on Jan 19, 1991 added PSS method for largest
c     eigenvalue
c     +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c 
c      implicit  none
c
#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

      integer NP
      parameter (NP = 128)

      integer neqns, ncolor, indptr(ncolor), indarray(neqns),
     1        size(ncolor), mask(*)
      integer ncomp, maxlan, lanstp, blksiz, 
     1        ierr , msglvl, output
      integer xadj(*), adj(*), icomp(*), comptr(*), csize(ncolor)
      real    eval(ncolor), diag(neqns), evec(neqns), u(neqns), 
     1        r(*), alf(maxlan,ncolor), bet(maxlan,ncolor), 
     2        betc(maxlan,ncolor), z(maxlan,ncolor)
    
      integer sched, noffproc
      common /parti/sched, noffproc

c
c-- Indirection Array 
c

c
c     -------------------
c     ... local variables
c     -------------------
c
      integer   i, ic(NP), iseed, j, jmax, iflag, k
c 

      real      alfa(NP), beto(NP), be2o(NP), disc(NP), eps,  
     1          res(NP), maceps, rr(NP), s(NP), 
     2          dif(NP), du(NP), dl(NP), xc(NP), xu(NP), xl(NP), 
     3          bd(NP), eig(NP), xa(NP), w(NP), xb(NP), da(NP), db(NP), 
     4          sumi(NP), tiny, tol(NP)
c
      integer   mynode, teqns, iwork, dflag(NP)
      real      dot, nrm2
      external  axpy, scal, copy, dot, nrm2, swap,
     1          opk, ortho, givens, pivot
c
c______________________________________________________________
c


      ierr   = 0
      iflag  = 0
      maceps = 1.0e-12
      tiny   = maceps**3
c
c     ----------------------
c     ... initialize Lanczos
c     ----------------------
c
      iseed = 1
      

      do j = 1, ncolor
        dflag(j) = 0
        bet(1,j)  = 0.0e0
        betc(1,j) = 0.0e0
        teqns = size(j)
        call PARTI_gisum(teqns, 1, iwork)
        do i = 1, size(j)
         r(indarray(indptr(j)+i)) = 0.0e0
         u(indarray(indptr(j)+i))=0.498e0*float(-teqns+2*(csize(j)+i))
        enddo
      end do


      do j = 1, noffproc
         r(neqns+j) = 0.0
      end  do
      call ortho(neqns, u, indarray, indptr, size, ncolor)



      do j = 1, ncolor
        res(j) = nrm2(size(j), u, indarray(indptr(j)+1))
      end do


      do j = 1, ncolor
        res(j) = 1.0e0 / res(j)
        call scal (size(j), res(j), u, indarray(indptr(j)+1))
      end do
        call copy (neqns, u, evec)


c     ----------------
c     ... Lanczos loop
c     ----------------
      lanstp = 0
      do i = 1, ncolor
        sumi(i)   = 0.0e0
      end do

      do 5000 j = 1, maxlan-1

        do i = 1, ncolor
	 call scal(size(i), -bet(j,i), r, indarray(indptr(i)+1))
        end do

	 call opk(xadj, adj, diag, u, r, neqns,mask)

        do i = 1, ncolor
 	 alf(j,i) = dot (size(i), u, r, indarray(indptr(i)+1))
 	 call axpy (size(i), -alf(j,i), u, r,indarray(indptr(i)+1))
        end do

c
c        --------------------------------------
c        ... orthogonalize against 1, 1, ...
c        --------------------------------------
c
	 call ortho(neqns, r, indarray, indptr, size, ncolor)
c
         do i = 1, ncolor
	   bet(j+1, i) = nrm2 (size(i), r, indarray(indptr(i)+1))
	   if (bet(j+1,i) .le. maceps) iflag = 1
c
	 call scal (size(i), 1.0e0/bet(j+1,i), r, indarray(indptr(i)+1))
         call swap (size(i), u, r, indarray(indptr(i)+1))
         end do


c
c        ---------------------------------------
c        ... end of Lanczos recurrence
c            u contains q sub j+1
c            r contains q sub j
c            alf(j) and bet(j) have been computed
c            now initialize for eigen anlysis
c        ----------------------------------------
c
         do i = 1, ncolor
           alfa(i)    = alf(j,i)
           beto(i)    = bet(j,i)
	   be2o(i)    = beto(i) * beto(i)
	   betc(j,i)    = be2o(i)
           if(iPARTI_mynode().eq.0) then
           end if
           tol(i) = 0.5 * eps * abs(eig(i))
         end do
           if (j .eq. 1) then
             goto 5000
           end if
       if (j .eq. 2) then 
c
c           ------------------------
c           ... initialize xu and xl
c           ------------------------
c
         do i = 1, ncolor

            dif(i) = (alf(1,i)-alfa(i))*0.5
            xc(i)  = alfa(i) + dif(i) + sqrt(dif(i)**2+be2o(i))
            xu(i)  = xc(i) + eps * abs(xc(i))
            du(i)  = alfa(i) - xu(i) - be2o(i)/(alf(1,i)-xu(i))
            xl(i)  = 2.0*xc(i) - xu(i)
	    eig(i) = xc(i)
         end do
         else if (j .gt. 2) then
c
c           -----------------------------------
c           ... normal case for j > 2, use
c               rational function approximation
c           -----------------------------------
c
         

         do i = 1, ncolor
            du(i)  = alfa(i) - xu(i) - be2o(i)/du(i)
            ic(i)  = 0
            if (du(i).le.0.0) dflag(i) = 1
         end do

         do i = 1, ncolor
          if(dflag(i) .ne. 1) then
            xl(i)  = xu(i)
            dl(i)  = du(i)
            xa(i)  = xl(i) + bd(i)
            call pivot(alf(1,i), betc(2,i), j, xa(i), da(i))
            ic(i)  = ic(i) + 1
            eig(i) = xc(i)
            xb(i)  = xl(i) + 0.5*(xa(i) - xl(i))
            xu(i)  = beto(i) + amax1(xl(i), alfa(i))
           end if
         end do


c
c           ------------------------
c           ... iterate for the zero
c           ------------------------
c
         do 100 k = 1, ncolor
           if (dflag(k) .ne. 1) then  
            do 600 i = 1, 40

               call pivot (alf(1,k), betc(2,k), j, xb(k), db(k))
               ic(k) = ic(k) + 1
               if (db(k).gt.0.0) then
		  xl(k) = xb(k)
	       else 
                  xu(k) = xb(k)
                  du(k) = db(k)
               endif
               w(k) = xb(k)-xa(k)
	       if ( w(k) .le. 0.0e0 ) then
		 xc(k) = xb(k)
               else
		s(k)   = -w(k)-(db(k)*(xb(k)-eig(k))-da(k)*
     1                              (xa(k)-eig(k)))/w(k)
		rr(k)  = db(k) * (xb(k)-eig(k))
		disc(k)= s(k) * s(k) - 4.0e0 * rr(k)
		if (disc(k) .le. tiny) then
		  xc(k) = xb(k)
                else
                 xc(k) = xb(k)+2.0e0*rr(k)/
     1                        (s(k)+sign(sqrt(disc(k)),s(k)))
		endif
               endif
               if (xc(k).lt.(xl(k)+.5*tol(k)).or.
     1               xc(k).gt.(xu(k)-.5*tol(k))) then
                   xc(k) = xl(k)+0.5*(xu(k)-xl(k))
               endif
               xa(k) = xb(k)
               da(k) = db(k)
               xb(k) = xc(k)
            if((xu(k)-xl(k)) .lt. eps) goto 100
  600       continue

            ierr = 2
          end if
  100       continue             
        end if

  800    do k = 1, ncolor
          if (du(k) .eq. 0.0) du(k) = eps
         end do
c
c        -------------------------------
c        ... use givens recursion to get
c            eigenvector and bound
c        -------------------------------
c
         do k = 1, ncolor
c         if(mynode() .eq.  0) then
c         if (j.eq.3)write(*, *) j, ' alf', alf(1,k), alf(2,k),xc(k)
c         if (j.eq.3)write(*, *) j, ' bet', bet(1,k), bet(2,k),xc(k)
c         end if

	 call givens (alf(1,k), bet(1,k), j, xc(k), z(1,k))
         res(k)  = abs(z(j,k)) * bet(j+1,k)
	 bd(k)   = res(k)
         sumi(k) = sumi(k) + float(ic(k))
c   
c        ----------------------------------
c        ... write smallest eigenvalue and
c            check for convergence
c        ---------------------------------
c
	 jmax = j
	 lanstp = j

         end do



         do k = 1, ncolor
	   if (res(k) .gt. eps) go to 5000
         end do

        go to 6000 
 5000 continue
      ierr = 1
c
 6000 do k = 1, ncolor
       eig(k)  = xc(k)
       sumi(k) = sumi(k) / float(jmax)
      end do

c
c     ------------------------------------------
c     ...  compute eigenvector after convergence
c     ------------------------------------------
c

      call copy (neqns, evec, u)

      do 6100 i = 1, ncolor
        do j = 1, size(i)
	 evec(indarray(indptr(i)+j)) = z(1,i) * u(indarray(indptr(i)+j))
	 r(indarray(indptr(i)+j))    = 0.0e0
        end do
        eval(i) = alf(1,i)
        bet(1,i) = 0.0e0
        betc(1,i) = 0.0e0
 6100   continue
c
c     ------------------------------
c     ... run Lanczos loop again for
c         computation of vector
c     -----------------------------
c
      do 8000 j = 1, jmax-1
c
       do k = 1, ncolor
	 call scal (size(k), -bet(j,k), r, indarray(indptr(k)+1))
       end do
	 call opk (xadj, adj, diag, u, r, neqns, mask)

       do k = 1, ncolor
	  alf(j,k) = dot (size(k), u, r, indarray(indptr(k)+1))
          call axpy (size(k), -alf(j,k), u, r, indarray(indptr(k)+1))
       end do
c
c        --------------------------------------
c        ... orthogonalize against 1, 1, ...
c        --------------------------------------
c
	 call ortho (neqns, r, indarray, indptr, size, ncolor) 
c
       do k = 1, ncolor
	 bet(j+1,k) = nrm2 (size(k), r, indarray(indptr(k)+1))
	 call scal (size(k), 1.0e0/bet(j+1,k), r, indarray(indptr(k)+1))
         call swap (size(k), u, r, indarray(indptr(k)+1))
	 call axpy (size(k), z(j+1,k), u, evec, indarray(indptr(k)+1))

       end do
c
 8000 continue


      return
 9600 ierr = - ierr
      if (msglvl .gt. 0) write (output, 9002) ierr
      return
c
 8050 format ( /,"average of pivotcalls", f15.2 )
 9001 format (1x, i4, 4x, e25.14, 4x, e25.14 )
 9002 format ( /, 1x, 'eigenvector computation error, ierr = ', i5 )
 9004 format (1x, 'small beta j+1 encountered  ', e20.14 )
      end
c
c     ----------------------
c
c     ---------------------
c
      subroutine opk(xadj, adj, diag, a,b,n, mask)
c
c********************************************************************
c**** multiplication by the Laplacian matrix in sparspak adjacency **
c**** format
c**** multiplication in the format b gets b + Aa ********************
c**** switch to original negative definite Laplacian because
c**** of new eigenvalue solver as of 1/19/91     ********************
c********************************************************************
c
      implicit  none
      integer sched, noffproc, n
      common /parti/sched, noffproc

      real a(n),b(n)
      integer xadj(*), adj(*), mask(*)
      real diag(n)
c
      integer i, j, k
      do 100 i = 1, n
         b(i) =  b(i) + diag(i) * a(i) 
  100 continue

      do 150 i = 1, noffproc
        b(n+i) = 0.0
  150 continue

      do 300 i = 1, n
CDIR$ IVDEP
         do 200 j = xadj(i), xadj(i+1) - 1
            if (mask(j) .gt. 0) then
              k =  adj(j)
              b(k) = b(k) + a(i)
            end if
  200    continue
  300 continue

      call fscatter_add(b(n+1), b, sched)

      return
      end
c
c     ----------------------
c
c     ---------------------
c
      subroutine ortho (neqns, r, indarray, indptr, size, ncolor) 
c
#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

      integer ncolor
      integer size(ncolor), indarray(neqns), indptr(ncolor)
      real    r(neqns)
c
c     -------------------
c     ... local variables
c     -------------------
c
      integer   i, j
      real      sum, rwork
      integer   teqns, iwork
c
c-----------------------------------------------------------------
c
c
c     -------------------------------------------
c     ... for simple component just orthogonalize
c     -------------------------------------------
c
      do j = 1, ncolor
           teqns = size(j)
           call PARTI_gisum(teqns, 1, iwork)
           srn = 1.0e0 / float(teqns)
	   sum = 0.0e0
	   do 100 i = 1, size(j)
	    sum = sum + r(indarray(indptr(j)+i))
 100   continue

       call PARTI_gssum(sum, 1, rwork)

	   sum = sum * srn
	   do 200 i = 1, size(j)
	    r(indarray(indptr(j)+i)) = r(indarray(indptr(j)+i)) - sum
 200     continue
       end do
      return
      end
c
c     ----------------------
c
c     ---------------------
c
      subroutine pivot(a,b,n,x,d)
      dimension a(*),b(*)
      d=1.
      if (n.le.0) return
      d=a(1)-x
      if(d.eq.0.) d=b(1)*1.e-6 
      if(n.eq.1) return
      do 10 i=2,n
         d=a(i)-x-b(i-1)/d
         if(d.eq.0.) d=b(i-1)*1.e-6
10       continue
      return
      end
c
c     -----------------
c
c     ----------------
c
      subroutine scal (n, a, x, indir)
      real a, x(1)
      integer indir(1)
      integer i, n
      do 100 i = 1, n
	 x(indir(i)) = a * x(indir(i))
  100 continue
      return
      end
c
