!  File:      smbfct.f
!  Date:      12/31/01  (downloaded fall 2001)
!  Author:    C. Paciorek - taken off the web
!  Purpose:   Cholesky symbolic factorization of pos def matrix stored in sparse format
!             namely, sets up sparse structure and col/row info of the Cholesky
!  Compile:   g77 -c -ffree-form
!  Usage:     
!  Arguments: 
!  Requires:  
!  Reference: sparspak.f90 from the web -- George & Liu 1981, see below

subroutine smbfct( n, xadj, iadj, perm, invprm, ixlnz, nofnz, &
      xnzsub, nzsub, maxsub, rchlnk, mrglnk ,marker)
!
!*******************************************************************************
!
!! SMB_FACTOR performs symbolic factorization on a permuted linear system.
!
!  taken from sparspak.f90 by CJP (9/4/01)
!
!
!  Discussion:
!
!    The routine also sets up the compressed data structure for the system.
!
!  Modified:
!
!    30 April 2000
!
!  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.
!
!    Input, integer PERM(N), the permutation vector.
!
!    Input, integer INVPRM(N), the inverse permutation vector.
!
!  updated parameters-
!  maxsub-size of the subscript array nzsub.  on return,
!  it contains the number of subscripts used
!
!  output parameters-
!  ixlnz-index into the nonzero storage vector lnz.
!  (xnzsub,nzsub)-the compressed subscript vectors.
!  nofnz-the number of nonzeros found.
!
!  working parameters-
!  mrglnk-a vector of size N.  at the kth step,
!  mrglnk(k),mrglnk(mrglnk(k)) ,
!  is a list containing all those columns l(*,j)
!  with j less than k,such that its first off-
!  diagonal nonzero is l(k,j).  thus,the
!  nonzero structure of column l(*,k) can be found
!  by merging that of such columns l(*,j) with
!  the structure of a(*,k).
!  rchlnk-a vector of size N.  it is used to accumulate
!  the structure of each column l(*,k).  at the
!end of the kth step,
!  rchlnk(k),rchlnk(rchlnk(k)),
!  is the list of positions of nonzeros in column k
!  of the factor l.
!
  integer n
!
  integer iadj(*)
  integer i
  integer inz
  integer invprm(n)
  integer ixlnz(n+1)
  integer j
  integer jstop
  integer jstrt
  integer k
  integer knz
  integer kxsub
  integer lmax
  integer m
  integer marker(n)
  integer maxsub
  integer mrgk
  integer mrglnk(n)
  integer mrkflg
  integer nabor
  integer node
  integer nofnz
  integer np1
  integer nzbeg
  integer nzend
  integer nzsub(*)
  integer perm(n)
  integer rchlnk(n)
  integer rchm
  integer xadj(n+1)
  integer xnzsub(*)
!
!  Initialization.
!
  nzbeg = 1
  nzend = 0
  ixlnz(1) = 1

  do k = 1, n
    mrglnk(k) = 0
    marker(k) = 0
  end do
!
!  MARKER is used to test if mass symbolic elimination can be performed.  
!  That is, it is used to check whether the structure of the current 
!  column K being processed is completely determined by the single
!  column MRGLNK(K).
!
!  marker = 0  ! this removed by CJP 9/25 (compilation error)
!  line with marker(k) = 0 added 9/25 based on code in George & Liu
!
!  For each column KNZ counts the number
!  of nonzeros in column K accumulated in rchlnk.
!
  np1 = n + 1

  do k = 1, n

    knz = 0
    mrgk = mrglnk(k)
    mrkflg = 0
    marker(k) = k
    if ( mrgk /= 0 ) then
      marker(k) = marker(mrgk)
    end if
    xnzsub(k) = nzend
!       write (*,*) ' k, xnzsub ', k, xnzsub(k)
   node = perm(k)
    jstrt = xadj(node)
    jstop = xadj(node+1)-1
    if ( jstrt > jstop ) then
      goto 160
    end if
