      subroutine pfrtb(b,dr,nfrt,k,jj,eq,aengy)

c     Purpose: Backsubstitution macro for frontal program

c     Inputs:
c        b(*)     - Front reduced equations
c        nfrt     - Size of current front
c        k        - Location of diagonal
c        jj       - Equation number in global equations
c        eq       - Current equation

c     Outputs
c        b(*)     - Remaining front reduced equations
c        dr(*)    - Solution increment
c        aengy    - Energy increment

      implicit none

	  integer k,kk, jj, nfrt
      real*8 b(*),dr(*),eq(*),dot,aengy

c     Expand b array

      kk = nfrt
500   if(kk.le.k) go to 501
      b(kk) = b(kk-1)
      kk = kk - 1
      go to 500
501   b(k) = 0.0

c     Extract pivot and solve, also compute energy

      aengy = aengy + dr(jj)**2/eq(k)
      dr(jj)=(dr(jj)-dot(eq,b,nfrt))/eq(k)
      b(k) = dr(jj)

      end 

      subroutine pfrtd(a,b,r,ig,lg,solv,eq,nfrt,k)

c     Purpose: Triangular decomposition for frontal program

c     Inputs:
c        a(*)      - Frontal matrix before reduction
c        b(*)      - RHS of frontal equ;ations
c        r         - Residual for current equation
c        ig(*)     - Pointer to diagonals in a(*)
c        lg(*)     - Global equation numbers for each frontal equation
c        solv      - Flag, solve equations if true; otherwise just factor

c     Outputs:
c        eq(*)     - Eliminated equation
c        nfrt      - Size of front
c        k         - Location of diagonal in eq(*)

      implicit none

      logical  solv
	  integer  j,jj, k,kk,kl,km,kp, l, nfrt
      integer  ig(*),lg(*)
      real*8   a(*),b(*),eq(*),pivot,term,r

	  real*8         dimx,dimn
      integer                  nv,npl
      common /nfrta/ dimx,dimn,nv,npl

c     Extract equation

      do 301 j=1,nfrt
        l = max(j,k)
        l=ig(l)-l+min(j,k)
        eq(j)=a(l)
301   continue

      pivot=eq(k)

      if(pivot.eq.0.0d0) then
        call pconsd(eq,nfrt,0.0d0)
        write(*,2000)
        pivot = 1.0
        eq(k) = 1.0
      else
        dimx = max(dimx,abs(pivot))
        if(dimn.eq.0.0d0) dimn = abs(pivot)
        dimn = min(dimn,abs(pivot))
      endif
      km = k - 1
      kp = k + 1
      kk = 1
      if(km.gt.0) then
        do 302 jj = 1,km
          if(eq(jj).ne.0.0d0) then
            term = -eq(jj)/pivot
            if(solv) b(jj) = b(jj) + term*r
            call saxpb(eq,a(kk),term,jj,a(kk))
          endif
          kk = ig(jj) + 1
302     continue
      endif
      kk = ig(k) + 1
      if(kp.le.nfrt) then
        do 303 jj = kp,nfrt
          lg(jj-1) = lg(jj)
          kl = kk - jj
          term = -eq(jj)/pivot
          if(solv) b(jj-1) = b(jj) + term*r
          call saxpb(eq(kp),a(kk+k),term,jj-k,a(kl+k))
          if(km.gt.0) then
            kl = kl + 1
            call saxpb(eq,a(kk),term,km, a(kl))
          endif
          kk = ig(jj) + 1
303     continue
      endif
      call pconsd(a(kk-nfrt),nfrt,0.0d0)
      b(nfrt)  = 0.0
      lg(nfrt) = 0

2000  format(' WARNING -- Zero pivot, check boundary codes.'/
     &       '            Pivot set to 1.0 and solution continued.')

      end 

      subroutine pfrtf(b,dr,nfrt,k,jj,eq)

c     Purpose: Forward elimination macro for front program

c     Inputs:
c        b(*)     - Front residual
c        nfrt     - Size of current front
c        k        - Location of diagonal
c        jj       - Equation number in global equations
c        eq       - Current equation

c     Outputs
c        b(*)     - Reduced front equation

      implicit none

	  integer  nfrt, ii, jj, k,km,kp
      real*8   b(*),dr(*),eq(*),r

      dr(jj) = dr(jj) + b(k)
      r = dr(jj)/eq(k)
      km = k - 1
      kp = k + 1

      if(km.gt.0) then
        do 402 ii = 1,km
          b(ii) = b(ii) - eq(ii)*r
402     continue
      endif

      if(kp.le.nfrt) then
        do 404 ii = kp,nfrt
          b(ii-1) = b(ii) - eq(ii)*r
404     continue
      endif

      b(nfrt) = 0.0

      end 
