subroutine addcom ( isub, jsub, value, invprm, diag, xlnz,
     &  ixlnz, nzsub, xnzsub, neqn )

c*********************************************************************72
c
cc ADDCOM adds values to a matrix stored in compressed storage scheme.
c
c  Discussion:
c
c    ADDCOM is called once the equations and variables have
c    been reordered, to add numbers to the matrix A, which is
c    assumed to be stored according to the compressed storage
c    scheme used by the ND and QMD methods.  
c
c    The arrays DIAG and XLNZ, which are used as storage for array entries,
c    should be zeroed out before this routine is first called for a
c    particular problem.
c
c  Modified:
c
c    01 January 2009
c
c  Author:
c
c    Alan George, Joseph Liu
c
c  Reference:
c
c    Alan George, Joseph Liu,
c    Computer Solution of Large Sparse Positive Definite Systems,
c    Prentice Hall, 1981,
c    ISBN: 0131652745,
c    LC: QA188.G46.
c
c  Parameters:
c
c    Input, integer ISUB, JSUB, the row and column of the matrix to
c    which the value is to be added.  Note that the matrix is assumed
c    to be symmetric, and only the lower triangle is stored.  Hence,
c    if ISUB < JSUB, the routine ignores the request.
c
c    Input, double precision VALUE, the number to be added to A(ISUB,JSUB).
c
c    Input, integer INVPRM(NEQN), the inverse ordering, which should have
c    been created by calling PERM_INVERSE once the reordering is set.
c
c    Input/output, double precision DIAG(NEQN), the diagonal elements of the
c    matrix.
c
c    Input/output, double precision XLNZ(*), the nonzero subdiagonal elements
c    of the matrix.
c
c    Input, integer IXLNZ(NEQN+1), NZSUB(*), XNZSUB(NEQN), data structures which
c    define the compressed storage scheme.
c
c    Input, integer NEQN, the order of the matrix.
c
      integer neqn

      double precision diag(neqn)
      integer i
      integer invprm(neqn)
      integer isub
      integer ixlnz(neqn+1)
      integer j
      integer jsub
      integer k
      integer kstop
      integer kstrt
      integer ksub
      integer nzsub(*)
      double precision value
      double precision xlnz(*)
      integer xnzsub(neqn)
c
c  Figure out the current locations of the given row and column.
c
      i = invprm(isub)
      j = invprm(jsub)
c
c  if the entry is on the diagonal, update DIAG and return.
c
      if ( i .eq. j ) then
        diag(i) = diag(i) + value
        return
      endif
c
c  if the entry is above the diagonal, don't store it at all.
c
      if ( i .lt. j ) then
        return
       end if

      kstrt = ixlnz(j)
      kstop = ixlnz(j+1) - 1
 
      if ( kstop .lt. kstrt ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'ADDCOM - Fatal error!'
        write ( *, '(a)' ) '  The IXLNZ array is incorrect.'
        write ( *, '(a)' ) ' '
        write ( *, '(a,i8)' ) '  ISUB = ', isub
        write ( *, '(a,i8)' ) '  JSUB = ', jsub
        write ( *, '(a)' ) ' '
        write ( *, '(a,i8)' ) '  I = ', i
        write ( *, '(a,i8)' ) '  J = ', j
        write ( *, '(a)' ) ' '
        write ( *, '(a,i8)' ) '  IXLNZ(J) =     ', ixlnz(j)
        write ( *, '(a,i8)' ) '  IXLNZ(J+1)-1 = ', ixlnz(j+1) - 1
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  However, we must have that'
        write ( *, '(a)' ) '  IXLNZ(J) <= IXLNZ(J+1)-1.'
        stop
      end if
 
      ksub = xnzsub(j)
 
      do k = kstrt, kstop
 
        if ( nzsub(ksub) .eq. i ) then
          xlnz(k) = xlnz(k) + value
          return
        end if
 
        ksub = ksub + 1
 
      end do
 
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'ADDCOM - Fatal error!'
      write ( *, '(a)' ) '  No storage for entry ISUB, JSUB'
      write ( *, '(a,i8)' ) '  ISUB = ', isub
      write ( *, '(a,i8)' ) '  JSUB = ', jsub

      stop
      end
      subroutine addrcm ( isub, jsub, value, invprm, diag, xenv, env,
     &  ierror, neqn )

c*********************************************************************72
c
cc ADDRCM adds values to a matrix stored in the RCM scheme.
c
c  Discussion:
c
c    ADDRCM can be called, once the equations and variables have
c    been reordered for the RCM method, in order to enter values
c    into the matrix A.  Note that values are added, so that
c    the matrix storage should be zeroed out before calling
c    ADDRCM.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    31 December 2008
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Alan George, Joseph Liu,
c    Computer Solution of Large Sparse Positive Definite Systems,
c    Prentice Hall, 1981,
c    ISBN: 0131652745,
c    LC: QA188.G46.
c
c  Parameters:
c
c    Input, integer ISUB, the original row to which the entry is to be added.
c
c    Input, integer JSUB, the original column to which the entry is to be added.
c
c    Input, double precision VALUE, the value to be added to entry A(ISUB,JSUB).
c
c    Input, integer INVPRM(NEQN), the reverse ordering produced by
c    subroutine INVERSE.
c
c    Input/output, double precision DIAG(NEQN), the diagonal elements of 
c    the matrix.
c
c    Input/output, integer XENV(NEQN+1), double precision ENV(*), the envelope 
c    structure of the matrix.
c
c    Output, integer IERROR, error flag.
c    0 for no errors detected.
c    5 for error (isub,jsub) lies outside envelope.
c
c    Input, integer NEQN, the number of equations.
c
      integer neqn

      double precision diag(neqn)
      double precision env(*)
      integer i
      integer ierror
      integer invprm(neqn)
      integer isub
      integer j
      integer jsub
      integer k
      double precision value
      integer xenv(neqn+1)

      ierror = 0
      i = invprm(isub)
      j = invprm(jsub)

      if ( i .lt. j ) then
        return
       end if

      if ( i .eq. j ) then
        diag(i) = diag(i) + value
        return
      endif
 
      k = xenv(i+1) - i + j
 
      if ( k .lt. xenv(i) ) then
        ierror = 5
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'ADDRCM - Fatal error!'
        write ( *, '(a)' ) '  Indices outside of envelope.'
        write ( *, '(a)' ) '  ISUB = ', isub, ' JSUB = ', jsub
        stop
      end if

      env(k) = env(k) + value
 
      return
      end
      subroutine addrhs ( invprm, isub, neqn, rhs, value )

c*********************************************************************72
c
cc ADDRHS adds a quantity to a specific entry of the right hand side.
c
c  Discussion:
c
c    ADDRHS is called after the variables and equations have been
c    reordered.  addrhs allows the user to add the number value
c    to the i-th entry of the right hand side, which has been
c    permuted to a new location.  please zero out rhs once before
c    calling this routine to build up a right hand side.
c
c  Modified:
c
c    01 January 2009
c
c  Parameters:
c
c    Input, integer INVPRM(NEQN), the inverse ordering vector which must 
c    have been set up before calling ADDRHS.  INVPRM is constructed
c    by subroutine INVERSE.
c
c  isub   - the original index, before reordering, of the location
c           in rhs to which the number value is to be added.
c
c  neqn   - the number of entries in rhs and invprm.
c
c  rhs    - the right hand side vector.
c
c  value  - the number to be added to rhs.
c
      integer neqn

      integer i
      integer invprm(neqn)
      integer isub
      double precision rhs(neqn)
      double precision value

      i = invprm(isub)
 
      if ( i .lt. 0 ) then
        return
      end if

      if ( neqn .lt. i ) then
        return
      end if

      rhs(i) = rhs(i) + value
 
      return
      end
      subroutine addrqt ( isub, jsub, value, invprm, 
     &  diag, xenv, env, xnonz, nonz, nzsubs, neqn )

c*********************************************************************72
c
cc ADDRQT adds values to a matrix stored in the implicit block storage scheme.
c
c  Discussion:
c
c    addrqt adds a number to the (i,j)-th position of a matrix
c    stored using the implicit block storage scheme.  since addrqt
c    only adds new values to those currently in storage, the
c    space used to store the matrix must be initialized to 0
c    before numerical values are supplied.  addrqt can be used with the
c    rqt and 1wd methods.
c
c  Modified:
c
c    01 January 2009
c
c  Reference:
c
c    Alan George, Joseph Liu,
c    Computer Solution of Large Sparse Positive Definite Systems,
c    Prentice Hall, 1981,
c    ISBN: 0131652745,
c    LC: QA188.G46.
c
c  input
c
c    Input, integer ISUB, the row index of the number to be added.
c
c    Input, integer JSUB, the column index of the number to be added.
c
c    diag   - an array containing the diagonal elements of the matrix.
c
c    value  - the number to be added.
c
c    invprm - invprm(i) is the new position of the variable whose
c    original number is i.  invp can be set up by
c    calling subroutine inverse after the reordering
c    routines are called.
c
c    xenv, env - a pair of arrays which contain the envelope
c    structure of the diagonal blocks.
c
c    xnonz, nonz, nzsubs - levels structure containing the off-
c    block diagonal parts of the rows of the lower
c    triangle of the original matrix.
c
c  output
c
c    also, diag, env or nonz has been modified by having the number
c    stored in value added in the (i,j)-th position.
c
      integer neqn

      double precision diag(neqn)
      double precision env(*)
      integer i
      integer invprm(neqn)
      integer isub
      integer j
      integer jsub
      integer k
      integer kstop
      integer kstrt
      double precision nonz(*)
      integer nzsubs(*)
      double precision value
      integer xenv(neqn+1)
      integer xnonz(neqn)

      i = invprm(isub)
      j = invprm(jsub)
 
      if ( i .eq. j ) then
        diag(i) = diag(i) + value
        return
      end if
 
      if ( i .lt. j ) then
        return
      end if
c
c  The value goes within the diagonal envelope.
c
      k = xenv(i+1) - i + j
 
      if ( k .ge. xenv(i) ) then
        env(k) = env(k) + value
        return
      endif
c
c  The value goes outside the diagonal blocks.
c
      kstrt = xnonz(i)
      kstop = xnonz(i+1)-1
 
      do k = kstrt, kstop
 
        if ( nzsubs(k) .eq. j ) then
          nonz(k) = nonz(k) + value
          return
        end if
 
      end do
 
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'ADDRQT - Fatal error!'
      write ( *, '(a)' ) '  lack of storage!'
      write ( *, * ) '  i=',i
      write ( *, * ) '  isub=',isub
      write ( *, * ) '  j=',j
      write ( *, * ) '  jsub=',jsub

      stop
      end
      subroutine adj_env_size ( n, adj_row, adj_num, adj, perm, 
     &  perm_inv, env_size )

c*********************************************************************72
c
cc ADJ_ENV_SIZE computes the envelope size for an adjacency structure.
c
c  Discussion:
c
c    The matrix is assumed to be symmetric.
c
c    The variables are assumed to have been permuted.
c
c  Modified:
c
c    24 February 2007
c
c  Author:
c
c    Alan George, Joseph Liu
c
c  Reference:
c
c    Alan George, Joseph Liu,
c    Computer Solution of Large Sparse Positive Definite Systems,
c    Prentice Hall, 1981,
c    ISBN: 0131652745,
c    LC: QA188.G46.
c
c  Parameters:
c
c    Input, integer N, the order of the matrix.
c
c    Input, integer ADJ_ROW(N+1).  Information about row I is stored
c    in entries ADJ_ROW(I) through ADJ_ROW(I+1)-1 of ADJ.
c
c    Input, integer ADJ_NUM, the number of adjacency entries.
c
c    Input, integer ADJ(ADJ_NUM), the adjacency structure. 
c    For each row, the column indices of the nonzero entries.
c
c    Input, integer PERM(N), integer PERM_INV(N), the permutation
c    and inverse permutation.
c
c    Output, integer ENV_SIZE, the number of cells in the envelope.
c    You could think of this number as the sum of the
c    bandwidths of each row.
c
      implicit none

      integer adj_num
      integer n

      integer add
      integer adj(adj_num)
      integer adj_row(n+1)
      integer col
      integer env_size
      integer i
      integer j
      integer perm(n)
      integer perm_inv(n)
      integer row

      env_size = 0
 
      do i = 1, n

        row = perm(i)

        add = 0

        do j = adj_row(row), adj_row(row+1) - 1
          col = perm_inv(adj(j))
          if ( col .lt. i ) then
            add = max ( add, i - col )
          end if
        end do

        env_size = env_size + add

      end do

      return
      end
      subroutine adj_print ( n, nadj, xadj, adj )

c*********************************************************************72
c
cc ADJ_PRINT prints the adjacency information stored in XADJ and ADJNCY.
c
c  Discussion:
c
c    The list has the form:
c
c    row   nonzeros
c
c    1       2   5   9
c    2       7   8   9   15   78   79   81  86  91  99
c          100 103
c    3      48  49  53
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    18 December 2002
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    input, integer N, the number of equations.
c
c    input, integer NADJ, the dimension of adjncy.
c
c    input, integer XADJ(N+1), organizes the entries of adjncy
c    into rows.  The entries for row i are in entries xadj(i)
c    through XADJ(I+1)-1.
c
c    input, integer ADJ(NADJ), the adjacency structure, which contains,
c    for each row, the column indices of the nonzero entries.
c
       implicit none

       integer nadj
       integer n

       integer adj(nadj)
       integer i
       integer j
       integer jhi
       integer jlo
       integer jmax
       integer jmin
       integer xadj(n+1)

       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'ADJ_PRINT'
       write ( *, '(a)' ) '  show adjacency structure of sparse matrix.'
       write ( *, '(a,i6)' ) '  the matrix order is ', n
       write ( *, '(a,i6)' ) '  the number of entries is ', nadj
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) '  row       nonzeros '
       write ( *, '(a)' ) ' '

       do i = 1, n

         jmin = xadj(i)
         jmax = xadj(i+1) - 1
         do jlo = jmin, jmax, 10
           jhi = min ( jlo + 9, jmax )
           if ( jlo .eq. jmin ) then
             write ( *, '(i6,6x,10i6)' ) i, (adj(j),j=jlo,jhi)
           else
             write ( *, '(6x,6x,10i6)' ) ( adj(j),j=jlo,jhi)
           end if
         end do

       end do

       return
      end
      subroutine adj_show ( adj, iband, invprm, nadj, n, perm, xadj )

