!  File:      genqmd.f
!  Date:      12/31/01  (downloaded fall 2001)
!  Author:    C. Paciorek - taken off the web
!  Purpose:   function for reordering rows and columns of a pos def matrix to achieve increased
!             sparseness of Cholesky factor
!  Compile:   g77 -c -ffree-form
!  Usage:     
!  Arguments: 
!  Requires:  
!  Reference: sparspak.f90 from the web -- George & Liu 1981, see below

     subroutine genqmd( n, xadj, iadj, perm, invprm, deg,marker,&
    rchset, nbrhd,qsize, qlink, nofsub )
!
!*******************************************************************************
!
!! GENQMD implements the quotient minimum degree algorithm.  
!  
!  namely, it reorders rows and columns to achieve increased sparseness
!  of the Cholesky factor - CJP 12/31/01
!
!  Code taken from sparspak.f90 on the web (CJP 9/4/01)
!
!  Discussion:
!
!    The routine uses the implicit representation of the elimination 
!    graphs by quotient graphs, and the notion of indistinguishable nodes.
!
!    The adjacency vector IADJ will be destroyed.
!
!  Reference:
!
!    Alan George and J W Liu,
!    Computer Solution of Large Sparse Positive Definite Systems,
!    Prentice Hall, 1981.
!
!  Parameters:
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer XADJ(N+1).  Information about row I is stored
!    in entries XADJ(I) through XADJ(I+1)-1 of IADJ.
!
!    Input, integer IADJ(*), the adjacency structure. 
!    For each row, it contains the column indices of the nonzero entries.
!
!  output parameters-
!  perm-the minimum degree ordering.
!  invprm-the inverse of perm.
!
!  working parameters-
!  marker-a marker vector, where marker(i) is
!  negative means node i has been merged with
!  another node and thus can be ignored.
!  rchset-vector used for the reachable set.
!  nbrhd-vector used for the neighborhood set.
!  qsize-vector used to store the size of
!  indistinguishable supernodes.
!  qlink-vector to store indistinguishable nodes,
!  i,qlink(i),qlink(qlink(i)) are the
!  members of the supernode represented by i.
                                !
      integer n
  !
      integer iadj(*)
      integer deg(n)
      integer inode
      integer invprm(n)
      integer ip
      integer irch
      integer j
      integer marker(n)
      integer mindeg
      integer nbrhd(*)
      integer ndeg
      integer nhdsze
      integer node
      integer nofsub
      integer np
      integer num
      integer nump1
      integer nxnode
      integer perm(n)
      integer qlink(n)
      integer qsize(n)
      integer rchset(*)
      integer rchsze
      integer search
      integer thresh
      integer xadj(n+1)
