subroutine preord(n, isym, ja, a, list, index, ierr)
c
c       the routines fillin, factor, and solve in this package are
c       prototype routines based on sparse elimination algorithms
c       and data structures described in:
c
c       General Sparse Elimination Requires No Permanent Integer Storage
c       by Randolph E. Bank and R. Kent Smith, vol. 8, pp.574-584, SIAM
c       Journal of Scientific and Statistical Computing, 1987.
c
c       a more complete discussion of these procedures can be found there.
c       to keep the code as simple as possible, we have assumed that the
c       identity ordering is used. if a reordering algorithm such as
c       minimum degree is to be used, one must either physically
c       reorder the matrix, or systematically modify references
c       to the arrays ja, a, x, and b to reflect the permutation matrix.
c       this is a shortcoming of the code, which was written for mainly
c       for illustrative purposes, and not of the underlying algorithms.
c
c       fillin and factor assume a particular order for the row indices
c       in ja. this routine is included to physically reorder ja and
c       a in the required way. it is not an essential part of the package,
c       and should be called only if necessary.
c
           integer ja(*), list(*), index(*)
           real a(*), ltemp, utemp
c
c       put row indices in ja in decreasing order
c
c       input: n, isym, ja, a 
c       output: ja, a, ierr
c       workspace: list(n), index(n)
c
c       remarks:
c
c           1. isym = 0 means nonsymmetric storage is used for a
c              isym = 1 means symmetric storage is used for a
c
c           2. ierr = 0 is a normal return
c              ierr = i means a row index greater than i was
c                       encountered in column i
c
c           3. if a is symmetric, the reordering of the lower
c              triangle in the do 40 loop is redundant
c
        ierr = 0
c
c       compute offset for a
c
        if (isym .eq. 1) then
           lshift = 0
        else
           lshift = ja(n+1) - ja(1)
        end if
c
c       the main loop
c  
        do 50 i = 1, n
           if (ja(i)+1 .ge. ja(i+1))  go to 50
c
c       sort seed indices in increasing order
c
           list(i) = i
           do 30 indexj = ja(i), ja(i+1)-1
              j = ja(indexj)
              if (j .ge. i)  go to 100
              index(j) = indexj
              last = i
              k = list(i)
   10         if (j .lt. k)  go to 20
              last = k
              k = list(k)
              go to 10
   20         list(last) = j
              list(j) = k
   30      continue
c
c       in place reordering of ja and a
c
           k = list(i)
           do 40 indexj = ja(i+1)-1, ja(i), -1
              j = ja(indexj)
              if (j .eq. k) go to 40
c
              indexk = index(k)
              ja(indexj) = k
              ja(indexk) = j
              index(j) = indexk
c
              ltemp = a(indexj+lshift)
              utemp = a(indexj)
              a(indexj+lshift) = a(indexk+lshift)
              a(indexj) = a(indexk)
              a(indexk+lshift) = ltemp
              a(indexk) = utemp
c
   40         k = list(k)
   50   continue
c
        return
c
c       error return
c
  100   ierr = i
        return
        end
        subroutine fillin (n, ja, jl, list)
           integer ja(*), jl(*), list(*), uptr
c
c       compute jl
c
c       input: n, ja
c       output: jl
c       workspace: list(n)
c
c       remarks:
c
c           1. for convenience, jl(n) is the size of the strict
c              upper triangular part of the factored matrix
c
c           2. the row indices for each column in ja are assumed
c              to be in decreasing order
c
c       initialize
c
        do 10 i = 1, n
           jl(i) = 0
   10      list(i) = 0
c
c       the main loop
c
        jl(n+1) = n + 2
        do 50 i = 1, n
           if (ja(i) .ge. ja(i+1))  go to 50
c
c       loop over seed indices in decreasing order
c
           list(i) = i
           length = 1
           do 30 iseed = ja(i), ja(i+1)-1
              k = ja(iseed)
              uptr = jl(iseed-1)
c 
c       add a new entry to list
c
   20         list(k) = list(i)
              list(i) = k
              uptr = uptr + 1
              length = length + 1
              if (jl(k) .eq. 0) jl(k) = i
              k = jl(k)
              if (list(k) .eq. 0) go to 20
c
   30         jl(iseed) = uptr
c
c       clean up loop for list array
c
           k = i
           do 40 j = 1, length
              ksave = k
              k = list(k)
   40         list(ksave) = 0
c
   50   continue
c
c       compute size of upper triangle
c
        lenjl = n + 1 + ja(n+1) - ja(1)
        jl(n) = jl(lenjl) - jl(n+1)
c
        return
        end
        subroutine factor(n,isym,ja,a,jl,u,index,expres,ierr)
           integer ja(*), jl(*), index(*), ashift, ushift ,
     1        expres(*), expk
           real a(*), u(*), dsum, usum, lsum
c
c       compute u
c
c       input: n, isym, ja, a, jl
c       output: u, ierr
c       workspace: index(n), expres(n)
c
c       remarks:
c
c           1. isym = 0 means nonsymmetric storage for a and u
c              isym = 1 means symmetric storage for a and u
c
c           2. ierr = 0 is a normal return
c              ierr = i means the i-th diagonal element of the
c                       factored matrix was zero
c
c           3. for symmetric matrices, the computation of lsum in the
c              do 60 loop, and the computation of u(indexk+ushift) in
c              the do 80 loop are redundant
c
c           4. the express vector modification was added in September, 
c              1986.
c
        ierr = 0