c*********************************************************************72
c
cc ADJ_SHOW displays a symbolic picture of a matrix.
c
c  Discussion:
c
c    The matrix is defined by the adjacency information in xadj and adj,
c    with a possible permutation through perm and invprm.  the routine
c    also computes the bandwidth and the size of the envelope.
c
c    if no permutation has been done, you must set invprm(i) = perm(i) = i
c    before calling this routine.  otherwise, you must call perm_inverse to
c    get the inverse permutation invprm before calling showmat.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    01 January 2009
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Alan George, Joseph Liu,
c    Computer Solution of Large Sparse Positive Definite Systems,
c    Prentice Hall, 1981,
c    ISBN: 0131652745,
c    LC: QA188.G46.
c
c  Parameters:
c
c    Input, integer ADJ(NADJ), the adjacency structure.
c    for each row, it contains the column indices of the nonzero entries.
c
c    Output, integer IBAND, the bandwidth of the matrix.
c
c    Input, integer PERM(N), INVPRM(N), the permutation and inverse permutation.
c
c    Input, integer NADJ, the number of adjacency entries in ADJ.
c
c    Input, integer N, the order of the matrix.
c
c    Input, integer XADJ(N+1).  information about row I is stored
c    in entries XADJ(I) through XADJ(I+1)-1 of ADJ.
c
      implicit none

      integer n_max
      parameter ( n_max = 100 )

      integer nadj
      integer n

      integer add
      integer adj(nadj)
      character band(n_max)
      integer col
      integer i
      integer iband
      integer invprm(n)
      integer j
      integer jhi
      integer jlo
      integer k
      integer nonz
      integer perm(n)
      integer xadj(n+1)

      iband = 0
      nonz = 0

      if ( n_max .lt. n ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'ADJ_SHOW - Fatal error!'
        write ( *, '(a)' ) '  N is too large!'
        write ( *, '(a,i6)' ) '  Maximum legal value is ', n_max
        write ( *, '(a,i6)' ) '  Your input value was ', n
        stop
      end if

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'ADJ_SHOW:'
      write ( *, '(a)' ) '  display nonzero structure of matrix.'
      write ( *, '(a)' ) ' '

      do i = 1, n

        do k = 1, n
          band(k) = ' '
        end do

        band(i) = 'x'

        do j = xadj(perm(i)), xadj(perm(i)+1)-1
          col = invprm(adj(j))
          if ( col .lt. i ) then
            nonz = nonz + 1
          end if
          iband = max ( iband, i-col )
          band(col) = 'x'
        end do

        write ( *, '(i6,1x,100a1)' ) i, ( band(j), j = 1, n )

      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a,i6)' ) 'lower bandwidth = ', iband
      write ( *, '(a,i6,a)' ) 
     &  'lower envelope contains ', nonz, ' nonzeros.'

      return
      end
      subroutine bshufl ( xadj, adjncy, perm, nblks, xblk,
     &  bnum, mask, subg, xls )

c*********************************************************************72
c
cc BSHUFL renumbers the nodes of each block to reduce its envelope.  
c
c  Discussion:
c
c    This routine renumbers the nodes of each block
c    so as to reduce its envelope.
c    nodes in a block with no neighbors in previous
c    blocks are renumbered by subrcm before the others.
c
c  Modified:
c
c    01 January 2009
c
c  Reference:
c
c    Alan George, Joseph Liu,
c    Computer Solution of Large Sparse Positive Definite Systems,
c    Prentice Hall, 1981,
c    ISBN: 0131652745,
c    LC: QA188.G46.
c
c  input parameters 
c
c  (xadj, adjncy) - the graph adjacency structure.
c
c  (nblks, xblk ) - the tree partitioning.
c
c  updated parameter 
c  perm - the permutation vector. on return, it contains
c  the new permutation.
c
c  working vectors -
c  bnum - stores the block number of each variable.
c  mask - mask vector used to prescribe a subgraph.
c  subg - vector used to contain a subgraph.
c  xls - index vector to a level structure.
c
      integer adjncy(*)
      integer bnum(*), mask(*), perm(*)
      integer subg(*), xblk(*), xls(*)
      integer xadj(*), i, ip, istop, istrt, j
      integer jstrt, jstop, k, nabor, nblks, nbrblk
      integer node, nsubg

      if ( nblks .le. 0 ) then
        return
      end if
c
c  initialization.  find the block number for each
c  variable and initialize the vector mask.
c
      do k = 1, nblks
        istrt = xblk(k)
        istop = xblk(k+1) - 1
        do i = istrt, istop
          node = perm(i)
          bnum(node) = k
          mask(node) = 0
        end do
      end do
c
c  for each block, find those nodes with no neighbors
c  in previous blocks and accumulate them in subg.
c  they will be renumbered before others in the block.
c
      do k = 1, nblks

        istrt = xblk(k)
        istop = xblk(k+1) - 1
        nsubg = 0

        do i = istrt, istop

          node = perm(i)
          jstrt = xadj(node)
          jstop = xadj(node+1) - 1

          if ( jstrt .le. jstop ) then

            do j = jstrt, jstop
              nabor = adjncy(j)
              nbrblk = bnum(nabor)
              if (nbrblk .lt. k) go to 400
            end do
            nsubg = nsubg + 1
            subg(nsubg) = node
            ip = istrt + nsubg - 1
            perm(i) = perm(ip)

          end if

400       continue

        end do
c
c  call subrcm to renumber the subgraph stored
c  in (nsubg, subg).
c
        if ( 0 .lt. nsubg ) then
          call subrcm ( xadj, adjncy, mask, nsubg,
     &      subg, perm(istrt),xls)
        endif
        
      end do

      return
      end
      subroutine copysi ( n, a, b )

c*********************************************************************72
c
cc COPYSI copies the N integer elements from the vector A to B.
c
c  Modified:
c
c    01 January 2009
c
c  Parameters:
c
c    Input, integer N, the number of entries to copy.
c
c    Input, integer A(N), the vector to copy.
c
c    Output, integer B(N), the copied vector.
c
      integer n

      integer a(n)
      integer b(n)
      integer i

      do i = 1, n
        b(i) = a(i)
      end do

      return
      end
      subroutine degree ( root, xadj, adjncy, mask, deg, ccsize, ls )

c*********************************************************************72
c
cc DEGREE computes node degrees in a connected component, for the RCM method.
c
c  Discussion:
c
c    This routine computes the degrees of the nodes in the connected component 
c    specified by MASK and ROOT.
c
c    Nodes for which MASK is zero are ignored.
c
c  Modified:
c
c    01 January 2009
c
c  Reference:
c
c    Alan George, Joseph Liu,
c    Computer Solution of Large Sparse Positive Definite Systems,
c    Prentice Hall, 1981,
c    ISBN: 0131652745,
c    LC: QA188.G46.
c
c  Parameters:
c
c    Input, integer ROOT, the node that defines the component.
c
c    (xadj, adjncy) - adjacency structure pair.
c
c    Input, integer MASK(N), specifies a section subgraph.
c
c    Output, integer DEGREE(N), the degrees of the nodes in the component.
c
c    Output, integer CCSIZE, the size of the component.
c
c    Workspace, integer LS(N), used to store the nodes of the
c    component level by level.
c
      integer adjncy(*)
      integer ccsize
      integer deg(*)
      integer i
      integer ideg
      integer j
      integer jstop
      integer jstrt
      integer lbegin
      integer ls(*)
      integer lvlend
      integer lvsize
      integer mask(*)
      integer nbr
      integer node
      integer root
      integer xadj(*)
c
c  Initialization.
c
c  the array xadj is used as a temporary marker to
c  indicate which nodes have been considered so far
c
      ls(1) = root
      xadj(root) = -xadj(root)
      lvlend = 0
      ccsize = 1
c
c  lbegin is the pointer to the beginning of the curi
c  level, and lvlend points to the end of this level.
c
100   continue

      lbegin = lvlend + 1
      lvlend = ccsize
c
c  find the degrees of nodes in the current level,
c  and at the same time, generate the next level.
c
      do i = lbegin, lvlend
        node = ls(i)
        jstrt = -xadj(node)
        jstop = iabs(xadj(node+ 1) ) - 1
        ideg = 0
        if ( jstop .lt. jstrt ) go to 300
        do 200 j = jstrt, jstop
          nbr = adjncy(j)
          if ( mask(nbr) .eq. 0 ) go to 200
          ideg = ideg + 1
          if ( xadj(nbr) .lt. 0 ) go to 200
          xadj(nbr) = -xadj(nbr)
          ccsize = ccsize + 1
          ls(ccsize) = nbr
200     continue
300     deg(node) = ideg
      end do
c
c  Compute the current level width.
c  if it is nonzero , generate another level.
c
      lvsize = ccsize - lvlend
      if ( lvsize .gt. 0 ) go to 100
c
c  Reset xadj to its correct sign and return.
c
      do i = 1, ccsize
        node = ls(i)
        xadj(node) = - xadj(node)
      end do

      return
      end
      subroutine elslv ( neqns, xenv, env, diag, rhs )

c*********************************************************************72
c
cc ELSLV solves a lower triangular system stored in the envelope format.
c
c  Discussion:
c
c    This routine solves a lower triangular system l x - rhs.
c
c    The factor l is stored in the envelope format.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c    Input, integer NEQNS, the number of equations.
c
c  (xenv, env) - array pair for the envelope of l.
c  diag - array for the diagonal of l.
c
c  updated parameter 
c  rhs - on input, it contains the right hand vector.
c  on return, it contains the solution vector.
c
      double precision diag(*)
      double precision env(*)
      double precision rhs(*)
      double precision s
      integer xenv(*), i, iband, ifirst, k, kstop
      integer kstrt, l, last, neqns
c
c  find the position of the first nonzero in rhs and put it in ifirst.
c
      ifirst = 0
100   ifirst = ifirst + 1
      if ( rhs(ifirst) .ne. 0.0D+00 ) go to 200
      if ( ifirst .lt. neqns ) go to 100
      return
200   last = 0
c
c  last contains the position of the most recently
c  computed nonzero component of the solution.
c
      do i = ifirst, neqns

        iband = xenv(i+1) - xenv(i)
        if ( iband .ge. i ) iband = i - 1
        s = rhs(i)
        l = i - iband
        rhs(i) = 0.0D+00
c
c  row of the envelope is empty, or corresponding
c  components of the solution are all zeros.
c
        if ( iband .eq. 0 .or. last .lt. l ) go to 400
        kstrt = xenv(i+1) - iband
        kstop = xenv(i+1) - 1
        do k = kstrt, kstop
          s = s - env(k)*rhs(l)
          l = l + 1
        end do

400     continue

        if ( s .ne. 0.0D+00 ) then
          rhs(i) = s / diag(i)
          last = i
        end if

      end do

      return
      end
      subroutine esfct ( neqns, xenv, env, diag, flag )

c*********************************************************************72
c
cc ESFCT factors a positive definite envelope matrix into L * L'.
c
c  Discussion:
c
c    The algorithm used is the standard bordering method.
c
c  Modified:
c
c    31 December 2008
c
c  Author:
c
c    Alan George, Joseph Liu
c
c  Reference:
c
c    Alan George, Joseph Liu,
c    Computer Solution of Large Sparse Positive Definite Systems,
c    Prentice Hall, 1981,
c    ISBN: 0131652745,
c    LC: QA188.G46.
c
c  Parameters:
c
c    Input, integer NEQNS, the number of equations.
c
c    Input, integer XENV(N+1), the envelope index vector.
c
c    Input/output, double precision ENV(*), on input, the envelope of A, 
c    and on output the envelope of the factor L.
c
c    Input/output, double precision DIAG(NEQNS), on input the diagonal of A, 
c    and on output the diagonal of L.
c
c    Output, integer FLAG, error flag.
c    0, no error, the factorization was carried out.
c    1, the matrix is not positive definite.
c
      implicit none

      integer neqns

      double precision diag(neqns)
      double precision env(*)
      integer flag
      integer i
      integer iband
      integer ifirst
      integer ixenv
      integer j
      integer jstop
      double precision temp
      integer xenv(neqns+1)

      flag = 0

      if ( diag(1) .le. 0.0D+00 ) then
        write ( *, '(a)' ) ' '
        write ( *, * ) 'ESFCT: Fail for DIAG(1)'
        flag = 1
        return
      end if

      diag(1) = sqrt ( diag(1) )
c
c  Loop over rows 2:NEQNS of the matrix
c
      do i = 2, neqns

        ixenv = xenv(i)
        iband = xenv(i+1) - ixenv
c
c  Compute row I of the triangular factor.
c
        temp = diag(i)

        if ( iband .ne. 0 ) then

          ifirst = i - iband

          call elslv ( iband, xenv(ifirst), env, 
     &      diag(ifirst), env(ixenv) )

          jstop = xenv(i+1) - 1

          do j = ixenv, jstop
            temp = temp - env(j) * env(j)
          end do

        end if

        if ( temp .le. 0.0D+00 ) then
          write ( *, '(a)' ) ' '
          write ( *, * ) 'ESFCT: Fail for DIAG(', i ,')'
          flag = 1
          return
        end if

        diag(i) = sqrt ( temp )

      end do

      return
      end
      subroutine euslv ( neqns, xenv, env, diag, rhs )

c*********************************************************************72
c
cc EUSLV solves an upper triangular system u x - rhs. 
c
c  Discussion:
c
c    the factor u is stored in the envelope format.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters -
c
c    Input, integer NEQNS, the number of equations.
c
c  (xenv, env) - array pair for the envelope of u.
c  diag - array for the diagonal of u.
c
c  updated parameter -
c
c  rhs - on input, it contains the right hand side.
c  on output, it contains the solution vector.
c
      double precision diag(*)
      double precision env(*), rhs(*), s
      integer xenv(*), i, iband, k, kstop, kstrt, l
      integer neqns

      i = neqns + 1

100   continue

      i = i - 1 

      if ( i .eq. 0 ) then
        return
      end if

      if ( rhs(i) .eq. 0.0D+00 ) then
        go to 100
      end if

      s = rhs(i) / diag(i)
      rhs(i) = s

      iband = xenv(i+1) - xenv(i)

      if ( iband .ge. i ) then
        iband = i - 1
      end if

      if ( iband .eq. 0 ) then
        go to 100
      end if

      kstrt = i - iband
      kstop = i - 1
      l = xenv(i+1) - iband

      do k = kstrt, kstop
        rhs(k) = rhs(k) - s * env(l)
        l = l + 1
      end do

      go to 100
      end
      subroutine fn1wd ( root, xadj, adjncy, mask, nsep, sep, nlvl, 
     &  xls, ls )

c*********************************************************************72
c
cc FN1WD finds one-way dissectors of
c  a connected component specified by mask and root.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters
c
c  root - a node that defines (along with mask) the
c  component to be processed.
c
c  (xadj, adjncy) - the adjacency structure.
c
c  output parameters 
c  nsep - number of nodes in the one-way dissectors.
c  sep - vector containing the dissector nodes.
c
c  updated parameter 
c  mask - nodes in the dissector have their mask values
c  set to zero.
c
c  working parameters
c  (xls, ls) - level structure used by the routine fnroot.
c
      integer adjncy(*)
      integer ls(*)
      integer mask(*), sep(*), xls(*)
      integer xadj(*), i, j, k, kstop, kstrt, lp1beg, lp1end
      integer lvl, lvlbeg, lvlend, nbr, nlvl, node
      integer nsep, root
      double precision deltp1, fnlvl, width

      call fnroot ( root, xadj, adjncy, mask,
     &  nlvl, xls, ls )

      fnlvl = dble ( nlvl )
      nsep = xls(nlvl + 1) - 1
      width = dble ( nsep ) / fnlvl
      deltp1 = 1.0D+00 
     &  + sqrt ( ( 3.0D+00 * width + 13.0D+00 ) / 2.0D+00 )

      if ( nsep .ge. 50 .and. deltp1 .le. 0.5D+00 * fnlvl ) then
        go to 300
      end if
