      subroutine psolve(u,a,b,dr,m,xm,s,ld,ig,idl,nst,nrs,afac,solv,
     &                  dyn,c1,ipd,rnorm,aengy,ifl)

c     Purpose:  Frontal assembly and solution of equations
c               N.B. interface for psolve same as for profile solutions

c     Inputs:
c        u(*)      - Current nodal solution
c        a(*)      - Storage for frontal tangent matrix
c        b(*)      - Frontal right hand side/solution vector
c        dr(*)     - Residual
c        m(*)      - Buffer for front equations eliminated
c        xm(*)     - Mass matrix (diagonal only)
c        s(nst,*)  - Element tangent matrix
c        ld(nst)   - Element equation numbers
c        ig(*)     - Pointers to end of frontal matrix row/columns
c        idl(*)    - Element ordering for frontal solution
c        nst       - Size of element arrays
c        nrs       - Number of right hand sides to solve
c        afac      - Flag, factor matrix if true
c        solv      - Flag, perform forward/backward solution if true
c        dyn       - Flag, transient solution if true
c        c1        - Factor to add mass matrix to tangent
c        ipd       - Precision of "real" variables

c     Outputs:
c        dr(*)     - Solution increment
c        rnorm     - Residual norm
c        aengy     - Energy of solution increment
c        ifl       - File name number for frontal equation storage on disk

      implicit  none

      integer   maxa
      parameter(maxa = 32000)

      logical   afac,solv,dyn,fa
      integer   nst,nrs,ipd,ifl,ihsiz
      integer   i,ie,ii, jj, k, n,ne,np,nal,nfrt
      real*8    aengy, c1, r,rxm
      real*8    rnorm

      integer   m(*),ld(*),ig(*),idl(*)
      real*8    u(*),a(*),b(*),dr(*),xm(*),s(*)

      integer         numnp,numel,nummat,nen,neq
      common /cdata/  numnp,numel,nummat,nen,neq

      logical         fl    ,pfr
      common /fdata/  fl(11),pfr

      integer         maxf
      common /frdata/ maxf

      integer         ior,iow
      common /iofile/ ior,iow

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

      character*12    tfile
      common /temfl1/ tfile(6)

      integer         itrec   ,nw1,nw2
      common /temfl2/ itrec(4),nw1,nw2

      integer         ihfac,ibuf
      common /temfl3/ ihfac,ibuf

      data fa/.false./