!
!  Initialize degree vector and other working variables.
!
!  DEG(I) negative means node I has been numbered.
!
      mindeg = n
      nofsub = 0
      
      do node = 1, n
         perm(node) = node
         invprm(node) = node
         marker(node) = 0
         qsize(node) = 1
         qlink(node) = 0
         ndeg = xadj(node+1)-xadj(node)
         deg(node) = ndeg
         mindeg = min ( mindeg, ndeg )
      end do
      
      num = 0
                                !
                                !  Perform threshold search to get a node of minimum degree.
                                !  Variable search points to where search should start.
                                !
 20   continue
      
      search = 1
      thresh = mindeg
      mindeg = n
      
 30   continue
      
      nump1 = num+1
      search = max ( search, nump1 )
      
      do j = search, n
         
         node = perm(j)
         
         if ( marker(node) >= 0 ) then
            ndeg = deg(node)
            if ( ndeg <= thresh ) then
               goto 50
            end if
            mindeg = min ( mindeg, ndeg )
         end if
            
      end do
      
      goto 20
                                !
                                !  The node has minimum degree. 
                                !  Find its reachable sets by calling qmdrch.
                                !
 50   continue
      
      search = j
      nofsub = nofsub+deg(node)
      marker(node) = 1
      
      call qmdrch ( node, xadj, iadj, deg, marker, rchsze, rchset, &
      nhdsze, nbrhd, n )
                                !
                                !  Eliminate all nodes indistinguishable from node.
                                !  They are given by node, qlink(node),
                                !
      nxnode = node
      
 60   continue
      
      num = num + 1
      np = invprm(nxnode)
      ip = perm(num)
      perm(np) = ip
      invprm(ip) = np
      perm(num) = nxnode
      invprm(nxnode) = num
      deg(nxnode) = -1
      nxnode = qlink(nxnode)
      
      if ( nxnode > 0 ) then
         goto 60
      end if
      
      if ( rchsze <= 0 ) then
         goto 80
      end if
                                !
                                !  Update the degrees of the nodes in the reachable
                                !  set and identify indistinguishable nodes.
                                !
      call qmdupd ( xadj, iadj, rchsze, rchset, deg, qsize, qlink, &
      marker, rchset(rchsze+1), nbrhd(nhdsze+1), n )
                                !
                                !  Reset marker value of nodes in reach set.
                                !  Update threshold value for cyclic search.
                                !  Call qmdqt to form new quotient graph.
                                !
      marker(node) = 0
      
      do irch = 1, rchsze
         
         inode = rchset(irch)
         
         if ( marker(inode) >= 0 ) then
         
         marker(inode) = 0
         ndeg = deg(inode)
         mindeg = min ( mindeg, ndeg )
         
         if ( ndeg <= thresh ) then
            mindeg = thresh
            thresh = ndeg
            search = invprm(inode)
         end if
         
      end if
         
      end do
      
      if ( nhdsze > 0 ) then
        call qmdqt ( node, xadj, iadj, marker, rchsze,rchset,nbrhd, n)
      end if
      
 80   continue
      
      if ( num < n ) then
         goto 30
      end if
      
      return
      end subroutine genqmd
      
      subroutine qmdqt ( root, xadj, iadj, marker, rchsze, rchset,& 
       nbrhd,n )
                                !
                                !*******************************************************************************
                                !
                                !! QMDQT performs the quotient graph transformation after a node has
                                !  been eliminated.
                                !
                                !
                                !  Reference:
                                !
                                !    Alan George and J W Liu,
                                !    Computer Solution of Large Sparse Positive Definite Systems,
                                !    Prentice Hall, 1981.
                                !
                                !  Parameters:
                                !
                                !    Input, integer ROOT, the node just eliminated.  It becomes the
                                !    representative of the new supernode.
                                !
                                !    Input, integer XADJ(N+1).  Information about row I is stored
                                !    in entries XADJ(I) through XADJ(I+1)-1 of IADJ.
                                !
                                !    Input, integer IADJ(*), the adjacency structure. 
                                !    For each row, it contains the column indices of the nonzero entries.
                                !
                                !  input parameters-
                                !
                                !  (rchsze,rchset)-the reachable set of root in the
                                !  old quotient graph.
                                !  nbrhd-the neighborhood set which will be merged
                                !  with root to form the new supernode.
                                !  marker-the marker vector.
                                !
                                !  updated parameter-
                                !  iadj-becomes the adjacency structure of the quotient graph.
                                !
      integer n
                                !
      integer iadj(*)
      integer inhd
      integer irch
      integer j
      integer jstrt
      integer jstop
      integer link
      integer marker(n)
      integer nabor
      integer nbrhd(*)
      integer node
      integer rchset(*)
      integer rchsze
      integer root
      integer xadj(n+1)
                                !
      irch = 0
      inhd = 0
      node = root
      
 10   continue
      
      jstrt = xadj(node)
      jstop = xadj(node+1)-2
                                !
                                !  Place reach nodes into the adjacent list of node.
                                !
      do j = jstrt, jstop
         irch = irch+1
         iadj(j) = rchset(irch)
         if ( irch >= rchsze ) then
            goto 40
         end if
      end do
                                !
                                !  Link to other space provided by the nbrhd set.
                                !
      link = iadj(jstop+1)
      node = -link
      
      if ( link >= 0 ) then
         inhd = inhd + 1
         node = nbrhd(inhd)
         iadj(jstop+1) = -node
      end if
      
      goto 10
                                !
                                !  All reachable nodes have been saved.  End the adj list.
                                !  Add root to the nbr list of each node in the reach set.
                                !
 40   continue
      
      iadj(j+1) = 0
      
      do irch = 1, rchsze
         
         node = rchset(irch)
         
         if ( marker(node) >= 0 ) then
         
         jstrt = xadj(node)
         jstop = xadj(node+1)-1
         
         do j = jstrt, jstop
            
            nabor = iadj(j)
            
            if ( marker(nabor) < 0 ) then
               iadj(j) = root
               goto 60
            end if
            
         end do
         
      end if
         
 60      continue
         
      end do
      
      return
      end subroutine qmdqt
      subroutine qmdrch ( root, xadj, iadj, deg, marker, rchsze,&
       rchset, nhdsze, nbrhd, n )
                                !
                                !*******************************************************************************
                                !
                                !! QMDRCH determines the reachable set of a node through a given subset.  
                                !
                                !
                                !  Discussion:
                                !
                                !    The adjacency structure is assumed to be stored in a quotient graph format.
                                !
                                !  Reference:
                                !
                                !    Alan George and J W Liu,
                                !    Computer Solution of Large Sparse Positive Definite Systems,
                                !    Prentice Hall, 1981.
                                !
                                !  Parameters:
                                !
                                !    Input, integer ROOT, the given node, which is not in the subset.
                                !
                                !    Input, integer XADJ(N+1).  Information about row I is stored
                                !    in entries XADJ(I) through XADJ(I+1)-1 of IADJ.
                                !
                                !    Input, integer IADJ(*), the adjacency structure. 
                                !    For each row, it contains the column indices of the nonzero entries.
                                !
                                !  input parameters-
                                !
                                !  deg-the degree vector.  deg(i) lt 0 means the node
                                !  belongs to the given subset.
                                !
                                !  output parameters-
                                !  (rchsze,rchset)-the reachable set.
                                !  (nhdsze,nbrhd)-the neighborhood set.
                                !
                                !  updated parameters-
                                !  marker-the marker vector for reach and nbrhd sets.
                                !  gt 0 means the node is in reach set.
                                !  lt 0 means the node has been merged with
                                !  others in the quotient or it is in nbrhd set.
                                !
      integer n
                                !
      integer iadj(*)
      integer deg(n)
      integer i
      integer istop
      integer istrt
      integer j
      integer jstop
      integer jstrt
      integer marker(n)
      integer nabor
      integer nbrhd(*)
      integer nhdsze
      integer node
      integer rchset(*)
      integer rchsze
      integer root
      integer xadj(n+1)
                                !
                                !  Loop through the neighbors of ROOT in the quotient graph.
                                !
      nhdsze = 0
      rchsze = 0
      istrt = xadj(root)
      istop = xadj(root+1)-1
      
      do i = istrt, istop
         
         nabor = iadj(i)
         if ( nabor == 0 ) then
         return
      end if
         
      if ( marker(nabor) /= 0 ) then
         goto 50
      end if
     !        
     !        Include NABOR in the reachable set.
     !        
      if ( deg(nabor) >= 0 ) then
         
         rchsze = rchsze+1
         rchset(rchsze) = nabor
         marker(nabor) = 1
         goto 50
         
      end if
     !        
     !        NABOR has been eliminated.  Find nodes reachable from it.
     !        
      marker(nabor) = -1
      nhdsze = nhdsze+1
      nbrhd(nhdsze) = nabor
      
 20   continue
      
      jstrt = xadj(nabor)
      jstop = xadj(nabor+1)-1
      
      do j = jstrt, jstop
         
         node = iadj(j)
         nabor = -node
         
         if ( node < 0 ) then
            goto 20
         end if
         
         if ( node == 0 ) then
            goto 50
         end if
         
         if ( marker(node) == 0 ) then
            rchsze = rchsze+1
            rchset(rchsze) = node
            marker(node) = 1
         end if
         
      end do
         
 50      continue
         
      end do
      
      return
      end subroutine qmdrch
      subroutine qmdupd ( xadj, iadj, nlist, list, deg, qsize, qlink,&
       marker, rchset, nbrhd, n )
                                !
                                !*******************************************************************************
                                !
                                !! QMDUPD performs degree update for a set of nodes in the minimum
                                !  degree algorithm.
                                !
                                !
                                !  Reference:
                                !
                                !    Alan George and J W Liu,
                                !    Computer Solution of Large Sparse Positive Definite Systems,
                                !    Prentice Hall, 1981.
                                !
                                !  Parameters:
                                !
                                !    Input, integer XADJ(N+1).  Information about row I is stored
                                !    in entries XADJ(I) through XADJ(I+1)-1 of IADJ.
                                !
                                !    Input, integer IADJ(*), the adjacency structure. 
                                !    For each row, it contains the column indices of the nonzero entries.
                                !
                                !  input parameters-
                                !
                                !  (nlist,list)-the list of nodes whose degree has to
                                !  be updated.
                                !
                                !  updated parameters-
                                !  deg-the degree vector.
                                !  qsize-size of indistinguishable supernodes.
                                !  qlink-linked list for indistinguishable nodes.
                                !  marker-used to mark those nodes in reach/nbrhd sets.
                                !
                                !  working parameters-
                                !  rchset-the reachable set.
                                !  nbrhd- the neighborhood set.
                                !
      integer n
                                !
      integer iadj(*)
      integer deg(n)
      integer deg0
      integer deg1
      integer il
      integer inhd
      integer inode
      integer irch
      integer j
      integer jstop
      integer jstrt
      integer list(*)
      integer mark
      integer marker(n)
      integer nabor
      integer nbrhd(*)
      integer nhdsze
      integer nlist
      integer node
      integer qlink(n)
      integer qsize(n)
      integer rchset(*)
      integer rchsze
      integer xadj(n+1)
                                !
                                !  Find all eliminated supernodes that are adjacent to some nodes in the
                                !  given list.  Put them into (nhdsze,nbrhd).  DEG0 contains the number of
                                !  nodes in the list.
                                !
      if ( nlist <= 0 ) then
         return
      end if
      
      deg0 = 0
      nhdsze = 0
      
      do il = 1, nlist
         
         node = list(il)
         deg0 = deg0+qsize(node)
         jstrt = xadj(node)
         jstop = xadj(node+1)-1
         
         do j = jstrt, jstop
         
            nabor = iadj(j)
         
            if ( marker(nabor) == 0 .and. deg(nabor) < 0 ) then
               marker(nabor) = -1
               nhdsze = nhdsze+1
               nbrhd(nhdsze) = nabor
            end if
            
         end do
         
      end do
                                !
                                !  Merge indistinguishable nodes in the list by calling QMDMRG.
                                !
      if ( nhdsze > 0 ) then
         call qmdmrg ( xadj, iadj, deg, qsize, qlink, marker, deg0,& 
        nhdsze, nbrhd, rchset, nbrhd(nhdsze+1), n )
      end if
                                !
                                !  Find the new degrees of the nodes that have not been merged.
                                !
      do il = 1, nlist
         
         node = list(il)
         mark = marker(node)
         
         if ( mark == 0 .or. mark == 1 ) then
         
         marker(node) = 2
         
         call qmdrch ( node, xadj, iadj, deg, marker, rchsze, rchset,&
         nhdsze, nbrhd, n )
         
         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
         
      end if
         
      end do
      
      return
      end subroutine qmdupd
      
      subroutine qmdmrg ( xadj, iadj, deg, qsize, qlink, marker, deg0,&
       nhdsze, nbrhd, rchset, ovrlp, n )
                                !
                                !*******************************************************************************
                                !
                                !! QMDMRG merges indistinguishable nodes in the minimum degree ordering
                                !  algorithm, and computes the new degrees of these new supernodes.
                                !
                                !
                                !  Reference:
                                !
                                !    Alan George and J W Liu,
                                !    Computer Solution of Large Sparse Positive Definite Systems,
                                !    Prentice Hall, 1981.
                                !
                                !  Parameters:
                                !
                                !    Input, integer XADJ(N+1).  Information about row I is stored
                                !    in entries XADJ(I) through XADJ(I+1)-1 of IADJ.
                                !
                                !    Input, integer IADJ(*), the adjacency structure. 
                                !    For each row, it contains the column indices of the nonzero entries.
                                !
                                !  input parameters-
                                !
                                !  deg0-the number of nodes in the given set.
                                !  (nhdsze,nbrhd)-the set of eliminated supernodes
                                !  adjacent to some nodes in the set.
                                !
                                !  updated parameters-
                                !  deg-the degree vector.
                                !  qsize-size of indistinguishable nodes.
                                !  qlink-linked list for indistinguishable nodes.
                                !  marker-the given set is given by those nodes with
                                !  marker value set to 1.  those nodes with degree
                                !  updated will have marker value set to 2.
                                !
                                !  working parameters-
                                !  rchset-the reachable set.
                                !  ovrlp- temp vector to store the intersection of two
                                !  reachable sets.
                                !
      integer n
                                !
      integer iadj(*)
      integer deg(n)
      integer deg0
      integer deg1
      integer head
      integer inhd
      integer iov
      integer irch
      integer j
      integer jstrt
      integer jstop
      integer link
      integer lnode
      integer mark
      integer marker(n)
      integer mrgsze
      integer nabor
      integer nbrhd(*)
      integer nhdsze
      integer node
      integer novrlp
      integer ovrlp(*)
      integer qlink(n)
      integer qsize(n)
      integer rchset(*)
      integer rchsze
      integer root
      integer xadj(n+1)
                                !
      if ( nhdsze <= 0 ) then
         return
      end if
      
      do inhd = 1, nhdsze
         root = nbrhd(inhd)
         marker(root) = 0
      end do
                                !
                                !  Loop through each eliminated supernode in the set (nhdsze,nbrhd).
                                !
      do inhd = 1, nhdsze
         
         root = nbrhd(inhd)
         marker(root) = - 1
         rchsze = 0
         novrlp = 0
         deg1 = 0
         
 20      continue
         
         jstrt = xadj(root)
         jstop = xadj(root+1)-1
     !        
     !        Determine the reachable set and its intersection with the input
     !        reachable set.
     !        
         do j = jstrt, jstop
         
         nabor = iadj(j)
         root = -nabor
         
         if ( nabor < 0 ) then
            goto 20
         end if
         
         if ( nabor == 0 ) then
            goto 70
         end if
         
         mark = marker(nabor)
         
         if ( mark == 0 ) then
            
            rchsze = rchsze+1
            rchset(rchsze) = nabor
            deg1 = deg1+qsize(nabor)
            marker(nabor) = 1
            
         else if ( mark == 1 ) then
            
            novrlp = novrlp+1
            ovrlp(novrlp) = nabor
            marker(nabor) = 2
            
         end if
         
      end do
     !           
     !           From the overlapped set, determine the nodes that can be merged.
     !           
 70      continue
         
         head = 0
         mrgsze = 0
         
         do iov = 1, novrlp
         
         node = ovrlp(iov)
         jstrt = xadj(node)
         jstop = xadj(node+1)-1
         
         do j = jstrt, jstop
            
            nabor = iadj(j)
            
            if ( marker(nabor) == 0 ) then
               marker(node) = 1
               goto 110
            end if
            
         end do
                                !
                                !  NODE belongs to the new merged supernode.
                                !  Update the vectors QLINK and QSIZE.
                                !
         mrgsze = mrgsze + qsize(node)
         marker(node) = -1
         lnode = node
         
 90      continue
         
         link = qlink(lnode)
         
         if ( link > 0 ) then
            lnode = link
            goto 90
         end if
         
         qlink(lnode) = head
         head = node
         
 110     continue
         
      end do
         
      if ( head > 0 ) then
         qsize(head) = mrgsze
         deg(head) = deg0+deg1-1
         marker(head) = 2
      end if
                                !        
                                !        Reset marker values.
                                !        
      root = nbrhd(inhd)
      marker(root) = 0
      
      do irch = 1, rchsze
         node = rchset(irch)
         marker(node) = 0
      end do
      
      end do
      
      return
      end subroutine qmdmrg
      
      
      
      
      
      
      
      
      