c
c  the component is too small, or the level structure
c  is very long and narrow. return the whole component.
c
      do i = 1, nsep
        node = ls(i)
        sep(i) = node
        mask(node) = 0
      end do

      return
c
c  find the parallel dissectors.
c
300   nsep = 0
      i = 0
400   i = i + 1
      lvl = int ( dble ( i ) * deltp1 + 0.5D+00 )

      if ( lvl .ge. nlvl ) then
        return
      end if

      lvlbeg = xls(lvl)
      lp1beg = xls(lvl + 1)
      lvlend = lp1beg - 1
      lp1end = xls(lvl + 2) - 1

      do j = lp1beg, lp1end
        node = ls(j)
        xadj(node) = - xadj(node)
      end do
c
c  nodes in level lvl are chosen to form dissector.
c  include only those with neighbors in lvl+1 level.
c  xadj is used temporarily to mark nodes in lvl+1
c
      do j = lvlbeg, lvlend
        node = ls(j)
        kstrt = xadj(node)
        kstop = iabs(xadj(node+1)) - 1
        do 600 k = kstrt, kstop
          nbr = adjncy(k)
          if ( xadj(nbr) .gt. 0 ) go to 600
          nsep = nsep + 1
          sep(nsep) = node
          mask(node) = 0
          go to 700
600     continue

700     continue

      end do

      do j = lp1beg, lp1end
        node = ls(j)
        xadj(node) = - xadj(node)
      end do

      go to 400
      end
      subroutine fnbenv ( xadj, adjncy, perm, invp, nblks, xblk,
     &  xenv, envsze, smask, marker, rchset )

c*********************************************************************72
c
cc FNBENV finds the exact envelope structure of the diagonal blocks of 
c  the cholesky factor of a permuted partitioned matrix.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c  (xadj, adjncy) - adjacency structure of the graph.
c  (perm, invp) - the permutation vector and its inverse.
c  (nblks, xblk) - the partitioning.
c
c  output parameters _
c  xenv - the envelope index vector.
c  envsze - the size of the envelope.
c
c  working parameters 
c  smask - marks nodes that have been considered.
c  marker - is used by routine reach.
c  rchset - is used by the subroutine reach.
c  stores both reachable and neighborhood sets.
c
      integer adjncy(*)
      integer invp(*)
      integer marker(*)
      integer perm(*)
      integer rchset(*)
      integer smask(*), xblk(*)
      integer xadj(*), xenv(*), blkbeg, blkend, i
      integer ifirst, inhd, k, envsze, nblks, neqns
      integer newnhd, nhdsze, node, rchsze
c
c  Initialization.
c
      neqns = xblk(nblks+1) - 1
      envsze = 1

      do i = 1, neqns
        smask(i) = 0
        marker(i) = 1
      end do
c
c  loop over the blocks 
c
      do k = 1, nblks

        nhdsze = 0
        blkbeg = xblk(k)
        blkend = xblk(k+1) - 1

        do i = blkbeg, blkend
          node = perm(i)
          marker(node) = 0
        end do
c
c  loop through the nodes in current block ...
c
        do i = blkbeg, blkend
          node = perm(i)
          call reach ( node, xadj, adjncy, smask,
     &      marker, rchsze, rchset(blkbeg),
     &      newnhd, rchset(nhdsze+1) )
          nhdsze = nhdsze + newnhd
          ifirst = marker(node)
          ifirst = invp(ifirst)
          xenv(i) = envsze
          envsze = envsze + i - ifirst
         end do
c
c  reset marker values of nodes in nbrhd set.
c
        do inhd = 1, nhdsze
          node = rchset(inhd)
          marker(node) = 0
        end do
c
c  reset marker and smask values of nodes in
c  the current block.
c
        do i = blkbeg, blkend
          node = perm(i)
          marker(node) = 0
          smask(node) = 1
        end do

      end do

      xenv(neqns+1) = envsze
      envsze = envsze - 1

      return
      end
      subroutine fndsep ( root, xadj, adjncy, mask, nsep, sep,
     &  xls, ls )

c*********************************************************************72
c
cc FNDSEP finds a small separator for a connected component.
c
c  Discussion:
c
c    The connected component is specified by MASK.
c
c  Modified:
c
c    24 July 2012
c
c  Reference:
c
c    Alan George, Joseph Liu,
c    Computer Solution of Large Sparse Positive Definite Systems,
c    Prentice Hall, 1981,
c    ISBN: 0131652745,
c    LC: QA188.G46.
c
c  Parameters:
c
c    Input, integer ROOT, the node that determines the masked component.
c
c    Input, integer XADJ(*), integer ADJNCY(*), the adjacency structure pair.
c
c    Input/output, integer MASK(*), nodes in the separator will have their mask
c    values set to zero.
c
c    Output, integer NSEP, the number of variables in the separator.
c
c    Output, integer SEP(*), the separator nodes.
c
c    Workspace, integer XLS(*), LS(*), the level structure pair for level 
c    structure found by FNROOT.
c
      implicit none

      integer adjncy(*)
      integer i
      integer j
      integer jstop
      integer jstrt
      integer ls(*)
      integer mask(*)
      integer midbeg
      integer midend
      integer midlvl
      integer mp1beg
      integer mp1end
      integer nbr
      integer nlvl
      integer node
      integer nsep
      integer root
      integer sep(*)
      integer xls(*)
      integer xadj(*)
c
c  Determine the level structure associated with ROOT.
c
      call rootls ( root, xadj, adjncy, mask, nlvl, xls, ls )
c
c  If the number of levels is less than 3, return the
c  whole component as the separator.
c
      if ( nlvl .lt. 3 ) then

        nsep = xls(nlvl+1) - 1
        do i = 1, nsep
          node = ls(i)
          sep(i) = node
          mask(node) = 0
        end do

        return

      end if
c
c  Find the middle level of the rooted level structure.
c
      midlvl = ( nlvl + 2 ) / 2
      midbeg = xls(midlvl)
      mp1beg = xls(midlvl + 1)
      midend = mp1beg - 1
      mp1end = xls(midlvl+2)
c
c  The separator is obtained by including only those
c  MIDDLE-level nodes with neighbors in the MIDDLE+1
c  level.  XADJ is used temporarily to mark those
c  nodes in the MIDDLE+1 level.
c
      do i = mp1beg, mp1end
        node = ls(i)
        xadj(node) = - xadj(node)
      end do

      nsep = 0
      
      do i = midbeg, midend
        node = ls(i)
        jstrt = xadj(node)
        jstop = iabs ( xadj(node+1) ) - 1

        do j = jstrt, jstop

          nbr = adjncy(j)

          if ( xadj(nbr) .le. 0 ) then
            nsep = nsep + 1
            sep(nsep) = node
            mask(node) = 0
            go to 10
          end if

        end do

10      continue

      end do
c
c  Reset XADJ to its correct sign.
c
      do i = mp1beg, mp1end
        node = ls(i)
        xadj(node) = - xadj(node)
      end do

      return
      end
      subroutine fnenv ( neqns, xadj, adjncy, perm, invp, xenv, envsze, 
     &  bandw )

c*********************************************************************72
c
cc FNENV finds the envelope structure of a permuted matrix.
c
c  Modified:
c
c    01 January 2009
c
c  Parameters:
c
c    Input, integer NEQNS, the number of equations.
c
c    (xadj, adjncy) - array pair containing the adjacency
c    structure of the graph of the matrix.
c
c    perm,invp - arrays containing permutation data about
c    the reordered matrix.
c
c  output parameters 
c
c    xenv - index vector for the level structure
c    to be used to store the lower (or upper)
c    envelope of the reordered matrix.
c
c    envsze - is equal to xenv(neqns+1) - 1.
c
c    bandw - bandwidth of the reordered matrix.
c
      integer adjncy(*)
      integer invp(*)
      integer perm(*)
      integer xadj(*), xenv(*), bandw, i, iband
      integer ifirst, iperm, j, jstop, jstrt, envsze
      integer nabor, neqns

      bandw = 0
      envsze = 1

      do i = 1, neqns

        xenv(i) = envsze
        iperm = perm(i)
        jstrt = xadj(iperm)
        jstop = xadj(iperm + 1) - 1
c
c  find the first nonzero in row i.
c
        ifirst = i

        do j = jstrt, jstop

          nabor = adjncy(j)
          nabor = invp(nabor)
          ifirst = min ( ifirst, nabor )

        end do

        iband = i - ifirst
        envsze = envsze + iband
        bandw = max ( bandw, iband )

      end do

      xenv(neqns+1) = envsze
      envsze = envsze - 1

      return
      end
      subroutine fnlvls ( root, xadj, adjncy, nodlvl,
     &  nlvl, xls, ls )

c*********************************************************************72
c
cc FNLVLS generates a rooted level structure for
c  a masked connected subgraph, rooted at a pseudo
c  peripheral node. the level numbers are recorded.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c  (xadj, adjncy) - the adjacency structure.
c
c  output parameters 
c  nlvl - number of levels in the level structure found
c  (xls, ls) - the level structure returned.
c
c  updated parameters 
c  root - on input, with the array nodlvl, specifies
c  the component whose pseudo-peripheral node is
c  to be found. on output, it contains that node
c  nodlvl - an inpht it cpfrlftes a sfnttan subgrapu
c  on return, it contains the node level numbers.
c
      integer adjncy(*)
      integer ls(*)
      integer nodlvl(*)
      integer xls(*)
      integer xadj(*), j, lbegin, lvl, lvlend, nlvl
      integer node, root

      call fnroot ( root, xadj, adjncy, nodlvl,
     &  nlvl, xls, ls )
     
      do lvl = 1, nlvl
        lbegin = xls(lvl)
        lvlend = xls(lvl + 1) - 1
        do j = lbegin, lvlend
          node = ls(j)
          nodlvl(node) = lvl
        end do
      end do

      return
      end
      subroutine fnofnz ( xadj, adjncy, perm, invp,
     &  nblks, xblk, xnonz, nzsubs, maxnz )

c*********************************************************************72
c
cc FNOFNZ finds the column subscripts of
c  the off-block-diagonal nonzeros in the lower triangle
c  of a partitioned matrix.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c  (xadj, adjncy) - adjacency structure pair for the graph.
c  (perm, invp) - the permutation vectors.
c  (nblks, xblk) - the block partitioning.
c
c  output parameters 
c  (xnonz, nzsubs) - the column subscripts of the nonzeros
c  of a to the left of the diagonal blocks are
c  stored row by row in continguous locations in the
c  array nzsubs. xnonz is the ~
c
c  index vector to it.
c
c  updated parameter 
c  maxnz - on input, it contains the size of the vector
c  nzsubs; and on output, the number of nonzeros
c
      integer adjncy(*)
      integer invp(*)
      integer nzsubs(*)
      integer perm(*), xblk(*)
      integer xadj(*), xnonz(*), blkbeg, blkend, i, j
      integer jperm, jxnonz, k, kstop, kstrt, maxnz
      integer nabor, nblks, nzcnt

      nzcnt = 1
c
c  loop over the blocks.
c
      do i = 1, nblks

        blkbeg = xblk(i)
        blkend = xblk(i+1) - 1
c
c  loop over the rows of the i-th block.
c
        do j = blkbeg, blkend

          xnonz(j) = nzcnt
          jperm = perm(j)
          kstrt = xadj(jperm)
          kstop = xadj(jperm+1) - 1
          if ( kstrt .gt. kstop ) go to 200
c
c  loop over the nonzeros of row j.
c
          do k = kstrt, kstop
            nabor = adjncy(k)
            nabor = invp(nabor)
c
c  check to see if it is to the left of the i-th diagonal block.
c
            if ( nabor .lt. blkbeg ) then
              if ( nzcnt .le. maxnz ) nzsubs(nzcnt) = nabor
              nzcnt = nzcnt + 1
            end if

          end do
c
c  sort the subscripts of row j
c
          jxnonz = xnonz(j)
          if ( nzcnt - 1 .le. maxnz ) then
            call sorts1 (nzcnt - jxnonz, nzsubs(jxnonz))
          endif

200       continue

        end do

      end do

      if ( 0 < nblks ) then
        xnonz(blkend+1) = nzcnt
      end if

      maxnz = nzcnt - 1

      return
      end
      subroutine fnroot ( root, xadj, adjncy, mask, nlvl, xls, ls)

c*********************************************************************72
c
cc FNROOT implements a modified version of the
c  scheme by gibbs, poole, and stockmeyer to find pseudo-
c  peripheral nodes it determines such a node for the
c  section subgraph specified by mask and root.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters -
c
c  (xadj, adjncy) - adjacency structure pair for the graph.
c  mask - specifies a section subgraph. nodes for which
c  mask is zero are ignored by fnroot.
c
c  updated parameter 
c  root - on input, it (along with mask) defines the
c
c
c  component for which a pseudo-peripheral node is
c  to be found. on output, it is the node obtained
c
c  output parameters 
c  nlvl - is the number of levels in the level structure
c  rooted at the node root.
c  (xls,ls) - the level structure array pair containing
c  the level structure found.
c
      integer adjncy(*)
      integer ls(*)
      integer mask(*), xls(*)
      integer xadj(*), ccsize, j, jstrt, k, kstop, kstrt
      integer mindeg, nabor, ndeg, nlvl, node, nunlvl, root
c
c  determine the level structure rooted at root.
c
      call rootls ( root, xadj, adjncy, mask, nlvl, xls, ls)
      ccsize = xls(nlvl+1) - 1

      if ( nlvl .eq. 1 .or. nlvl .eq. ccsize ) then
        return
      end if
c
c  pick a node with minimum degree from the last level.
c
100   continue

      jstrt = xls(nlvl)
      mindeg = ccsize
      root = ls(jstrt)
      if ( ccsize .eq. jstrt ) go to 400

      do j = jstrt, ccsize

        node = ls(j)
        ndeg = 0
        kstrt = xadj(node)
        kstop = xadj(node+1) - 1

        do k = kstrt, kstop
          nabor = adjncy(k)
          if ( mask(nabor) .gt. 0 ) then
            ndeg = ndeg + 1
          end if
        end do

        if ( ndeg .lt. mindeg ) then
          root = node
          mindeg = ndeg
        end if

      end do