c     Record length factors (ihsiz = # int words; ihfac = # bytes/int wd)

      ipd   = 2
      ihfac = 4
      ihsiz = maxa
      if(.not.afac) go to 400

c     Set control data and zero

      if(fl(4)) then
        nal= 1 + (maxf*(maxf+3))/2
        if((nal+maxf)*ipd.gt.ihsiz) stop 'front too large'
        ibuf = (min(ihsiz - nal,(maxf+2)*neq))*ipd
        itrec(1) = ibuf*ihfac + ihfac
c       open (4,file=tfile(1),status='new',access='direct',
c    1        form='unformatted',recl=itrec(1))
c       close(4)
        fl(4) = fa
      endif
      call pconsd(a,maxf*(maxf+1)/2,0.0d0)
      call pconsd(b,maxf,0.0d0)
      ig(1)=1
      ig(maxf+1)=0

      do 100 n=2,maxf
       ig(n)=ig(n-1)+n
       ig(maxf+n)=0
100   continue

c     Begin loop through elements to perform front solution

      np  = 0
      nv  = 0
      nfrt= 0
      dimx  = 0.0
      dimn  = 0.0
      rnorm = 0.0

c     Frontal elimination program

      if(ior.lt.0) write(*,2000)

      do 320 n = 1,numel

c       Pick next element

        ne = idl(n)

c       Compute element arrays

        call formfe(u,dr,.true.,solv,fa,fa,3,ne,ne,1)
        if(ior.lt.0 .and. mod(n,20).eq.0) write(*,2001) n,nfrt

c       Assemble element and determine eliminations

        call pfrtas(a,s,ld,ig,ig(maxf+1),ie,nfrt,nst)

        if(ie.ne.0) then

c         Eliminate equations

          do 310 i=1,ie
            k=ld(i)
            jj=-ig(maxf+k)
            if(solv) then
              rnorm = rnorm + dr(jj)**2
              r=b(k)+dr(jj)
              dr(jj) = r
            endif
            if(dyn) then
             ii = ig(k)
             a(ii) = a(ii) + c1*xm(jj)
            endif
            if(np+8+nfrt*ipd.gt.ibuf) then
             call pbuff(m,ibuf,np,nv,2,ifl)
            endif
            m(np+1) = nfrt
            m(np+2) = k
            m(np+3) = jj
            call pfrtd(a,b,r,ig,ig(maxf+1),solv,m(np+5),nfrt,k)

c           Align last word of buffer for backsubstitution reads

            np = np + 8 + nfrt*ipd
            m(np) = nfrt
            nfrt = nfrt - 1
310       continue
        endif

c     End of forward elimination and triangularization

320   continue

      rnorm = sqrt(rnorm)
      rxm   = 0.0
      if(dimn.ne.0.0d0) rxm = dimx/dimn
      if(ior.lt.0) write(*,2002) dimx,dimn,rxm
      write(iow,2002) dimx,dimn,rxm

c     Clear buffer one last time

      if(np.gt.0) call pbuff(m,ibuf,np,nv,2,ifl)

      go to 500

c     Forward and back substitutions for solution

400   if(solv) call pfrtfw(b,dr,m,ipd,ibuf,maxf,nv,neq,nrs,ifl)
500   if(solv) call pfrtbk(b,dr,m,ipd,ibuf,maxf,nv,neq,nrs,aengy,ifl)

2000  format('   Solution status'/' ')
2001  format('+  ->',i4,' elements completed.  Current front width is',
     &           i4)
2002  format(5x,'Condition check: D-max',1p1e11.4,'; D-min',1p1e11.4,
     &          '; Ratio',1p1e11.4)

      end 

      subroutine pbuff(m,ibuf,ilast,nv,is,ifl)

c     Purpose: Input/output routine for frontal program

c     Inputs:
c        m(*)     - Buffer storage
c        ibuf     - Buffer length (in integer words)
c        ilast    - Uset to go backward through data files
c        nv       - Record number for I/O
c        is       - Switch: =1 for read; =2 for write
c        ifl      - File name number for frontal equation storage on disk

      implicit none

      character*12 fname
      integer  ibuf,ilast,nv,is,ifl
      integer  m(ibuf)

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

      character*12    tfile
      common /temfl1/ tfile(6)

      integer         itrec   ,nw1,nw2
      common /temfl2/ itrec(4),nw1,nw2

      integer         i,j

      if(is.eq.2) nv = nv + 1
      if(nv.lt.10) then
        write(fname,'(a10,i1)') 'Frontal.00',nv
      elseif(nv.lt.100) then
        write(fname,'( a9,i2)') 'Frontal.0',nv
      else
        write(fname,'( a8,i3)') 'Frontal.',nv
      endif
      open(ifl,file=fname,form='unformatted',status='unknown')

c     Read record 'nv' from the file

      if(is.eq.1) then

        do i = 1,ibuf,8000
          read(ifl,err=901)  ilast,(m(j),j=i,min(ibuf,i+7999))
        end do

c     Write record 'nv' from the file

      elseif(is.eq.2) then

        if(nv.eq.1) npl = ilast
        do i = 1,ibuf,8000
          write(ifl,err=901) ilast,(m(j),j=i,min(ibuf,i+7999))
        end do
        ilast = 0

      endif
      close(ifl)
      return

c     Error messages

901   call pend('PBUFF ')
      stop

4000  format(' **ERROR** PBUFF records do not match problem?')

      end 

      subroutine pfrtas(a,s,ld,ig,lg,ie,nfrt,nst)

c     Purpose: Assembly and elimination determination for frontal program

c     Inputs:
c        s(nst,*) - Element matrix to assemble
c        ld(*)    - Equation numbers for element (negative tag to elim)
c        ig(*)    - Pointers to end of columns in frontal matrix
c        lg(*)    - Front/global equation numbers
c        nfrt     - Size of front before eliminations
c        nst      - Size of element arrays

c     Outputs:
c        a(*)     - Assembled frontal matrix
c        ld(*)    - Locations to eliminate
c        ie       - Number of equations to eliminate after assembly step

      implicit none

      integer  ie, nfrt, nst
      integer  i,ii, j, k, l, m

      real*8   a(*),s(nst,*)
      integer  ld(*),ig(*),lg(*)

c     Convert ld to front order

      do 203 j=1,nst
        ii=abs(ld(j))
        i = 0

c       Check if ii is already in list

        if(ii.ne.0) then
          if(nfrt.ne.0) then
            do 200 i=1,nfrt
              if(ii.eq.abs(lg(i))) go to 202
200         continue
          endif

c         Assign ii to next available entry and increase front width

          nfrt = nfrt + 1
          i = nfrt

c         Replace destination value by new value

202       lg(i)=ld(j)

c         Set ld for assembly

          ld(j)=i
        endif
203   continue

c     Assemble element into front

      do 205 j = 1,nst
        k = abs(ld(j))
        if(k.gt.0) then
          l = ig(k) - k
          do 204 i = 1,nst
            m = abs(ld(i))
            if(m.gt.0 .and. m.le.k) a(m+l) = a(m+l) + s(i,j)
204       continue
        endif
205   continue

c     Set up equations to be eliminated

      ie=0
      do 206 i = nfrt,1,-1
        if(lg(i).lt.0) then
          ie = ie + 1
          ld(ie) = i
        endif
206   continue

      end 

      subroutine pfrtbk(b,dr,m,ipd,ibuf,maxf,nv,neq,nev,aengy,ifl)

c     Purpose: Backsubstitution for frontal solution

c     Inputs:
c        b(*)     - Frontal right hand side/solution vector
c        dr(*)    - Forward reduced solution vector
c        m(*)     - Buffer storage array
c        ipd      - Real*8 precision
c        ibuf     - Buffer length
c        maxf     - Maximum front allowed
c        nv       - Number of buffer blocks
c        neq      - Number of equations to solve
c        nev      - Number of right hand sides
c        ifl      - File name number for frontal equation storage on disk

c     Outputs:
c        dr(*)    - Solution incrrement
c        aengy    - Energy of solution increment

      implicit none

      integer  ipd,ibuf,maxf,nv,neq,nev,ifl
      integer  i, jj, k, n,np,nfrt
      real*8   aengy

      real*8   b(maxf,*),dr(neq,*)
      integer  m(*)

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

      call pconsd(b,maxf*nev,0.0d0)
      aengy = 0.0

c     Recover block

      np = npl
      do 503 n = nv,1,-1
        if(nv.gt.1) call pbuff(m,ibuf,np,n,1,ifl)
501     nfrt = m(np)
        np = np - 8 - nfrt*ipd
        k  = m(np+2)
        jj = m(np+3)
        do 502 i = 1,nev
          call pfrtb(b(1,i),dr(1,i),nfrt,k,jj,m(np+5),aengy)
502     continue
        if(np.gt.0) go to 501
503   continue

      end 

      subroutine pfrtfw(b,dr,m,ipd,ibuf,maxf,nv,neq,nev,ifl)

c     Purpose: Forward solution for frontal resolutions

c     Inputs:
c        b(*)     - Frontal right hand side/solution vector
c        dr(*)    - Residual vector
c        m(*)     - Buffer storage array
c        ipd      - Real*8 precision
c        ibuf     - Buffer length
c        maxf     - Maximum front allowed
c        nv       - Number of buffer blocks
c        neq      - Number of equations to solve
c        nev      - Number of right hand sides
c        ifl      - File name number for frontal equation storage on disk

c     Outputs:
c        dr(*)    - Forward reduced solution

      implicit none

      integer  ipd,ibuf,maxf,nv,neq,nev,ifl
      integer  i,ilast, jj, k, n,np,nfrt

      real*8   b(maxf,1),dr(neq,1)

      integer  m(*)

      call pconsd(b,maxf*nev,0.0d0)

      do 403 n = 1,nv
        call pbuff(m,ibuf,ilast,n,1,ifl)
        np = 1
401     nfrt = m(np)
        k    = m(np+1)
        jj   = m(np+2)
        do 402 i = 1,nev
          call pfrtf(b(1,i),dr(1,i),nfrt,k,jj,m(np+4))
402     continue

c       Align last word of buffer for backsubstitution reads

        np = np + 8 + m(np)*ipd
        if(np.lt.ilast) go to 401
403   continue

      end 

      subroutine prefrt(il,idl,ix,maxf,ndf,nen,nen1,numel,numnp)

c     Purpose: Prefrontal routine to flag last occurance of nodes

c     Inputs:
c        il(*)      - Working vector
c        idl(*)     - Elemenet elimination order
c        ix(nen1,*) - Element nodal connections
c        ndf        - Number dof/node
c        nen        - Maximum number of nodes/element
c        nen1       - Dimension of ix array
c        numel      - Number of elements in mesh
c        numnp      - Number of nodes in mesh

c     Outputs:
c        ix(nen1,*) - Element nodal connections (tag last occurances)
c        maxf       - Maximum front estimate

      implicit  none

      integer   maxf,ndf,nen,nen1,numel,numnp
      integer   i,ii, jj, n,nu,nowf

      integer   il(*),idl(*),ix(nen1,*)

      integer         ior,iow
      common /iofile/ ior,iow

c     Preset check array

      do 100 n=1,numnp
        il(n)=0
100   continue

c     Set last occurance of nodes

      do 102 nu=numel,1,-1
        n = idl(nu)
        do 101 i=nen,1,-1
          ii=abs(ix(i,n))
          if((ii.ne.0).and.(il(ii).eq.0)) then
            il(ii)=n
            ii=-ii
          endif
          ix(i,n)=ii
101     continue
102   continue

c     Get estimate to maximum frontwith

      maxf=0
      nowf=0
      do 107 nu=1,numel
        n = idl(nu)
        do 105 i=1,nen
          ii=ix(i,n)
          if(ii.ne.0) then
            jj=abs(ii)
            if(il(jj).ne.0) nowf=nowf+ndf
            maxf=max(maxf,nowf)
            il(jj)=0
          endif
105     continue
        do 106 i = 1,nen
          if(ix(i,n).lt.0) nowf = max(0,nowf-ndf)
106     continue
107   continue
      if(ior.lt.0) write(*,3000) maxf

3000  format('   -> Maximum Front Estimate =',i4,' d.o.f.')

      end 
c
      subroutine profil (jd,idl,id,ix,ndf,nen1)

c     Purpose: Compute front profile of global arrays

c     Inputs:
c        jd(*)      - Working array
c        idl(*)     - Element elimination order
c        id(ndf,*)  - Boundary condition array
c        ix(nen1,*) - Element nodal connection list
c        ndf        - Number dof/node
c        nen1       - Dimension of ix array

c     Outputs:
c        id(ndf,*)  - Equation numbers for each dof
c        ix(nen1,*) - Element nodal connection list tagged for frontal
c                     eliminations.

      implicit  none

      integer   ndf,nen1, nneq
      integer   j, n

      integer   jd(*),idl(*),id(*),ix(nen1,*)

      integer        numnp,numel,nummat,nen,neq
      common /cdata/ numnp,numel,nummat,nen,neq

      integer         ior,iow
      common /iofile/ ior,iow

      integer         maxf
      common /frdata/ maxf

c     Set up the equation numbers

      neq = 0
      nneq = ndf*numnp
      do 10 n = 1,nneq
        j = id(n)
        if(j.eq.0) then
          neq = neq + 1
          id(n) = neq
        else
          id(n) = 0
        endif
10    continue

c     Set default element elimination order
c     N.B. Can replace with optimized ordering.

      do 11 n = 1,numel
        idl(n) = n
11    continue

c     Compute front width

      call prefrt(jd,idl,ix,maxf,ndf,nen,nen1,numel,numnp)

c     Max front width must produce a triangle which fits in adata
c     plus some buffer space for reduced equations.

      if(maxf.gt.150) then
        write(*,2001) maxf
        if(ior.lt.0) write(iow,2001) maxf
      endif

2001  format(' *ERROR* Front requires too much storage =',i8)

      end

      subroutine rsolve(b,dr,m,ipd,ipr,maxf,nv,neq,nev,aengy,ifl)

c     Purpose: Resolution for frontal solution

c     Inputs:
c        b(*)      - Frontal RHS/solution
c        dr(neq,*) - Residual
c        m(*)      - Buffer for fronal equations
c        ipd       - Precision of 'real*8' variablels
c        ipr       - Precision of 'real*4' variablels
c        maxf      - Maximum front estimate
c        nv        - Number of frontal buffer blocks
c        neq       - Number of equations to solve
c        nev       - Number of RHS
c        ifl       - File name number for frontal equation storage on disk
        
c     Outputs:
c        dr(neq,*) - Residual
c        aengy     - Energy of solution

      implicit  none

      integer   ipd,ipr,maxf,nv,neq,nev,ifl
      real*8    aengy

      real*8    b(maxf,*),dr(neq,*)

      integer   m(*)

      integer         itrec   ,nw1,nw2
      common /temfl2/ itrec(4),nw1,nw2

      integer         ihfac,ibuf
      common /temfl3/ ihfac,ibuf

      call pfrtfw(b,dr,m,ipd,ibuf,maxf,nv,neq,nev,ifl)
      call pfrtbk(b,dr,m,ipd,ibuf,maxf,nv,neq,nev,aengy,ifl)

      end

      subroutine delfrt()

      implicit none

      character fname*12
      logical   lfil
      integer   n

c     Delete all files from frontal solver

      do n = 1,9999

c       Set name

        if(n.lt.10) then
          write(fname,'(a10,i1)') 'Frontal.00',n
        elseif(n.lt.100) then
          write(fname,'( a9,i2)') 'Frontal.0',n
        else
          write(fname,'( a8,i3)') 'Frontal.',n
        endif

c       Check and delete

        inquire(file=fname,exist=lfil)
        if(lfil) then
          open (4,file=fname,status='old')
          close(4,status='delete')
        else
          return
        endif
      end do

      end