c
c       compute offsets for a and u
c
        if (isym .eq. 1) then
           ashift = 0
           ushift = 0
           lenu = n + 1 + jl(n)
        else
           ashift = ja(n+1) - ja(1)
           ushift = jl(n)
           lenu = n + 1 + 2 * jl(n)
        end if
c
c       initialize
c
        do 10 i = 1, lenu
   10      u(i) = 0.0e0
        expk = n
        do 20 i = n, 1, -1
           index(i) = 0
           if (jl(i) .eq. i+1) then
              expres(i) = expk
           else
              expres(i) = 0
              expk = i
           end if
   20   continue
        expres(n) = 0
c
c       the main loop
c
        do 100 i = 1, n
           dsum = 0.0
           if (ja(i) .ge. ja(i+1)) go to 90
c
c      loop over seed indices for column i in increasing order
c
           do 60 iseed = ja(i+1)-1, ja(i), -1
              k = ja(iseed)
c
c       move off diagonal entries from a to u
c
              u(jl(iseed-1)) = a(iseed)
              u(jl(iseed-1)+ushift) = a(iseed+ashift)
c
              do 50 indexk = jl(iseed-1), jl(iseed)-1
                 index(k) = indexk
                 expk = expres(k)
                 if (expk .gt. 0) then
                    if (expres(expk) .eq. 0) expres(expk) = -k
                 end if
                 if (ja(k) .ge. ja(k+1))  go to 50
c
c       loop over seed indices for column k in decreasing order
c
                 lsum = u(indexk+ushift)
                 usum = u(indexk)
c
                 do 40 kseed = ja(k), ja(k+1)-1
                    j = ja(kseed)
                    jstart = jl(kseed-1)
c
c       find start index for inner product loop
c
   25               if (jstart .ge. jl(kseed)) go to 40
                    if (index(j) .eq. 0) then
                       if (expres(j) .gt. 0) then
                          jold = j
                          j = expres(j) 
                          if (expres(j) .lt. 0) j = -expres(j)
                          jstart = jstart + j - jold
                       else
                          j = jl(j)
                          jstart = jstart + 1
                       end if
                       go to 25
                    end if
c
c       form inner products
c
                    do 30 indexj = jstart, jl(kseed)-1
                       lsum = lsum - u(indexj) * u(index(j)+ushift)
                       usum = usum - u(indexj+ushift) * u(index(j))
   30                  j = jl(j)
   40            continue
c
c       update u(indexk) and u(indexk+ushift)
c
                 u(indexk+ushift) = lsum
                 u(indexk) = usum
   50            k = jl(k)
   60      continue
c
c        clean up loop for index array, compute diagonal
c
           do 80 iseed = ja(i), ja(i+1)-1
              k = ja(iseed)
              do 70 indexk = jl(iseed-1), jl(iseed)-1
c
                 usave = u(indexk+ushift)
                 u(indexk) = u(indexk) * u(k)
                 u(indexk+ushift) = usave * u(k)
                 dsum = dsum + u(indexk) * usave
c
                 index(k) = 0
                 expk = expres(k)
                 if (expk .gt. 0) expres(expk) = 0
c
   70            k = jl(k)
   80      continue
c
c       compute the diagonal
c
   90      u(i) = a(i) - dsum
           if (u(i) .eq. 0.0e0)  go to 200
           u(i) = 1.0e0 / u(i)
  100   continue
c
        return
c
c       error return
c
  200   ierr = i
        return
        end

        subroutine solve (n, isym, ja, jl, u, x, b)
           integer ja(*), jl(*), ushift, lshift
           real u(*), x(*), b(*)
c
c       compute the solution x of a * x = b
c
c       input: n, isym, ja, jl, u, b
c       output: x
c
c       remark: isym = 1 means symmetric storage is used for u
c               isym = 0 means nonsymmetric storage is used for u
c               isym = -1 means nonsymmetric storage is used, and the
c                        problem is to be solved with the transpose of a
c
c       compute offsets for u
c
        if (isym .eq. 0) then
           lshift = jl(n)
        else
           lshift = 0
        end if
        if (isym .eq. -1) then
           ushift = jl(n)
        else
           ushift = 0
        end if
c
c       lower triangular system
c
        do 30 i = 1, n
           sum = 0.0e0
           if (ja(i) .ge. ja(i+1)) go to 30
           do 20 iseed = ja(i), ja(i+1)-1
              k = ja(iseed)
              do 10 j = jl(iseed-1), jl(iseed)-1
                 sum = sum + u(j+lshift) * x(k)
   10            k = jl(k)
   20      continue
   30      x(i) = b(i) - sum
c
c       diagonal system
c
        do 40 i = 1, n
   40      x(i) = u(i) * x(i)
c
c       upper triangular system
c
        do 70 i = n, 1, -1
           if (ja(i) .ge. ja(i+1)) go to 70
           do 60 iseed = ja(i), ja(i+1)-1
              k = ja(iseed)
              do 50 j = jl(iseed-1), jl(iseed)-1
                 x(k) = x(k) - u(j+ushift) * x(i)
   50            k = jl(k)
   60      continue
   70   continue
c
        return
        end