c
c  and generate its rooted level structure.
c
400   continue

      call rootls ( root, xadj, adjncy, mask, nunlvl, xls, ls)
      if (nunlvl .le. nlvl) return
      nlvl = nunlvl
      if (nlvl .lt. ccsize ) go to 100
      return
      end
      subroutine fnspan ( xadj, adjncy, nodlvl, nspan, set,
     &  level, nadjs, adjs, leaf )

c*********************************************************************72
c
cc FNSPAN finds the span of a given subset
c  in a given level subgraph in a level structure.
c  the adjacent set of the span in the lower level is
c  also determined. if the span has an unnumbered node
c  in the higher level, an unnumbered leaf node (i.e. one
c  with no neighbor in next level) will be returned.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c  (xadj, adjncy) - the adjacent structure.
c  level - level number of the current set.
c
c  updated parameters 
c  (nspan, set) - the input set. on return, it contains
c  the resulting span set.
c  nodlvl - the level number vector. nodes considered
c  will have their nodlvl changed to zero.
c
c  output parameters 
c  (nadjs, adjs) - the adjacent set of the span in the
c  lower level.
c  leaf - if the span has an unnumbered higher level nodi
c  leaf returns an unnumbered leaf node in the le'
c  structure, otherwise, leaf is zero.
c
      integer adjncy(*)
      integer adjs(*)
      integer nodlvl(*), set(*)
      integer xadj(*), i, j, jstop, jstrt, leaf, level
      integer lvl, lvlm1, nadjs, nbr, nbrlvl, node
      integer nspan, setptr
c
c  initialization ...
c
      leaf = 0
      nadjs = 0
      setptr = 0

100   continue

      setptr = setptr + 1
      if ( setptr .gt. nspan ) then
        return
      end if
c
c  for each node in the partially spanned set ...
c
      node = set(setptr)
      jstrt = xadj(node)
      jstop = xadj(node + 1) - 1
      if ( jstop .lt. jstrt ) go to 100
c
c  for each neighbor of node, test its nodlvl value
c
      do 500 j = jstrt, jstop
        nbr = adjncy(j)
        nbrlvl = nodlvl(nbr)
        if (nbrlvl .le. 0) go to 500
        if (nbrlvl - level) 200, 300, 600
c
c  nbr is in level-1, add it to adjs.
c
200     nadjs = nadjs + 1
        adjs(nadjs) = nbr
        go to 400
c
c  nbr is in level, add it to the span set.
c
300     nspan = nspan + 1
        set(nspan) = nbr
400     nodlvl(nbr) = 0
500   continue
      go to 100
c
c  nbr is in level+1. find an unnumbered leaf node by
c  tracing a path up the level structure. then
c  reset the nodlvl values of nodes in adjs.
c
600   leaf = nbr
      lvl = level + 1
700   jstrt = xadj(leaf)
      jstop = xadj(leaf+1) - 1
      do 800 j = jstrt, jstop
        nbr = adjncy(j)
        if ( nodlvl(nbr) .le. lvl ) go to 800
          leaf = nbr
          lvl = lvl + 1
          go to 700
800   continue

      do i = 1, nadjs
        node = adjs(i)
        nodlvl(node) = level - 1
      end do

      return
      end
      subroutine fntadj ( xadj, adjncy, perm, invp,
     &  nblks, xblk, father, bnum )

c*********************************************************************72
c
cc FNTADJ determines the quotient tree adjacency structure of a graph. 
c
c  Discussion:
c
c    The structure is represented by the father vector.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c  (xadj, adjncy) - adjacency structure pair for the
c  (perm, invp) - the permutation vectors.
c  (nblks, xblk) - the tree partitioning. 
c
c  output parameters 
c  father - the father vector of the quotient tree.
c
c  working parameters 
c  bnum - temporary vector to store the block number
c  of each variable.
c
      integer adjncy(*)
      integer bnum(*)
      integer father(*)
      integer invp(*)
      integer perm(*), xblk(*)
      integer xadj(*), i, istop, istrt, j, jstop, jstrt
      integer k, nabor, nblks, nbm1, nbrblk, node
c
c  initialize the block number vector.
c
      do k = 1, nblks
        istrt = xblk(k)
        istop = xblk(k+1) - 1
        do i = istrt, istop
          node = perm(i)
          bnum(node) = k
        end do
      end do
c
c  for each block
c
      father(nblks) = 0
      nbm1 = nblks - 1

      do k = 1, nbm1

        istrt = xblk(k)
        istop = xblk(k+1) - 1
c
c  find its father block in the tree structure.
c
        do i = istrt, istop

          node = perm(i)
          jstrt = xadj(node)
          jstop = xadj(node+1) -1

	      do j = jstrt, jstop
            nabor = adjncy(j)
            nbrblk = bnum(nabor)
            if ( nbrblk .gt. k ) go to 500
          end do

          father(k) = 0

        end do

        go to 600
500     father(k) = nbrblk
600     continue

      end do

      return
      end
      subroutine fntenv ( xadj, adjncy, perm, invp,
     &  nblks, xblk, xenv, envsze )

c*********************************************************************72
c
cc FNTENV determines the envelope index
c  vector for the envelope of the diagonal blocks of a
c  tree partitioned system.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters -
c
c  (xadj, adjncy) - adjacency structure pair for the graph.
c  (perm, invp) - the permutation vectors.
c  (nblks, xblk) - the tree partitioning.
c
c  output parameters -
c  xenv - the envelope index vector.
c  envsze - the size of the envelope found.
c
      integer adjncy(*)
      integer invp(*)
      integer perm(*), xblk(*)
      integer xadj(*), xenv(*), blkbeg, blkend
      integer i, ifirst, j, jstop, jstrt, k, kfirst
      integer envsze, nblks, nbr, node

      envsze = 1
c
c  loop through each block in the partitioning .
c
      do 400 k = 1, nblks

        blkbeg = xblk(k)
        blkend = xblk(k+1) - 1
c
c  kfirst stores the first node in the k-th block
c  that has a neighbour in the previous blocks.
c
        kfirst = blkend

        do 300 i = blkbeg, blkend

          xenv(i) = envsze
          node = perm(i)
          jstrt = xadj(node)
          jstop = xadj(node+1) - 1
          if ( jstop .lt. jstrt ) go to 300
c
c  ifirst stores the first nonzero in the
c  i-th row within the k-th block.
c
          ifirst = i
          do j = jstrt, jstop
            nbr = adjncy(j)
            nbr = invp(nbr)
            if ( nbr .lt. blkbeg ) go to 100
            if ( nbr .lt. ifirst ) ifirst = nbr 
            go to 200
100         if ( kfirst .lt. ifirst ) ifirst = kfirst
            if (i .lt. kfirst ) kfirst = i
200         continue
          end do

          envsze = envsze + i - ifirst
300     continue
400   continue
      xenv(blkend+1) = envsze
      envsze = envsze - 1
      return
      end
      subroutine gen1wd ( neqns, xadj, adjncy, mask,
     &  nblks, xblk, perm, xls, ls )

c*********************************************************************72
c
cc GEN1WD finds a one-way dissection partitioning
c  for a general graph. fnlwd is used for each connected
c  component.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c    Input, integer NEQNS, the number of equations.
c
c  (xadj, adjncy) - the adjacency structure pair.
c
c  output parameters 
c  (nblks, xblk) - the partitioning found.
c  perm - the one-way dissection ordering.
c
c  working vectors 
c  mask - is used to mark variables that have
c  been numbered during the ordering process.
c  (xls, ls) - level structure used by rootls
c
      integer adjncy(*)
      integer ls(*)
      integer mask(*), perm(*)
      integer xblk(*), xls(*)
      integer xadj(*), ccsize, i, j, k, lnum
      integer nblks, neqns, nlvl, node, nsep
      integer num, root

      do i = 1, neqns
        mask(i) = 1
      end do

      nblks = 0
      num = 0

      do i = 1, neqns

        if ( mask(i) .eq. 0 ) go to 400
c
c  find a one-way dissector for each component.
c
        root = i
        call fn1wd ( root, xadj, adjncy, mask,
     &    nsep, perm(num+1), nlvl, xls, ls )
        num = num + nsep
        nblks = nblks + 1
        xblk(nblks) = neqns - num + 1
        ccsize = xls(nlvl+1) - 1
c
c  number the remaining nodes in the component.
c  each component in the remaining subgraph forms
c  a new block in the partitioning.
c
        do 300 j = 1, ccsize
          node = ls(j)
          if ( mask(node) .eq. 0 ) go to 300
          call rootls ( node, xadj, adjncy, mask,
     &      nlvl, xls, perm(num+1) )
          lnum = num + 1
          num = num + xls(nlvl+1) - 1
          nblks = nblks + 1
          xblk(nblks) = neqns - num + 1
          do k = lnum, num
            node = perm(k)
            mask(node) = 0
          end do
          if ( num .gt. neqns ) go to 500
300     continue

400     continue          

      end do
c
c  since dissectors found first should be ordered last,
c  routine revrse is called to adjust the ordering
c  vector, and the block index vector.
c
500   continue

      call i4vec_reverse ( neqns, perm )
      call i4vec_reverse ( nblks, xblk )
      xblk(nblks+1) = neqns + 1

      return
      end
      subroutine gennd ( neqns, xadj, adjncy, mask,
     &  perm, xls, ls )

c*********************************************************************72
c
cc GENND finds a nested dissection ordering for a general graph.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c    Input, integer NEQNS, the number of equations.
c
c  (xadj, adjncy) - adjacency structure pair.
c
c  output parameters 
c  perm - the nested dissection ordering.
c  
c  working parameters 
c  mask - is used to mask off variables that have
c  been numbered during the orderng process,
c  (xls, ls) - this level structure pair is used as
c  temporary storage by fnroot.
c
      integer adjncy(*)
      integer ls(*)
      integer mask(*), perm(*)
      integer xls(*)
      integer xadj(*), i, neqns, nsep, num, root

      do i = 1, neqns
        mask(i) = 1
      end do

      num = 0

      do 300 i = 1, neqns
c
c  for each masked component ...
c
200     if ( mask(i) .eq. 0 ) go to 300
        root = i
c
c  find a separator and number the nodes next,
c
        call fndsep ( root, xadj, adjncy, mask,
     &    nsep, perm(num+1), xls, ls )
        num = num + nsep
        if ( num .ge. neqns ) go to 400
        go to 200
300   continue
c
c  since separators found first should be ordered
c  last, routine revrse is called to adjust the
c  ordering vector,
c
400   continue

      call i4vec_reverse ( neqns, perm )

      return
      end      
      subroutine genqmd ( neqns, xadj, adjncy, perm, invp, deg,
     &  marker, rchset, nbrhd, qsize, qlink, nofsub )

c*********************************************************************72
c
cc GENQMD implements the minimum degre]
c  algorithm. it makes use of the implicit represei
c  ation of the elimination graphs by quotient grapi
c  and the notion of indistinguishable nodes.
c  caution - the adjacency vector adjncy will be
c  destroyed.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c    Input, integer NEQNS, the number of equations.
c
c  (xadj, adjncy) - the adjacency structure.
c
c  output parameters 
c  perm - the minimum degree ordering
c  invp - the inverse of perm.
c
c  jrking parameters 
c  deg - the degree vector. deg(i) is negative mean
c  node i has been numbered.
c  marker - a marker vector, where marker(i) is
c  negative means node i has been merged wit
c  another node and thus can be ignored.
c  rchset - vector used for the reachable set.
c  nbrhd - vector used for the neighborhood set.
c  qsize - vector used to store the size of
c  indistinguishable supernodes.
c  qlink - vector to store indistinguishable nodes,
c  i, qlink(i), qlink(qlink(i)) ... are the
c  members of the supernode represented by i.
c
      integer adjncy(*)
      integer invp(*)
      integer perm(*), deg(*), marker(*)
      integer rchset(*), nbrhd(*), qsize(*), qlink(*)
      integer xadj(*), inode, ip, irch, j, mindeg, ndeg
      integer neqns, nhdsze, node, nofsub, np, num, nump1
      integer nxnode, rchsze, search, thresh
c
c  initialize degree vector and other working variab
c
      mindeg = neqns
      nofsub = 0
      
      do node = 1, neqns
        perm(node) = node
        invp(node) = node
        marker(node) = 0
        qsize(node) = 1
        qlink(node) = 0
        ndeg = xadj(node+1) - xadj(node)
        deg(node) = ndeg
        if ( ndeg .lt. mindeg ) mindeg = ndeg
      end do

      num = 0
c
c  perform threshold search to get a node of min degree.
c  variable search points to where search should start.
c
200   continue

      search = 1
      thresh = mindeg
      mindeg = neqns

300   continue

      nump1 = num + 1
      if ( nump1 .gt. search ) search = nump1
      do 400 j = search, neqns
        node = perm(j)
        if ( marker(node) .lt. 0 ) go to 400
        ndeg = deg(node)
        if ( ndeg .le. thresh ) go to 500
        if ( ndeg .lt. mindeg ) mindeg = ndeg
400   continue

      go to 200
c
c  node has minimum degree. find its reachable sets by
c  calling qmdrch.
c
500   continue

      search = j
      nofsub = nofsub + deg(node)
      marker ( node ) = 1
      call qmdrch (node, xadj, adjncy, deg, marker,
     &  rchsze, rchset, nhdsze, nbrhd )
c
c  eliminate all nodes indistinguishable from node.
c  they are given by node, qlink(node), ....
c
      nxnode = node
600   continue
      num = num + 1
      np = invp(nxnode)
      ip = perm(num)
      perm(np) = ip
      invp(ip) = np
      perm(num) = nxnode
      invp(nxnode) = num
      deg(nxnode) = - 1
      nxnode = qlink(nxnode)
      if (nxnode .gt. 0) go to 600

      if ( rchsze .le. 0 ) go to 800
c
c  update the degrees of the nodes in the reachable
c  set and identify indistinguishable nodes.
c
      call qmdupd ( xadj, adjncy, rchsze, rchset, deg,
     &  qsize, qlink, marker,
     &  rchset(rchsze+1), nbrhd(nhdsze+1) )
c
c  reset marker value of nodes in reach set.
c  update threshold value for cyclic search.
c  also call qmdqt to form new quotient graph.
c
      marker(node) = 0

      do irch = 1, rchsze
        inode = rchset(irch)
        if ( marker(inode) .lt. 0 ) go to 700
        marker(inode) = 0
        ndeg = deg(inode)
        if ( ndeg .lt. mindeg ) mindeg = ndeg
        if ( ndeg .gt. thresh ) goto 700
        mindeg = thresh
        thresh = ndeg
        search = invp(inode)
700     continue
      end do

      if ( nhdsze .gt. 0 ) then
        call qmdqt ( node, xadj,
     &    adjncy, marker, rchsze, rchset, nbrhd )
      endif

800   if ( num .lt. neqns ) go to 300

      return
      end
      subroutine genrcm ( neqns, xadj, adjncy, perm, mask, xls )