!
!  Use RCHLNK to link through the structure of A(*,K) below diagonal.
!
    rchlnk(k) = np1

    do j = jstrt, jstop

      nabor = iadj(j)
      nabor = invprm(nabor)

      if ( nabor > k ) then

        rchm = k

   20       continue

        m = rchm
        rchm = rchlnk(m)
        if ( rchm <= nabor ) then
          goto 20
        end if

        knz = knz+1
        rchlnk(m) = nabor
        rchlnk(nabor) = rchm

        if ( marker(nabor) /= marker(k) ) then
          mrkflg = 1
        end if

      end if

    end do
!
!  Test for mass symbolic elimination
!
    lmax = 0

    if ( mrkflg /= 0 .or. mrgk == 0 ) then
      goto 40
    end if

    if ( mrglnk(mrgk) /= 0 ) then
      goto 40
    end if

    xnzsub(k) = xnzsub(mrgk) + 1
!      write (*,*) ' k, xnzsub ', k, xnzsub(k)
    knz = ixlnz(mrgk+1) - (ixlnz(mrgk)+1)
    goto 150
!
!  Link through each column I that affects l(*,k).
!
   40   continue

    i = k

   50   continue

    i = mrglnk(i)
    if ( i == 0 ) then
      goto 90
    end if

    inz = ixlnz(i+1)-(ixlnz(i)+1)
    jstrt = xnzsub(i)+1
    jstop = xnzsub(i)+inz

    if ( inz > lmax ) then
      lmax = inz
      xnzsub(k) = jstrt
 !     write (*,*) ' k, xnzsub ', k, xnzsub(k)
    end if
!
!  Merge structure of l(*,i) in nzsub into RCHLNK.
!
    rchm = k

    do j = jstrt, jstop

      nabor = nzsub(j)

   70     continue

      m = rchm
      rchm = rchlnk(m)
      if ( rchm < nabor ) then
        goto 70
      end if

      if ( rchm /= nabor ) then
        knz = knz+1
        rchlnk(m) = nabor
        rchlnk(nabor) = rchm
        rchm = nabor
      end if

    end do

    goto 50
!
!  Check if subscripts duplicate those of another column...
!
   90   continue

    if ( knz == lmax ) then
      goto 150
    end if
!
!  ...or if tail of column K-1 matches head of column K.
!
    if ( nzbeg > nzend ) then
      goto 130
    end if

    i = rchlnk(k)

    do jstrt = nzbeg, nzend

      if ( nzsub(jstrt) == i ) then
        goto 110
      end if

      if ( nzsub(jstrt) > i ) then
        goto 130
      end if

    end do

    goto 130

  110   continue

    xnzsub(k) = jstrt

    do j = jstrt, nzend
      if ( nzsub(j) /= i ) then
        goto 130
      end if
      i = rchlnk(i)
      if ( i > n ) then
        goto 150
      end if
    end do

    nzend = jstrt-1
!
!  Copy the structure of l(*,k) from rchlnk
!  to the data structure (xnzsub,nzsub).
!
  130   continue

    nzbeg = nzend+1
    nzend = nzend+knz

    if ( nzend > maxsub ) then
      write ( *, * ) ' '
      write ( *, * ) 'SMB_FACTOR - Fatal error!'
      write ( *, * ) '  Insufficient storage for nonzero entries.'
      write ( *, * ) '  MAXSUB  = ', maxsub
      write ( *, * ) '  Exceeded by NZEND = ', nzend
      stop
    end if

    i = k

    do j = nzbeg, nzend
      i = rchlnk(i)
      nzsub(j) = i
      marker(i) = k
    end do

    xnzsub(k) = nzbeg
!      write (*,*) ' k, xnzsub ', k, xnzsub(k)
    marker(k) = k
!
!  Update the vector mrglnk.  Note column L(*,K) just found
!  is required to determine column L(*,J), where
!  L(J,K) is the first nonzero in L(*,K) below diagonal.
!
  150   continue

    if ( knz > 1 ) then
      kxsub = xnzsub(k)
      i = nzsub(kxsub)
      mrglnk(k) = mrglnk(i)
      mrglnk(i) = k
    end if

  160   continue

    ixlnz(k+1) = ixlnz(k) + knz

  end do

  nofnz = ixlnz(n) - 1
  maxsub = xnzsub(n)
  xnzsub(n+1) = xnzsub(n)

  return
end