c*********************************************************************72
c
cc GENRCM finds the reverse Cuthill-Mckee ordering for a general graph.
c
c  Discussion:
c
c    For each connected component in the graph, GENRCM obtains the order
c    by calling the subroutine RCM.
c
c  Modified:
c
c    31 December 2008
c
c  Parameters:
c
c    Input, integer NEQNS, the number of equations.
c
c    Input, (xadj, adjncy) - array pair containing the adjacency
c    structure of the graph of the matrix.
c
c    Output, integer PERM(NEQNS), the RCM ordering.
c
c    Workspace, integer MASK(NEQNS), used to mark variables that have been
c    numbered during the ordering process.  It is initialized to 1, and set 
c    to zero as each item is numbered.
c
c    Workspace, xls - the index vector for a level structure. this
c    level structure is stored in the currently
c    unused spaces in the permutation vector perm.
c
      implicit none

      integer neqns

      integer adjncy(*)
      integer ccsize
      integer i
      integer mask(neqns)
      integer nlvl
      integer num
      integer perm(neqns)
      integer root
      integer xls(*)
      integer xadj(*)

      do i = 1, neqns
        mask(i) = 1
      end do

      num = 1

      do i = 1, neqns
c
c  for each masked connected component ...
c
        if ( mask(i) .ne. 0 ) then

          root = i
c
c  first find a pseudo-peripheral node root.
c
c  note that the level structure found by
c  fnroot is stored starting at perm(num).
c
          call fnroot ( root, xadj, adjncy, mask,
     &      nlvl, xls, perm(num) )
c
c  RCM is called to order the component
c  using root as the starting node.
c
          call rcm ( root, xadj, adjncy, mask,
     &      perm(num), ccsize, xls )

          num = num + ccsize

          if ( neqns .lt. num ) then
            return
          end if

        end if

      end do

      return
      end
      subroutine genrqt ( neqns, xadj, adjncy, nblks,
     &  xblk, perm, xls, ls, nodlvl )

c*********************************************************************72
c
cc GENRQT determines a
c  partitioned ordering for a possibly disconnected
c  graph using the refined quotient tree algorithm.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters
c
c    Input, integer NEQNS, the number of equations.
c
c  (xadj, adjncy) - the adjacency structure.
c
c	output parameters 
c	(nblks, xblk) - the quotient tree partitioning.
c	perm - the permutation vector.
c
c	working parameters 
c	(xls, ls) - this level structure pair is used by
c	fnroot to find a pseudo-peripheral node.
c	nodlvl - a temporary vector to store the level
c	number of each node in a level structure.
c
      integer adjncy(*)
      integer ls(*), nodlvl(*), perm(*)
      integer xblk(*), xls(*)
      integer xadj(*), i, ixls, leaf, nblks, neqns, nlvl
      integer root
c
c  initialization.
c
      do i = 1, neqns
        nodlvl(i) = 1
      end do

      nblks = 0
      xblk(1) = 1
c
c  for each connected component, find a rooted level
c  structure, and then call rqtree for its block order.
c
      do i = 1, neqns

        if ( 0 .lt. nodlvl(i) ) then
          root = i
          call fnlvls ( root, xadj, adjncy, nodlvl, nlvl, xls, ls )
          ixls = xls(nlvl)
          leaf = ls(ixls)
          call rqtree ( leaf, xadj, adjncy, perm,
     &      nblks, xblk, nodlvl, xls, ls )
         end if

      end do

      return
      end
      subroutine gsfct ( neqns, xlnz, lnz, xnzsub, nzsub, diag,
     &  link, first, temp, flag )

c*********************************************************************72
c
cc GSFCT performs the symmetric
c  factorization for a general sparse system, stor
c  the compressed subscript data format.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c    Input, integer NEQNS, the number of equations.
c
c  xlnz - index vector for lnz. xlnz(i) points to the
c  start of nonzeros in column i of factor L.
c
c  (xnzsub, nzsub) - the compressed subscript data
c  structure for factor l.
c
c  updated parameters 
c
c  Input/output, LNZ(*). on input, contains nonzeros of a, and on
c  return, the nonzeros of l.
c
c  Input/output, double precision DIAG(*).  the diagonal of l overwrites 
c  that of a.
c
c  Output, integer FLAG, the error flag. it is set to 1 if a zero
c  or negative square root occurs during the factorization.
c
c  parameters 
c  LINK, at step j. the list in
c  link(j), link(link(j)), ...........
c  consists of those columns that will modify
c  the column l(*,j).
c
c  FIRST - temporary vector to point to the first
c  nonzero in each column that will be used
c  next for modification.
c
c  TEMP - a temporary vector to accumulate modifica'
c
      double precision diag(*)
      double precision lnz(*), temp(*), diagj, ljk
      integer link(*), nzsub(*)
      integer first(*), xlnz(*), xnzsub(*)
      integer i, flag, ii, istop, istrt, isub, j
      integer k, kfirst, neqns, newk

      flag = 0
c
c  Initialize working vectors.
c
      do i = 1, neqns
        link(i) = 0
      end do

      do i = 1, neqns
        temp(i) = 0.0D+00
      end do
c
c  compute column l(*,j) for j - 1,...., neqns.
c
      do 600 j = 1, neqns
c
c  for each column l(*,k) that affects l(s,j).
c
        diagj = 0.0D+00
        newk = link(j)
200     k = newk
        if (k.eq.0) go to 400
        newk = link(k)
c
c  outer product modification of l($ j) by
c  l(e,k) starting at first(k) of l(;,k).
c
        kfirst = first(k)
        ljk = lnz(kfirst)
        diagj = diagj + ljk*ljk
        istrt = kfirst + 1
        istop = xlnz(k+1) - 1
        if ( istop .lt. istrt ) go to 200
c
c  before modification, update vectors first,
c  and link for future modification steps.
c
        first(k) = istrt
        i = xnzsub(k) + (kfirst-xlnz(k)) + 1
        isub = nzsub(i)
        link(k) = link(isub)
        link(isub) = k
c
c  the actual mod is saved in vector temp.
c
        do ii = istrt, istop
          isub = nzsub(i)
          temp(isub) = temp(isub) + lnz(ii)*ljk
          i = i + 1
        end do

        go to 200
c
c  apply the modifications accumulated in temp to column l(*,j).
c
400     continue

        diagj = diag(j) - diagj

        if ( diagj .le. 0.0D+00 ) then
          flag = 1
          return
        end if

        diagj = sqrt(diagj)
        diag(j) = diagj
        istrt = xlnz(j)
        istop = xlnz(j+1) - 1
        if ( istop .lt. istrt ) go to 600
        first(j) = istrt
        i = xnzsub(j)
        isub = nzsub(i)
        link(j) = link(isub)
        link(isub) = j

        do ii = istrt, istop
          isub = nzsub(i)
          lnz(ii) = ( lnz(ii) - temp(isub) ) / diagj
          temp(isub) = 0.0D+00
          i = i + 1
        end do

600   continue

      return
      end
      subroutine gsslv ( neqns, xlnz, lnz, xnzsub, nzsub,
     &  diag, rhs )

c*********************************************************************72
c
cc GSSLV solves a factored system in compressed subscript sparse format.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c    Input, integer NEQNS, the number of equations.
c
c  {xlnz, lnz) - structure of nonzeros in l.
c  (xnzsub, nzsub) - compressed subscript structure.
c  diag - diagonal components of l.
c
c  updated parameter. 
c
c  rhs - on input, it contains the rhs vector, and on
c  output, the solution vector.
c
      double precision diag(*)
      integer i
      integer ii
      integer istop
      integer istrt
      integer isub
      integer j
      integer jj
      double precision lnz(*)
      integer neqns
      integer nzsub(*)
      double precision rhs(*)
      double precision rhsj
      double precision s
      integer xlnz(*)
      integer xnzsub(*)
c
c  Forward substitution.
c
      do j = 1, neqns
      
        rhsj = rhs(j) / diag(j)
        rhs(j) = rhsj
        istrt = xlnz(j)
        istop = xlnz(j+1) - 1
        i = xnzsub(j)

        do ii = istrt, istop
          isub = nzsub(i)
          rhs(isub) = rhs(isub) - lnz(ii) * rhsj
          i = i + 1
        end do

      end do
c
c  Backward substitution.
c
      j = neqns

      do jj = 1, neqns

        s = rhs(j)
        istrt = xlnz(j)
        istop = xlnz(j+1) - 1

        i = xnzsub(j)

        do ii = istrt, istop
          isub = nzsub(i)
          s = s - lnz(ii) * rhs(isub)
          i = i + 1
        end do

        rhs(j) = s / diag(j)
        j = j - 1

      end do

      return
      end
      subroutine i4vec_reverse ( n, a )

c*********************************************************************72
c
cc I4VEC_REVERSE reverses the elements of an I4VEC.
c
c  Discussion:
c
c    An I4VEC is a vector of I4's.
c
c    In FORTRAN90, call I4VEC_REVERSE is equivalent to:
c
c      A(1:N) = A(N:1:-1)
c
c  Example:
c
c    Input:
c
c      N = 5,
c      A = ( 11, 12, 13, 14, 15 ).
c
c    Output:
c
c      A = ( 15, 14, 13, 12, 11 ).
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    25 June 2009
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, integer N, the number of entries in the array.
c
c    Input/output, integer A(N), the array to be reversed.
c
      implicit none

      integer n

      integer a(n)
      integer i
      integer t

      do i = 1, n / 2
        t        = a(i)
        a(i)     = a(n+1-i)
        a(n+1-i) = t
      end do

      return
      end
      subroutine perm_inverse ( neqn, perm, invprm ) 

c*********************************************************************72
c
cc PERM_INVERSE produces the inverse permutation.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    01 January 2009
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, integer NEQN, the number of equations.
c
c    Input, integer PERM(NEQN), the reordering of the variables and equations.
c
c    Output, integer INVPRM(NEQN), the inverse ordering, with the property 
c    that INVPRM(PERM(K))=K.
c
      integer neqn

      integer i
      integer invprm(neqn)
      integer k
      integer perm(neqn)

      do i = 1, neqn
        k = perm(i)
        invprm(k) = i
      end do

      return
      end
      subroutine permrv ( neqn, rhs, perm )

c*********************************************************************72
c
cc PERMRV should be called once the linear system has been solved and the
c  solution returned in rhs.  permrv then undoes the permutation of rhs, 
c  restoring the original ordering.  to do this, it needs the perm 
c  vector which defined the reordering used by the solver.
c
c  Modified:
c
c    01 January 2009
c
c  Parameters:
c
c    Input, integer NEQN, the number of equations.
c 
c  rhs    input/output, double precision rhs(neqn).  
c         on input, rhs contains the solution of the permuted linear 
c         system.
c         on output, rhs contains the solution of the original linear 
c         system.
c 
c  perm   input, integer perm(neqn), the permutation information.    
c         perm(i)=k means that the k-th equation and variable in the 
c         original ordering became the i-th equation and variable in the 
c         reordering.
c
      integer neqn

      integer i
      integer iput
      integer istart
      integer perm(neqn)
      double precision pull
      double precision put
      double precision rhs(neqn)
c
c  mark perm with negative signs which will be removed
c  as each permuted element is restored to its rightful place
c
      do i = 1, neqn
        perm(i) = - perm(i)
      enddo
c
c  search for the next element of perm which is the first
c  element of a permutation cycle
c
      istart = 0

   20 continue

      istart = istart + 1
      if ( istart .gt. neqn ) then
        return
      end if
      if(perm(istart).gt.0)go to 20
      if(iabs(perm(istart)).ne.istart)go to 30
      perm(istart)=iabs(perm(istart))
      go to 20
c
c  begin a cycle
c
   30 continue

      perm(istart)=iabs(perm(istart))
      iput=istart
      pull=rhs(iput)

   40 continue

      iput=iabs(perm(iput))
      put=rhs(iput)
      rhs(iput)=pull
      pull=put
      if(perm(iput).gt.0)go to 20
      perm(iput)=iabs(perm(iput))
      go to 40

      end
      subroutine qmdmrg ( xadj, adjncy, deg, qsize, qlink,
     &  marker, dego, nhdsze, nbrhd, rchset, ovrlp)

c*********************************************************************72
c
cc QMDMRG merges indistinguishable nodes in
c  the minimum degree ordering algorithm.
c  it also computes the new degrees of these
c  new supernodes.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c  (xadj, adjncy) - the adjacency structure.
c  deg0 - the number of nodes in the given set.
c  (nhdsze, nbrhd) - the set of eliminated supernodes
c  adjacent to some nodes in the set.
c
c  updated parameters 
c  deg - the degree vector.
c  qsize - size of indistinguishable nodes.
c  qlink - linked list for indistinguishable nodes.
c  marker - the given set is given by those nodes with
c  marker value set to 1. those nodes with degree
c  updated will have marker value set to 2.
c
c  working parameters 
c  rchset - the reachable set.
c  ovrlp - temp vector to store the intersection of two
c  reachable sets.
c
      integer adjncy(*), deg(*), qsize(*), qlink(*)
      integer marker(*), rchset(*), nbrhd(*), ovrlp(*)
      integer xadj(*), dego, deg1, head, inhd, iov, irch
      integer j, jstrt, jstop, link, lnode, mark, mrgsze
      integer nabor, nhdsze, node, novrlp, rchsze, root
c
c  initialization ...
c
      if ( nhdsze .le. 0 ) return

      do inhd = 1, nhdsze
        root = nbrhd(inhd)
        marker(root) = 0
      end do
c
c  loop through each eliminated supernode in the
c  (nhdsze, nbrhd).
c
      do 1400 inhd = 1, nhdsze
        root = nbrhd(inhd)
        marker (root) = - 1
        rchsze = 0
        novrlp = 0
        deg1 = 0
200     jstrt = xadj(root)
        jstop = xadj(root+1) - 1
c
c  determine the reachable set and its inters
c  ion with the input reachable set.
c
        do 600 j = jstrt, jstop
          nabor = adjncy(j)
          root = - nabor
          if (nabor) 200, 700, 300
300       mark = marker(nabor)
          if ( mark ) 600, 400, 500
400       rchsze = rchsze + 1
          rchset(rchsze) = nabor
          deg1 = deg1 + qsize(nabor)
          marker(nabor) = 1
          go to 600
500       if ( mark .gt. 1 ) goto 600
          novrlp = novrlp + 1
          ovrlp(novrlp) = nabor
          marker(nabor) = 2
600     continue
c
c  from the overlapped set, determine the node
c  that can be merged together.
c
700     head = 0
        mrgsze = 0
        do 1100 iov = 1, novrlp
          node = ovrlp(iov)
          jstrt = xadj(node)
          jstop = xadj(node+1) - 1
          do j = jstrt, jstop
            nabor = adjncy(j)
            if ( marker(nabor) .eq. 0 ) then
              marker(node) = 1
              go to 1100
            end if
          end do
c
c  node belongs to the new merged supernode.
c  update the vectors qlink and qsize.
c
          mrgsze = mrgsze + qsize(node)
          marker(node) = - 1
          lnode = node
900       link = qlink(lnode)
          if ( link .le. 0 ) goto 1000
          lnode = link
          go to 900
1000      qlink(lnode) = head
          head = node
1100    continue

        if ( head .le. 0 ) goto 1200
        qsize(head) = mrgsze
        deg(head) = dego + deg1 - 1
        marker (head) = 2
c
c  reset marker values.
c
1200    root = nbrhd(inhd)
        marker(root) = 0
        
        do irch = 1, rchsze
          node = rchset(irch)
          marker(node) = 0
        end do

1400  continue

      return
      end
      subroutine qmdqt ( root, xadj, adjncy, marker,
     &  rchsze, rchset, nbrhd )

c*********************************************************************72
c
cc QMDQT performs the quotient
c  transformation after a node has been elimin}
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c  root - the node just eliminated. it becomes
c  representative of the new supernode.
c  (xadj, adjncy) - the adjacency structure.
c  (rchsze, rchset) - the reachable set of root
c  old quotient graph.
c  nbrhd - the neighborhood set which will be b
c  with root to form the new supernode.
c  marker - the marker vector.
c
c  updated parameter 
c  adjncy - becomes the adjncy of the quotient
c
      integer adjncy(*)
      integer marker(*), rchset(*), nbrhd(*)
      integer xadj(*), inhd, irch, j, jstrt, jstop, link
      integer nabor, node, rchsze, root

      irch = 0
      inhd = 0
      node = root
100   jstrt = xadj(node)
      jstop = xadj(node+1) - 2
      if ( jstop .lt. jstrt ) go to 300
c
c  place reach nodes into the adjacent list 4
c
      do j = jstrt, jstop
        irch = irch + 1
        adjncy(j) = rchset(irch)
        if ( irch .ge. rchsze ) go to 400
      end do
c
c  link to other space provided by the neighborhood set.
c
300   link = adjncy(jstop+1)
      node = - link
      if ( link .lt. 0 ) go to 100
      inhd = inhd + 1
      node = nbrhd(inhd)
      adjncy(jstop+1) = - node
      go to 100
c
c  all reachable nodes have been saved. end the adj list.
c  add root to the nbr list of each node in the reach set.
c
400   adjncy(j+1) = 0
      do 600 irch = 1, rchsze
        node = rchset(irch)
        if ( marker(node) .lt. 0 ) go to 600
        jstrt = xadj(node)
        jstop = xadj(node+1) - 1
        do 500 j = jstrt, jstop
          nabor = adjncy(j)
          if ( marker(nabor) .ge. 0 ) go to 500
          adjncy(j) = root
          go to 600
500     continue
600   continue

      return
      end
      subroutine qmdrch ( root, xadj, adjncy, deg, marker,
     &  rchsze, rchset, nhdsze, nbrhd )

c*********************************************************************72
c
cc QMDRCH determines the reachable set
c  a node through a given subset. the adjacency stre
c  is assumed to be stored in a quotient graph forma1
c
c  Modified:
c
c    01 January 2009
c
c input parameters 
c
c  root - the given node not in the subset.
c  (xadj, adjncy) - the adjacency structure pair.
c  deg - the degree vector. deg( i ) lt 0 means the n(
c  belongs to the given subset.
c
c  output parameters 
c  (rchsze, rchset) - the reachable set.
c  (nhdsze, nbrhd) - the neighborhood set.
c
c  updated parameters 
c  marker - the marker vector for reach and nbrhd se1
c  gt 0 means the node is in reach set.
c  lt 0 means the node has been merged with
c  others in the quotient or it is in nbrhd sl
c
      integer adjncy(*)
      integer deg(*), marker(*)
      integer rchset(*), nbrhd(*)
      integer xadj(*), i, istrt, istop, j, jstrt, jstop
      integer nabor, nhdsze, node, rchsze, root
c
c  loop through the neighbors of root in the quotient graph.
c
      nhdsze = 0
      rchsze = 0
      istrt = xadj(root)
      istop = xadj(root+1) - 1
      if ( istop .lt. istrt ) return
      
      do 600 i = istrt, istop
        nabor = adjncy(i)
        if ( nabor .eq. 0 ) return
        if ( marker(nabor) .ne. 0 ) go to 600
        if ( deg(nabor) .lt. 0 ) go to 200
c
c  include nabor into the reachable set.
c
        rchsze = rchsze + 1
        rchset(rchsze) = nabor
        marker(nabor) = 1
        go to 600
c
c  nabor has been eliminated. find nodes
c  reachable from it.
c
200     marker(nabor) = -1
        nhdsze = nhdsze + 1
        nbrhd(nhdsze) = nabor
300     jstrt = xadj(nabor)
        jstop = xadj(nabor+1) - 1
        
        do 500 j = jstrt, jstop
          node = adjncy(j)
          nabor = - node
          if (node) 300, 600, 400
400       if ( marker(node) .ne. 0 ) go to 500
          rchsze = rchsze + 1
          rchset(rchsze) = node
          marker(node) = 1
500     continue
       
600   continue

      return
      end
      subroutine qmdupd ( xadj, adjncy, nlist, list, deg,
     &  qsize, qlink, marker, rchset, nbrhd)

c*********************************************************************72
c
cc QMDUPD performs degree update for a set
c  of nodes in the minimum degree algorithm.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c  (xadj, adjncy) - the adjacency structure.
c  (nlist, list) - the list of nodes whose degree has to
c  be updated.
c
c  updated parameters 
c  deg - the degree vector.
c  qsize - size of indistinguishable supernodes.
c  qlink - linked list for indistinguishable nodes.
c  marker - used to mark those nodes in reach/nbrhd sets.
c
c  working parameters 
c  rchset - the reachable set.
c  nbrhd - the neighborhood set.
c
      integer adjncy(*)
      integer list(*), deg(*), marker(*)
      integer rchset(*), nbrhd(*), qsize(*), qlink(*)
      integer xadj(*), deg0, deg1, il, inhd, inode, irch
      integer j, jstrt, jstop, mark, nabor, nhdsze, nlist
      integer node, rchsze 
c
c  find all eliminated supernodes that are adjm
c  to some nodes in the given list. put them in'
c  (nhdsze, nbrhd). deg0 contains the number of
c  nodes in the list.
c
      if ( nlist .le. 0 ) then
        return
      end if

      deg0 = 0
      nhdsze = 0
      do 200 il = 1, nlist
        node = list(il)
        deg0 = deg0 + qsize(node)
        jstrt = xadj(node)
        jstop = xadj(node+1) - 1
        do 100 j = jstrt, jstop
          nabor = adjncy(j)
          if ( marker(nabor) . ne . 0 . or .
     &      deg(nabor) .ge. 0 ) go to 100
          marker(nabor) = - 1
          nhdsze = nhdsze + 1
          nbrhd(nhdsze) = nabor
100     continue
200   continue
c
c  merge indistinguishable nodes in the list by
c  calling the subroutine qmdmrg.
c
      if ( nhdsze .gt. 0 ) then
        call qmdmrg ( xadj, adjncy, deg, qsize, qlink,
     &    marker, deg0, nhdsze, nbrhd, rchset,
     &    nbrhd(nhdsze+1) )
      endif
c
c  find the new degrees of the nodes that have not been
c  merged.
c
      do 600 il = 1, nlist

        node = list(il)
        mark = marker(node)

        if ( mark .gt. 1 .or. mark .lt. 0 ) go to 600

        marker(node) = 2

        call qmdrch ( node, xadj, adjncy, deg, marker,
     &    rchsze, rchset, nhdsze, nbrhd)

        deg1 = deg0

        do irch = 1, rchsze
          inode = rchset(irch)
          deg1 = deg1 + qsize(inode)
          marker(inode) = 0
        end do

        deg(node) = deg1 - 1

        do inhd = 1, nhdsze
          inode = nbrhd(inhd)
          marker(inode) = 0
        end do

600   continue

      return
      end
      subroutine rcm ( root, xadj, adjncy, mask, perm, ccsize, deg )

c*********************************************************************72
c
cc RCM renumbers a connected component by the reverse Cuthill McKee algorithm.
c
c  Discussion:
c
c    The connected component is specified by mask and root.
c    The numbering is done starting at the node ROOT, and proceeding
c    using the Reverse Cuthill-McKee algorithm.
c
c    An outline of the algorithm is as follows:
c
c    X(1) = ROOT.
c
c    for ( I = 1 to N-1)
c      Find all unlabeled neighbors of X(I), 
c      assign them the next available labels, in order of increasing degree.
c
c    When done, reverse the ordering.
c
c  Modified:
c
c    31 December 2008
c
c  Author:
c
c    Alan George, Joseph Liu
c
c  Reference:
c
c    Alan George, Joseph Liu,
c    Computer Solution of Large Sparse Positive Definite Systems,
c    Prentice Hall, 1981,
c    ISBN: 0131652745,
c    LC: QA188.G46.
c
c  input parameters 
c
c    Input, integer ROOT, the node that defines the connected component.
c    It is used as the starting point for the RCM ordering.
c
c    Input, integer XADJ(N+1).  Information about row I is stored
c    in entries XADJ(I) through XADJ(I+1)-1 of ADJNCY.
c
c    Input, integer ADJNCY(*), the adjacency structure. 
c    For each row, the column indices of the nonzero entries.
c
c    Input/output, integer MASK(N), a mask for the nodes.  Only those nodes 
c    with nonzero input mask values are considered by the routine.  The
c    nodes numbered by RCM will have their mask values set to zero.
c
c    Output, integer PERM(N), the RCM ordering.
c
c    Output, integer CCSZE, the size of the connected component
c    that has been numbered.
c
c    Workspace, integer DEG(N), a temporary vector used to hold the degree
c    of the nodes in the section graph specified by MASK and ROOT.
c
      integer adjncy(*)
      integer ccsize
      integer deg(*)
      integer fnbr
      integer i
      integer j
      integer jstop
      integer jstrt
      integer k
      integer l
      integer lbegin
      integer lnbr
      integer lperm
      integer lvlend
      integer mask(*)
      integer nbr
      integer node
      integer perm(*)
      integer root
      integer xadj(*)
c
c  Find the degrees of the nodes in the component specified by MASK and ROOT.
c
      call degree ( root, xadj, adjncy, mask, deg, ccsize, perm )

      mask(root) = 0

      if ( ccsize .le. 1 ) then
        return
      end if

      lvlend = 0
      lnbr = 1
c
c  lbegin and lvlend point to the beginning and
c  the end of the current level respectively.
c
100   continue

      lbegin = lvlend + 1
      lvlend = lnbr

      do i = lbegin, lvlend
c
c  For each node in current level.
c
        node = perm(i)
        jstrt = xadj(node)
        jstop = xadj(node+i) - 1
c
c  Find the unnumbered neighbors of node.
c
c  fnbr and lnbr point to the first and last
c  nundneumibnred neighbors respectively of the curr
c
        fnbr = lnbr + 1

        do j = jstrt, jstop
          nbr = adjncy(j)
          if ( mask(nbr) .ne. 0 ) then
            lnbr = lnbr + 1
            mask(nbr) = 0
            perm(lnbr) = nbr
          end if
        end do
c
c  sort the neighbors of node in increasing
c  order by degree. linear insertion is used
c
        if ( fnbr .lt. lnbr ) then

          k = fnbr

300       continue

          l = k
          k = k + 1
          nbr = perm(k)

400       if ( l .lt. fnbr ) go to 500
          lperm = perm(l)
          if(deg(lperm) .le. deg(nbr) ) go to 500
          perm(l+1)=lperm
          l=l-1
          go to 400

500       continue

          perm(l+1)=nbr
 
          if ( k .lt. lnbr ) then
            go to 300
          end if

        end if

      end do

      if ( lnbr .gt. lvlend ) then
        go to 100
      end if
c
c  we now have the cuthill mckee ordering.
c  reverse it.
c
      k = ccsize / 2
      l = ccsize

      do i=1, k
        lperm = perm(l)
        perm(l) = perm(i)
        perm(i) = lperm
        l = l - 1
      end do

      return
      end
      subroutine reach ( root, xadj, adjncy, smask, marker,
     &  rchsze, rchset, nhdsze, nbrhd )

c*********************************************************************72
c
cc REACH determines the reachable set of a node.
c
c  Discussion:
c
c    The routine determines the reachable set of a node y through a subset s
c    (i.e. reach(y,s) ) in a given subgraph. moreover
c    it returns the neighborhood set of y in s. i.e
c    nbrhd(y,s), the set of nodes in s that can be
c    reached from y through a subset of s.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters -
c
c    Input, integer ROOT, the given node not in the subset s.
c
c  (xadj, adjncy) - the adjacency structure pair.
c  smask - the mask vector for the set s.
c  - 0, if the node is not in s.
c  > 0, if the node is in s.
c
c  output parameters -
c  (nhdsze, nbrhd) - the neighborhood set.
c  (rchsze, rchset) - the reachable set.
c
c  updated parameters -
c  marker - the marker vector used to define the subgraph
c  nodes in the subgraph have marker value 0.
c  on return, the reachable and neighborhood node
c  sets have their marker values reset to root.
c
      integer adjncy(*)
      integer marker(*)
      integer nbrhd(*)
      integer rchset(*)
      integer smask(*)
      integer xadj(*), i, istop, istrt, j, jstop, jstrt
      integer nabor, nbr, nhdbeg, nhdptr, nhdsze, node
      integer rchsze, root
c
c  initialization ...
c
      nhdsze = 0
      rchsze = 0

      if ( marker(root) .le. 0 ) then
        rchsze = 1
        rchset(1) = root
        marker(root) = root
      end if

      istrt = xadj(root)
      istop = xadj(root+1) - 1
c
c  loop through the neighbors of root ...
c
      do 600 i = istrt, istop

        nabor = adjncy(i)
        if ( marker(nabor) .ne. 0 ) go to 600
        if ( smask(nabor) .gt. 0 ) go to 200
c
c  nabor is not in s. include it in the reach set.
c
        rchsze = rchsze + 1
        rchset(rchsze) = nabor
        marker(nabor) = root
        go to 600
c
c  nabor is in subset s. and has not been considerei
c  include it into the nbrhd set and find the nodes
c  reachable from root through this nabor.
c
200     nhdsze = nhdsze + 1
        nbrhd(nhdsze) = nabor
        marker(nabor) = root
        nhdbeg = nhdsze
        nhdptr = nhdsze
300     node = nbrhd(nhdptr)
        jstrt = xadj(node)
        jstop = xadj(node+1) - 1
        do 500 j = jstrt, jstop
          nbr = adjncy(j)
          if ( marker(nbr) .ne. 0 ) go to 500
          if ( smask(nbr) .eq. 0 ) go to 400
          nhdsze = nhdsze + 1
          nbrhd(nhdsze) = nbr
          marker(nbr) = root
          go to 500
400       rchsze = rchsze + 1
          rchset(rchsze) = nbr
          marker(nbr) = root
500     continue
        nhdptr = nhdptr + 1
        if ( nhdptr .le. nhdsze ) go to 300

600   continue

      return
      end
      subroutine rootls ( root, xadj, adjncy, mask, nlvl, xls, ls )

c*********************************************************************72
c
cc ROOTLS generates the level structure rooted at a given node.
c
c  Discussion:
c
c    Only those nodes for which MASK is nonzero will be considered.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c  root, the node at which the level structure is to be rooted.
c
c  (xadj, adjncy) - adjacency structure pair for the
c  given graph.
c
c  mask - is used to specify a section subgraph. nodes
c  with mask(i)=0 are ignored.
c
c  output parameters -
c  nlvl - is the number of levels in the level structure.
c  (xls, ls) - array pair for the rooted level structure.
c
      integer adjncy(*)
      integer ls(*), mask(*), xls(*)
      integer xadj(*), i, j, jstop, jstrt, lbegin
      integer ccsize, lvlend, lvsize, nbr, nlvl
      integer node, root
c
c  Initialization.
c
      mask(root) = 0
      ls(1) = root
      nlvl = 0
      lvlend = 0
      ccsize = 1
c
c  LBEGIN is the pointer to the beginning of the current
c  level, and lvlend points to the end of this level.
c
200   lbegin = lvlend + 1
      lvlend = ccsize
      nlvl = nlvl + 1
      xls(nlvl) = lbegin
c
c  generate the next level by finding all the masked
c  neighbors of nodes in the current level.
c
      do 400 i = lbegin, lvlend
        node = ls(i)
        jstrt = xadj(node)
        jstop = xadj(node + 1) - 1
        if ( jstop .lt. jstrt ) go to 400
        do 300 j = jstrt, jstop
          nbr = adjncy(j)
          if (mask(nbr) .eq. 0) go to 300
          ccsize = ccsize + 1
          ls(ccsize) = nbr
          mask(nbr) = 0
300     continue
400   continue
c
c  compute the current level width.
c  if it is nonzero, generate the next level.
c
      lvsize = ccsize - lvlend
      if (lvsize .gt. 0 ) go to 200
c
c  reset mask to one for the nodes in the level structure.
c
      xls(nlvl+1) = lvlend + 1
      do i = 1, ccsize
        node = ls(i)
        mask(node) = 1
      end do

      return
      end
      subroutine rqtree ( leaf, xadj, adjncy, perm,
     &  nblks, xblk, nodlvl, adjs, stack )

c*********************************************************************72
c
cc RQTREE finds a quotient tree ordering
c  for the component specified by leaf and nodlvl.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c  (xadj, adjncy) - the adjacency structure.
c
c  leaf - the input node that defines the connected
c  component. it is also a leaf node in the
c  rooted level structure passed to rqtree.
c  i.e. it has no neighbor in the next level.
c
c  output parameters 
c  perm - the permutation vector containing the ordering.
c  (nblks, xblk) - the quotient tree partitioning.
c
c  updated parameters 
c  nodlvl - the node level number vector. nodes in the
c  component have their nodlvl set to zero as
c  as they are nwbered.
c
c  working parameters 
c  adjs - temporary vector to store the adjacent set
c  of nodes in a particular level.
c  stack - temporary vector used to maintain the stack
c  of node subsets. it is organised as 
c  ( subset nodes, subset size, subset level )
c
      integer adjncy(*)
      integer adjs(*)
      integer nodlvl(*)
      integer perm(*)
      integer stack(*), xblk(*)
      integer xadj(*), blksze, ip, j, jp, leaf, level
      integer nadjs, nblks, node, npop, nuleaf
      integer num, toplvl, topstk
c
c  initialize the stack vector and its pointers.
c
      stack(1) = 0
      stack(2) = 0
      topstk = 2
      toplvl = 0
      num = xblk(nblks+1) - 1
c
c  form a leaf block, that is, one with no neighbors
c  in its next higher level.
c
100   level = nodlvl(leaf)
      nodlvl(leaf) = 0
      perm(num+1) = leaf
      blksze = 1
      call fnspan ( xadj, adjncy, nodlvl, blksze, perm(num+1),
     &  level, nadjs, adjs, nuleaf )
      if ( nuleaf .le. 0 ) go to 300
      jp = num

      do j = 1, blksze
        jp = jp + 1
        node = perm(jp)
        nodlvl(node) = level
      end do

      leaf = nuleaf
      go to 100
c
c  a new block has been found ...
c
300   nblks = nblks + 1
      xblk(nblks) = num + 1
      num = num + blksze
c
c  find the next possible block by using the adjacent
c  set in the lower level and the top node subset (if
c  appropriate) in the stack.
c
      level = level - 1
      if ( level .le. 0 ) go to 500
      call copysi ( nadjs, adjs, perm(num+1) )
      blksze = nadjs
      if ( level .ne. toplvl ) go to 400
c
c  the level of the node subset at the top of the
c  stack is the same as that of the adjacent set.
c  pop the node subset from the stack.
c
      npop = stack(topstk-1)
      topstk = topstk - npop - 2
      ip = num + blksze + 1
      call copysi ( npop, stack(topstk+1), perm(ip) )
      blksze = blksze + npop
      toplvl = stack(topstk)

400   continue

      call fnspan ( xadj, adjncy, nodlvl, blksze,
     &  perm(num+1), level, nadjs, adjs, nuleaf)

      if ( nuleaf .le. 0 ) go to 300
c
c  push the current node set into the stack.
c
      call copysi ( blksze, perm(num+1), stack(topstk+1) )
      topstk = topstk + blksze + 2
      stack(topstk-1) = blksze
      stack(topstk) = level
      toplvl = level
      leaf = nuleaf
      go to 100
c
c  before exit 
c
500   xblk(nblks+1) = num + 1
      return
      end
      subroutine setadj ( adj, irow, jcol, maxadj, nadj, neqn, xadj )

c*********************************************************************72
c
cc SETADJ sets up adjacency information in XADJ and ADJNCY.
c
c  the first call to setadj for a given problem should be with
c  irow=jcol=0, which is a flag telling setadj to clear out
c  the old information, and set xadj(i)=1 for i=1,neqn+1.
c
c  thereafter, just call setadj with values of irow and jcol
c  for each entry a(irow,jcol) of the matrix that is nonzero.
c  you only have to give the subdiagonal or superdiagonal
c  entries (irow.lt.jcol) or (jcol.lt.irow).  setadj automatically
c  'reflects' the information for the other half of the matrix.
c  repeated calls with the same values of irow and jcol do not
c  hurt.  no extra storage will be allocated.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    01 January 2009
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input/output, integer ADJ(MAXADJ), stores adjacency information.
c
c  irow   - the row index of a nonzero entry.
c           irow should be between 1 and neqn, except for the
c           initialization option, when irow=0 should be used.
c
c  jcol   - the column index of a nonzero entry.
c           jcol should be between 1 and neqn, except for the
c           initialization option, when jcol=0 should be used.
c
c  maxadj - the dimension of adjncy, the maximum amount of
c           storage you have allocated.
c
c    input/output, integer nadj, the number of adjaceny entries in adjncy.
c
c  neqn   - the number of equations represented by the matrix.
c           (the number of rows and columns).
c
c  xadj   - an integer vector, indexed by a row or column of the
c           matrix, and itself used to index adjncy.
c
      integer maxadj
      integer neqn

      integer adj(maxadj)
      integer i
      integer icall
      integer irow
      integer j
      integer jcol
      integer k
      integer kback
      integer khi
      integer klo
      integer nadj
      integer xadj(neqn+1)

      save icall

      data icall /0/

      if(irow.eq.0.or.jcol.eq.0.or.icall.eq.0)then
 
        icall = 1
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'SETADJ:'
        write ( *, '(a)' ) '  Initialize XADJ adjacency pointer array.'

        nadj = 0
        do i = 1, neqn + 1
          xadj(i) = 1
        end do
 
        return
 
      end if
 
      if ( irow .eq. jcol ) then
        return
      end if

      if(xadj(neqn+1).gt.maxadj+1)then
        write(*,*)' '
        write(*,*)'setadj - fatal error!'
        write(*,*)'  all available storage in adjncy has been used.'
        write(*,*)'  no more information can be stored!'
        stop
      endif
 
      if(irow.gt.neqn.or.irow.le.0.or.jcol.gt.neqn.or.jcol.le.0)then
        write(*,*)' '
        write(*,*)'setadj - fatal error!'
        write(*,*)'  illegal indices irow=',irow,' jcol=',jcol
        stop
      endif
 
      i=irow
      j=jcol

   20 continue

      klo=xadj(i)
      khi=xadj(i+1)-1
 
      do k=klo,khi

        if(adj(k).eq.j)then

          if(i.eq.irow)then
            i=jcol
            j=irow
            go to 20
          endif

          return

        endif

      enddo
 
      do k=xadj(i+1),xadj(neqn+1)
        kback=xadj(neqn+1)+xadj(i+1)-k
        adj(kback+1)=adj(kback)
      enddo
 
      adj(xadj(i+1))=j
 
      do k=i+1,neqn+1
        xadj(k)=xadj(k)+1
      end do
 
      nadj = xadj(neqn+1) - 1

      if(i.eq.irow)then
        i=jcol
        j=irow
        go to 20
      end if
 
      return
      end
      subroutine shomat ( adjncy, iband, ienv, invprm, maxadj, neqn,
     &  perm, xadj )

c*********************************************************************72
c
cc SHOMAT can display a symbolic picture of the matrix defined by
c  the adjacency information in xadj and adjncy, with a possible
c  permutation through perm and invprm.  shomat also computes
c  the bandwidth and the size of the envelope.
c
c  if no permutation has been done, you must set invprm(i)=perm(i)=i
c  before calling shomat.  otherwise, you must call invrse to
c  get the inverse permutation invprm before calling showmat.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    01 January 2009
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c  adjncy - an integer vector containing adjacency information.
c
c  iband  - compute by shomat, the bandwidth of the matrix.
c
c  ienv   - computed by shomat, the number of cells in the envelope.
c           you could think of this number as the sum of the
c           bandwidths of each row.
c
c  invprm - the inverse permutation vector.
c
c  maxadj - the dimension of adjncy.
c
c  neqn   - the number of equations.
c
c  xadj   - integer vector of dimension neqn+1, used to index
c           the vector adjncy.
c
      integer maxadj
      integer neqn

      integer adjncy(maxadj)
      character*1 band(100)
      integer i
      integer iadd
      integer iband
      integer ienv
      integer invprm(neqn)
      integer itemp
      integer j
      integer jhi
      integer jlo
      integer k
      integer perm(neqn)
      integer xadj(neqn+1)

      iband=0
      ienv=0
 
      if(neqn.gt.100)then
        write(*,*)' '
        write(*,*)'shomat - fatal error!'
        write(*,*)'  neqn is too large!'
        write(*,*)'  maximum legal value is 100.'
        write(*,*)'  your input value was ',neqn
        stop
      end if
 
      write(*,*)' '
      write(*,*)'shomat - display nonzero structure of matrix.'
      write(*,*)' '

        do i = 1, neqn

          do k = 1, neqn
            band(k) = ' '
          end do

          band(i) = 'x'

          iadd = 0
          jlo = xadj(perm(i))
          jhi = xadj(perm(i)+1)-1

          do j = jlo, jhi
            itemp = invprm(adjncy(j))
            if(i-itemp.gt.iband) iband=i-itemp
            band(itemp) = 'x'
            if(i-itemp.gt.iadd) iadd=i-itemp
          end do

          write(*,'(1x,i6,1x,100a1)')i,(band(j),j=1,neqn)
          ienv = ienv + iadd

        end do

      write(*,*)' '
      write(*,*)'lower bandwidth=',iband
      write(*,*)'lower envelope contains ',ienv,' entries.'

      return
      end
      subroutine smbfct ( neqns, xadj, adjncy, perm, invp,
     &  xlnz, maxlnz, xnzsub, nzsub, maxsub,
     &  rchlnk, mrglnk, marker, flag )

c*********************************************************************72
c
cc SMBFCT performs symbolic factorization on a permuted linear system.
c
c  Discussion:
c
c    The routine also sets up the compressed data structure for the system.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters -
c
c    Input, integer NEQNS, the number of equations.
c
c  (xadj, adjncy) - the adjacency structure.
c
c  (perm, invp) - the permutation vector and its inverse.
c
c  updated parameters 
c  maxsub - size of the subscript array nzsub. on return,
c  it contains the number of subscripts used
c
c  output parameters 
c  xlnz - index into the nonzero storage vector lnz.
c  (xnzsub, nzsub) - the compressed subscript vectors.
c  maxlnz - the number of nonzeros found.
c  flag - error flag. positive value indicates that.
c  nzsub array is too small.
c
c  working parameters 
c  mrglnk - a vector of size neqns. at the kth step,
c  mrglnk(k), mrglnk(mrglnk(k)) , .........
c  is a list containing all those columns l($,j)
c  with j less than k, such that its first off-
c  diagonal nonzero is l(k,j). thus, the
c  nonzero structure of column l(*,k) can be found
c  by merging that of such columns l(t,j) with
c  the structure of a(*,k).
c  k - a vector of size neqns. it is used
c  the structure of each column l(-,k). at the
c  end of the kth step,
c  rchlnk(k), rchlnk(rchlnk(k)), ........
c  is the list of positions of nonzeros in column k
c  of the factoi2 t
c
c  to accumulate
c
c  marker - an integer vector of length neqns. it is used
c  to test if mass symbolic elimination can be
c  performed. that is, it is used to check whether
c  the structure of the current column k being
c  processed is completely determined by the single
c  column mrglnk(k).
c
      integer adjncy(*), invp(*), mrglnk(*), nzsub(*)
      integer perm(*), rchlnk(*), marker(*)
      integer xadj(*), xlnz(*), xnzsub(*)
      integer flag, i, inz, j, jstop, jstrt, k, knz
      integer kxsub, mrgk, lmax, m, maxlnz, maxsub
      integer nabor, neqns, node, np1, nzbeg, nzend
      integer rchm, mrkflg

      flag = 0
c
c  initialization ...
c
      nzbeg = 1
      nzend = 0
      xlnz(1) = 1

      do k = 1, neqns
        mrglnk(k) = 0
        marker(k) = 0
      end do
c
c  for each column knz counts the number
c  of nonzeros in column k accumulated in rchlnk.
c
      np1 = neqns + 1

      do 1500 k = 1, neqns

        knz = 0
        mrgk = mrglnk(k)
        mrkflg = 0
        marker(k) = k
        if (mrgk .ne. 0 ) marker(k) = marker(mrgk)
        xnzsub(k) = nzend
        node = perm(k)
        jstrt = xadj(node)
        jstop = xadj(node+1) - 1
        if (jstrt.gt.jstop) go to 1500
c
c  use rchlnk to link through the structure of
c  a(*,k) below diagonal
c
        rchlnk(k) = np1

        do 300 j = jstrt, jstop
          nabor = adjncy(j)
          nabor = invp(nabor)
          if ( nabor .le. k ) go to 300
          rchm = k
200       continue
          m = rchm
          rchm = rchlnk(m)
          if ( rchm .le. nabor ) go to 200
          knz = knz+1
          rchlnk(m) = nabor
          rchlnk(nabor) = rchm
          if ( marker(nabor) .ne. marker(k) ) mrkflg = 1
300     continue
c
c  test for mass symbolic elimination ...
c
        lmax = 0
        if ( mrkflg .ne. 0 .or. mrgk .eq. 0 ) go to 350
        if ( mrglnk(mrgk) .ne. 0 ) go to 350
        xnzsub(k) = xnzsub(mrgk) + 1
        knz = xlnz(mrgk+1) - (xlnz(mrgk) + 1)
        go to 1400
c
c  link through each column i that affects l(*,k
c
350     i = k
400     i = mrglnk(i)
        if (i.eq.0) go to 800
        inz = xlnz(i+1) - (xlnz(i)+1)
        jstrt = xnzsub(i) + 1
        jstop = xnzsub(i) + inz
        if (inz.le.lmax) go to 500
        lmax = inz
        xnzsub(k) = jstrt
c
c  merge structure of l(*,i) in nzsub into rci
c
500     rchm = k

        do 700 j = jstrt, jstop
          nabor = nzsub(j)
600       m = rchm
          rchm = rchlnk(m)
          if (rchm.lt.nabor) go to 600
          if (rchm.eq.nabor) go to 700
          knz = knz+1
          rchlnk(m) = nabor
          rchlnk(nabor) = rchm
          rchm = nabor
700     continue

        go to 400
c
c  check if subscripts duplicate those of another column.
c
800     if (knz.eq.lmax) go to 1400
c
c  or if tail of k-1st column matches head of kth.
c
        if (nzbeg.gt.nzend) go to 1200
        i = rchlnk(k)
        
        do 900 jstrt=nzbeg,nzend
          if (nzsub(jstrt)-i) 900, 1000, 1200
900     continue

        go to 1200
1000    xnzsub(k) = jstrt

        do j = jstrt, nzend
          if (nzsub(j).ne.i) go to 1200
          i = rchlnk(i)
          if (i.gt.neqns) go to 1400
        end do

        nzend = jstrt - 1
c
c  copy the structure of l($,k) from rchlnk
c  to the data structure (xnzsub. nzsiir~
c
1200    nzbeg = nzend + 1
        nzend = nzend + knz
        if (nzend.gt.maxsub) go to 1600
        i = k
        do j=nzbeg,nzend
          i = rchlnk(i)
          nzsub(j) = i
          marker(i) = k
        end do

        xnzsub(k) = nzbeg
        marker(k) = k
c
c  update the vector mrglnk. note column l(e,k) just found
c  is required to determine column l(*,j), where
c  l(j,k) is the first nonzero in l(e,k) below diagonal.
c
1400    if (knz.le.1) go to 1500
        kxsub = xnzsub(k)
        i = nzsub(kxsub)
        mrglnk(k) = mrglnk(i)
        mrglnk(i) = k
1500    xlnz(k+1) = xlnz(k) + knz

      maxlnz = xlnz(neqns) - 1
      maxsub = xnzsub(neqns)
      xnzsub(neqns+1) = xnzsub(neqns)
      flag = 0
      return
c
c  error - insufficient storage for nonzero subscripts.
c
1600  flag=1
      return
      end
      subroutine sorts1 ( na, array )

c*********************************************************************72
c
cc SORTS1 uses linear insertion to sort integers into increasing order.
c
c  Modified:
c
c    01 January 2009
c
c  Parameters:
c
c    Input, integer NA, the array size.
c
c    Input/output, integer ARRAY(NA).  On output, the entries have been
c    sorted in increasing order.
c
      integer na

      integer array(na)
      integer k
      integer l
      integer node

      do k = 2, na
        node = array(k)
        l = k - 1
100     if(l.lt. 1) go to 200
        if ( array(l) .le. node) go to 200
        array(l+1) = array(l)
        l = l - 1
        go to 100
200	    array(l+1) = node
      end do

      return
      end
      subroutine subrcm ( xadj, adjncy, mask, nsubg, subg, perm, xls )

c*********************************************************************72
c
cc SUBRCM finds the rcm ordering for a given subgraph (possibly disconnected).
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c  (xadj, adjncy) - adjacency structure pair for the graph.
c
c  (nsubg, subg) - the given subgraph. nsubg is the
c  the size of the subgraph, and subg contains
c  the nodes in it.
c
c  output parameter 
c
c  perm - the permutation vector. it is also used
c  temporarily to store a level structure.
c
c  working parameters 
c
c  mask - mask vector with all zeros. it is used to
c  specify nodes in the subgraph.
c
c  xls - index to a level structure. note that the level
c  structure is stored in part of perm
c
      integer adjncy(*)
      integer ccsize
      integer i
      integer mask(*)
      integer nlvl
      integer node
      integer nsubg
      integer num
      integer perm(*)
      integer subg(*)
      integer xls(*)
      integer xadj(*)

      do i = 1, nsubg
        node = subg(i)
        mask(node) = 1
      end do

      num = 0

      do i = 1, nsubg

        node = subg(i)
c
c  For each connected component in the subgraph,
c  call FNROOT and RCM for the ordering.
c
        if ( 0 .lt. mask(node) ) then

          call fnroot ( node, xadj, adjncy, mask,
     &      nlvl, xls, perm(num+1) )

          call rcm ( node, xadj, adjncy, mask,
     &      perm(num+1), ccsize, xls )

          num = num + ccsize

          if ( nsubg .le. num ) then
            return
          end if

        end if

      end do

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c  Discussion:
c
c    This FORTRAN77 version is made available for cases where the
c    FORTRAN90 version cannot be used.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    12 January 2007
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    None
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, month(m), y, h, ':', n, ':', s, '.', mm, ampm

      return
      end
      subroutine tsfct ( nblks, xblk, father, diag, xenv, env,
     &  xnonz, nonz, nzsubs, temp, first, flag )

c*********************************************************************72
c
cc TSFCT performs the symmetric factorization of a tree-partitioned system.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c  (nblks, xblk, father) - the tree partitioning.
c
c  xenv - the envelope index vector.
c  (xnonz, nonz, nzsubs) - the off-diagonal nonzeros in
c  the original matrix.
c
c  updated parameters 
c  (diag, env) - storage arrays for the envelope of
c  the diagonal blocks of the matrix. on output,
c  contains the diagonal blocks of the factor.
c  flag - the error flag. it is set to 1 if a zero or
c  negative square root is detected during the
c  factorization.
c
c  working parameter 
c  temp - temporary array required to implement the
c  asymmetric version of the factorization.
c  first - temporary vector used to facilitate the
c  indexing to the vector nonz (or nzsubs)
c  for non-null subcolumns in off-diagonal
c  blocks.
c
      double precision diag(*)
      double precision env(*), nonz(*), temp(*), s
      integer father(*), nzsubs(*), xblk(*)
      integer first(*), xenv(*), xnonz(*)
      integer blksze, col, col1, colbeg, colend
      integer colsze, fnz, fnz1, i, flag, istrt, istop
      integer isub, j, jstop, jstrt, k, kenv, kenv0, kfathr
      integer nblks, neqns, row, rowbeg, rowend

      flag = 0
c
c  initialization.
c
      neqns = xblk(nblks+1) - 1

      do i = 1,neqns
        temp(i) = 0.0D+00
      end do

      do i = 1, neqns
        first(i) = xnonz(i)
      end do
c
c  loop through the blocks ...
c
      do 1600 k = 1, nblks
        rowbeg = xblk(k)
        rowend = xblk(k+1) - 1
        blksze = rowend - rowbeg + 1
        call esfct ( blksze, xenv(rowbeg), env,
     &    diag(rowbeg), flag )
        if ( flag .gt. 0 ) return
c
c  perform modification of the father diagonal block
c  a(father(k),father(k)) from the off-diagonal block
c  a(k,father(k)).
c
        kfathr = father(k)
        if ( kfathr .le. 0 ) go to 1600
        colbeg = xblk(kfathr)
        colend = xblk(kfathr+1) - 1
c
c  find the first and last non-null column in
c  the off-diagonal block. reset colbeg,colend.
c
        do col = colbeg, colend
          jstrt = first(col)
          jstop = xnonz(col+1) - 1
          if ( jstop .ge. jstrt .and.
     &      nzsubs(jstrt) .le. rowend ) go to 300
        end do

300     colbeg = col
        col = colend

        do col1 = colbeg, colend
          jstrt = first(col)
          jstop = xnonz(col+1) - 1
          if ( jstop .ge. jstrt .and.
     &      nzsubs(jstrt) .le. rowend ) go to 500
          col = col - 1
        end do

        colend = col
500     do 1300 col = colbeg, colend
          jstrt = first(col)
          jstop = xnonz(col+1) - 1
c
c  test for null subcolumn. fnz stores the
c  first nonzero subscript in the block column.
c
          if ( jstop .lt. jstrt ) go to 1300
          fnz = nzsubs(jstrt)
          if ( fnz .gt. rowend ) go to 1300
c
c  unpack a column in the off-diagonal block
c  and perform upper and lower solves on the
c  unpacked column.
c
          do j = jstrt, jstop
            row = nzsubs(j)
            if ( row .gt. rowend ) go to 700
            temp(row) = nonz(j)
          end do

700       colsze = rowend - fnz + 1
          call elslv ( colsze, xenv(fnz), env,
     &      diag(fnz), temp(fnz) )
          call euslv ( colsze, xenv(fnz), env,
     &      diag(fnz), temp(fnz) )
c
c  do the modification by looping through
c  the columns and forming inner products.
c
          kenv0 = xenv(col+1) - col
          do 1100 col1 = colbeg, colend
            istrt = first(col1)
            istop = xnonz(col1+1) - 1
c
c  check to see if subcolumn is null.
c
            fnz1 = nzsubs(istrt)
            if ( istop .lt. istrt .or.
     &        fnz1 .gt. rowend ) go to 1100
c
c  check if inner product should be done
c
            if ( fnz1 .lt. fnz ) go to 1100
            if ( fnz1 .eq. fnz .and.
     &        col1 .lt. col ) go to 1100

            s = 0.0D+00
            do i = istrt, istop
              isub = nzsubs(i)
              if ( isub .gt. rowend ) go to 900
              s = s + temp(isub) * nonz(i)
            end do
c
c  modify the env or the diag entry.
c
900         if ( col1 .eq. col ) go to 1000
            kenv = kenv0 + col1
            if ( col1 .gt. col )then
              kenv = xenv(col1+1) - col1 + col
            endif
            env(kenv) = env(kenv) - s
            go to 1100
1000        diag(col1) = diag(col1) - s
1100      continue
c
c  reset part of the temp vector to zero.
c
          do row = fnz, rowend
            temp(row) = 0.0D+00
          end do

1300    continue
c
c  update the first vector for columns in
c  father(k) block, so that it will index to
c  the beginning of the next off-diagonal
c  block to be considered.
c
        do 1500 col = colbeg, colend

          jstrt = first(col)
          jstop = xnonz(col+1) - 1
          if ( jstop .lt. jstrt ) go to 1500

          do 1400 j = jstrt, jstop
            row = nzsubs(j)
            if ( row .le. rowend ) go to 1400
            first(col) = j
            go to 1500
1400      continue

        first(col) = jstop + 1
1500    continue

1600  continue

      return
      end
      subroutine tsslv ( nblks, xblk, diag, xenv, env,
     &  xnonz, nonz, nzsubs, rhs, temp )

c*********************************************************************72
c
cc TSSLV solves a tree-partitioned factored system by implicit back substitution.
c
c  Modified:
c
c    01 January 2009
c
c  input parameters 
c
c  (nblks, xblk) - the partitioning.
c  (xenv, env) - envelope of the diagonal blocks.
c  (xnonz, nonz, nzsubs) - data structure for the off
c  block diagonal nonzeros.
c
c  updated parameters 
c  rhs - on input it contains the right hand vector
c  on output, the solution vector.
c
c  working vector 
c  temp - temporary vector used in back substitution.
c
      double precision diag(*)
      double precision env(*)
      double precision nonz(*)
      double precision rhs(*)
      double precision s
      double precision temp(*)

      integer nzsubs(*), xblk(*)
      integer xenv(*), xnonz(*), col, col1, col2, i, j
      integer jstop, jstrt, last, nblks, ncol, nrow, row
      integer row1, row2
c
c  forward substitution 
c
      do 400 i = 1, nblks

        row1 = xblk(i)
        row2 = xblk(i+1) - 1
        last = xnonz(row2+1)
        if ( i .eq. 1 .or. last .eq. xnonz(row1) ) go to 300
c
c  modify rhs vector by the product of the off-
c  diagonal block with the corresponding part of rhs.
c
        do 200 row = row1, row2

          jstrt = xnonz(row)
          if ( jstrt .eq. last ) go to 300
          jstop = xnonz(row+1) - 1
          if ( jstop .lt. jstrt ) go to 200
          s = 0.0D+00

          do j = jstrt, jstop
            col = nzsubs(j)
            s = s + rhs(col)*nonz(j)
          end do

          rhs(row) = rhs(row) - s

200     continue

300     nrow = row2 - row1 + 1
        call elslv ( nrow, xenv(row1), env, diag(row1),
     &    rhs(row1) )
        call euslv ( nrow, xenv(row1), env, diag(row1),
     &    rhs(row1) )

400   continue
c
c  backward solution 
c
      if ( nblks .eq. 1 ) return
      last = xblk(nblks) - 1
      do i = 1, last
        temp(i) = 0.0D+00
      end do

      i = nblks
      col1 = xblk(i)
      col2 = xblk(i+1)-1
600   if (i.eq.1)return
      last = xnonz(col2+1)
      if ( last .eq. xnonz(col1) ) go to 900
c
c  multiply off-diagonal block by the corresponding
c  part of the solution vector and store in temp.
c
      do 800 col = col1, col2

        s = rhs(col) 
        if (s.eq.0.0D+00 ) go to 800
        jstrt = xnonz(col)
        if ( jstrt .eq. last) go to 900
        jstop = xnonz(col+1) - 1
        if ( jstop .lt. jstrt ) go to 800

        do j = jstrt, jstop
          row = nzsubs(j)
          temp(row) = temp(row) + s*nonz(j)
        end do

800   continue

900   i = i - 1
      col1 = xblk(i)
      col2 = xblk(i+1) - 1
      ncol = col2 - col1 + 1
      call elslv ( ncol, xenv(col1), env,
     &  diag(col1), temp(col1) )

      call euslv ( ncol, xenv(col1), env, diag(col1),
     &  temp(col1) )

      do j = col1, col2
        rhs(j) = rhs(j) - temp(j)
      end do

      go to 600
      end