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

c     Purpose:  Active column assembly and solution of equations
c               N.B. interface for psolve same as for frontal solutions

c     Inputs:
c        u(*)      - Current nodal solution
c        a(*)      - Storage for tangent matrix
c        b(*)      - Not used for profile solution
c        dr(*)     - Residual
c        m(*)      - Not used for profile solution
c        xm(*)     - Mass matrix (diagonal only)
c        s(nst,*)  - Element tangent matrix
c        ld(nst)   - Element equation numbers
c        ig(*)     - Pointers to end of profile row/columns
c        idl(*)    - Not used for profile solution
c        nst       - Size of element arrays
c        nrs       -
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       - Not used for profile solution

      implicit  none

      integer    maxa
      parameter (maxa=32000)

      logical   afac,afl,solv,dyn,fa
      integer   nst,nrs,ipd,ifl
      integer   n,ne,nep
      integer   ibuf,ihsiz,ihfac,iz
      real*8    aengy,c1,dot,rnorm

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

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

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

      integer         ior,iow
      common /iofile/ ior,iow

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

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

      data   fa /.false./

c     Form and assemble matrix

c     Record length factors (may be too long for some machines)

      ihfac = 8 ! record length factor (8)
      ihsiz = maxa

      if(afac) then
        if(fl(6)) then
          nep = neq + ig(neq)
          afl = .true.
        else
          nep = neq
          afl = .false.
        endif
        if(fl(3).or.fl(4)) then
          ibuf = ig(neq)+neq
          if(fl(3)) then
            ibuf  = ibuf + ig(neq)
            fl(3) = fa
          endif
          if(ibuf.gt.ihsiz) stop 'profile too large'
          iz       = ibuf
          itrec(1) = iz*ihfac
          open (4,file=tfile(1),status='new',access='direct',
     &          form='unformatted',recl=itrec(1))
          close(4)
          fl(4) = fa
        endif

        call pconsd(a,ibuf,0.0d0)

c       Modify tangent form lumped mass effects

        if(dyn) then
          do 310 n = 1,neq
            a(n)  =  c1*xm(n)
310       continue
        endif

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

c       Compute and assemble element arrays

        do 320 n = 1,numel
          ne = n
          call formfe(u,dr,.true.,solv,fa,fa,3,ne,ne,1)
          if(ior.lt.0 .and. mod(n,50).eq.0) write(*,2001) n
          call dasbly(s,s,ld,ig,nst,afl,afac,fa,dr,a(nep+1),a(neq+1),a)
320     continue

        rnorm = sqrt(dot(dr,dr,neq))

        call datri(a(nep+1),a(neq+1),a,ig,neq,afl)
        call phstio(4,1,a,ibuf,2,tfile(1),itrec(1))

      endif

      if(solv) then
        if(.not.afac) call phstio(4,1,a,ibuf,1,tfile(1),itrec(1))
        do 330 n = 1,nrs
          ne = (n-1)*neq + 1
          call dasol(a(nep+1),a(neq+1),a,dr(ne),ig,neq, aengy)
330     continue
      endif

2000  format('+  Solution status'/' ')
2001  format('+  ->',i4,' Elements completed.')

      end

      subroutine dasbly(s,p,ld,jp,ns,alfl,aufl,bfl,  b,al,au,ad)

c     Purpose:  Assemble symmetric or unsymmetric arrays for 'dasol'

c     Inputs:
c        s(ns,ns)   - element matrix
c        p(ns)      - Element vector
c        ld(ns)     - Equation numbers for elements
c        jp(*)      - Pointer to end of row/columns in profile
c        alfl       - Flag, unsymmetric if true
c        aufl       - Flag, assemble tangent if true
c        bfl        - Flag, assembel vector if true

c     Outputs:
c        b(*)       - Assembled vector
c        al(*)      - Assembled lower part of matrix
c        au(*)      - Assembled upper part of matrix
c        ad(*)      - Assembled matrix diagonal

      implicit none
      logical  alfl,aufl,bfl
      integer  ns
      integer  i,ii, j,jc,jj

      integer  ld(ns),jp(*)
      real*8   al(*),au(*),ad(*),b(*),s(ns,ns),p(ns)

c     Loop through rows to perform assembly
      do 200 i = 1,ns
        ii = ld(i)
        if(ii.gt.0) then
          if(aufl) then

c          Loop through columns to perform assembly

            do 100 j = 1,ns
              if(ld(j).eq.ii) then
                ad(ii) = ad(ii) + s(i,j)
              elseif(ld(j).gt.ii) then
                jc = ld(j)
                jj = ii + jp(jc) - jc + 1
                au(jj) = au(jj) + s(i,j)
                if(alfl) al(jj) = al(jj) + s(j,i)
              endif
100         continue
          endif
          if(bfl) b(ii)  = b(ii)  + p(i)
        endif
200   continue

      end

      subroutine dasol(al,au,ad,b,jp,neq, energy)

c     Purpose:
c        Solution of symmetric equations stored in profile form
c        coefficient matrix must be decomposed into its triangular
c        factors using datri before using dasol.

c     Inputs:
c        al(*)   - Lower part of factored matrix
c        au(*)   - Upper part of factored matrix
c        ad(*)   - Reciprocal diagonal of factored matrix
c        b(*)    - Right hand side of equations
c        jp(*)   - Pointer to end of row/columns in profile
c        neq     - Number of equations

c     Outputs:
c        b(*)    - Solution increment
c        energy  - Engergy of incremental solution

      implicit none

      integer  neq
      integer  is, j,jh,jr
      real*8   energy,bd,dot

      integer   jp(*)
      real*8    al(*),au(*),ad(*),b(*)

      integer         ior,iow
      common /iofile/ ior,iow

c     Find first nonzero entry in right hand side

      energy = 0.0d0

      do 100 is = 1,neq
        if(b(is).ne.0.0d0) go to 200
100   continue

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

200   if(is.lt.neq) then

c     Reduce right hand side

        do 300 j = is+1,neq
          jr = jp(j-1)
          jh = jp(j) - jr
          if(jh.gt.0) then
            b(j) = b(j) - dot(al(jr+1),b(j-jh),jh)
          endif
300     continue
      endif

c     Multiply by inverse of diagonal elements

      do 400 j = is,neq
        bd = b(j)
        b(j) = b(j)*ad(j)
        energy = energy + bd*b(j)
400   continue

c     Backsubstitution

      if(neq.gt.1) then
        do 500 j = neq,2,-1
          jr = jp(j-1)
          jh = jp(j) - jr
          if(jh.gt.0) then
            call saxpb(au(jr+1),b(j-jh),-b(j),jh, b(j-jh))
          endif
500     continue
      endif

2000  format(' **WARNING** Zero right-hand-side vector')

      end

      subroutine datest(au,jh,daval)

c     Purpose: Test for rank of matrix with zero diagonals

c     Inputs:
c        au(*)      - Column for j-th equation
c        jh         - Length of column

c     Outputs:
c        daval      - Sum of absolute values of terms, if zero
c                     entire equation assumed to be zero and skipped.

      implicit none

      integer  j,jh
      real*8   au(jh),daval

      daval = 0.0d0
      do 100 j = 1,jh
         daval = daval + abs(au(j))
100   continue

      end
c
      subroutine datri(al,au,ad,jp,neq,flg)

c     Purpose: Triangular decomposition of a matrix stored in profile form

c     Inputs:
c        al(*)    - Lower part of matrix
c        au(*)    - Upper part of matrix
c        ad(*)    - Diagonal of matrix
c        jp(*)    - Pointer to end of row/columns
c        neq      - Number of equations to factor
c        flg      - Flag, equations unsymmetric if true

c     Outputs:
c        al(*)    - Lower part of factored matrix
c        au(*)    - Upper part of factored matrix
c        ad(*)    - Reciprocal diagonal of factored matrix

      implicit none

      logical  flg
      integer  neq
      integer  i,id,ie,ih,is,idh,ifig, j,jd,jh,jr,jrh
      real*8   dd,dimx,dimn,dfig,tol,daval,dot

      integer  jp(*)
      real*8   al(*),au(*),ad(*)

      integer         ior,iow
      common /iofile/ ior,iow

c     N.B.  tol should be set to approximate half-word precision.

      data tol/0.5d-07/

c     Set initial values for conditioning check

      dimx = 0.0d0
      dimn = 0.0d0
      do 50 j = 1,neq
        dimn = max(dimn,abs(ad(j)))
50    continue
      dfig = 0.0d0

c     Loop through columns to perform triangular decomposition

      jd = 1
      do 200 j = 1,neq
        jr = jd + 1
        jd = jp(j)
        jh = jd - jr
        if(jh.gt.0) then
          is = j - jh
          ie = j - 1

c         If diagonal is zero compute a norm for singularity test

          if(ad(j).eq.0.0d0) call datest(au(jr),jh,daval)
          do 100 i = is,ie
            jr = jr + 1
            id = jp(i)
            ih = min(id-jp(i-1),i-is+1)
            if(ih.gt.0) then
              jrh = jr - ih
              idh = id - ih + 1
              au(jr) = au(jr) - dot(au(jrh),al(idh),ih)
              if(flg) al(jr) = al(jr) - dot(al(jrh),au(idh),ih)
            endif
100       continue
        endif

c       Reduce diagonal

        dd = ad(j)
        if(jh.ge.0) then
          jr = jd - jh
          jrh = j - jh - 1
          call dredu(al(jr),au(jr),ad(jrh),jh+1,flg  ,ad(j))

c         Check for possible errors and print warnings

          if(abs(ad(j)).lt.tol*abs(dd))  write(iow,2000) j
          if(dd.lt.0.d0.and.ad(j).gt.0.d0) write(iow,2001) j
          if(dd.gt.0.d0.and.ad(j).lt.0.d0) write(iow,2001) j
          if(ad(j) .eq.  0.d0)             write(iow,2002) j
          if(dd.eq.0.d0.and.jh.gt.0) then
            if(abs(ad(j)).lt.tol*daval)   write(iow,2003) j
          endif
          if(ior.lt.0) then
            if(abs(ad(j)).lt.tol*abs(dd))  write(*,2000) j
            if(dd.lt.0.d0.and.ad(j).gt.0.d0) write(*,2001) j
            if(dd.gt.0.d0.and.ad(j).lt.0.d0) write(*,2001) j
            if(ad(j) .eq.  0.d0)             write(*,2002) j
            if(dd.eq.0.d0.and.jh.gt.0) then
              if(abs(ad(j)).lt.tol*daval)   write(*,2003) j
            endif
          endif
        endif

c       Store reciprocal of diagonal, compute condition checks

        if(ad(j).ne.0.d0) then
          dimx  = max(dimx,abs(ad(j)))
          dimn  = min(dimn,abs(ad(j)))
          dfig  = max(dfig,abs(dd/ad(j)))
          ad(j) = 1.d0/ad(j)
        endif
        if(ior.lt.0 .and. mod(j,50).eq.0) write(*,2005) j,neq
200   continue

c     Print conditioning information

      dd = 0.0d0
      if(dimn.ne.0.0d0) dd = dimx/dimn
      ifig = log10(dfig) + 0.6
      write(iow,2004) dimx,dimn,dd,ifig
      if(ior.lt.0) write(*,2004) dimx,dimn,dd,ifig

c     Formats

2000  format(' **WARNING** Loss of at least 7 digits in reducing',
     & ' diagonal of equation',i5)
2001  format(' **WARNING** Sign of diagonal changed when reducing',
     & ' equation',i5)
2002  format(' **WARNING** Reduced diagonal is zero for equation',i5)
2003  format(' **WARNING** Rank failure for zero unreduced',
     & ' diagonal in equation',i5)
2004  format(' Condition check: D-max',1pe11.4,'; D-min',1pe11.4,
     & '; Ratio',1pe11.4/' Maximum no. diagonal digits lost:',i3)
2005  format('+  ->',i5,' Equations of',i5,' reduced')

      end

      subroutine dredu(al,au,ad,jh,flg  ,dj)

c     Purpose: Reduce diagonal element in triangular decomposition

c     Inputs:
c        al(*)    - Loser row of equation
c        au(*)    - Upper column of equation
c        ad(*)    - Reciprocal diagonals for equations
c        jh       - length of row/column 
c        flg      - Flag, equations unsymmetric if true

c     Outputs:
c        dj       - Reduced diagonal

      implicit none

      logical  flg
      integer  j,jh

      real*8 al(jh),au(jh),ad(jh),dj,ud

c     Computation of column for unsymmetric matrices

      if(flg) then
        do 100 j = 1,jh
          au(j) = au(j)*ad(j)
          dj    = dj - al(j)*au(j)
          al(j) = al(j)*ad(j)
100   continue

c     Computation of column for symmetric matrices

      else
        do 200 j = 1,jh
          ud    = au(j)*ad(j)
          dj    = dj - au(j)*ud
          au(j) = ud
200     continue
      endif

      end

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

c     Purpose: Compute profile of global arrays

c     Inputs:
c        idl(ndf,*) - Temporary storage
c        id(*)      - Boundary condition codes
c        ix(nen1,*) - Element nodal connections
c        ndf        - Number dof/node
c        nen1       - Dimension of ix array

c     Outputs:
c        id(ndf,*)  - Equation numbers for each dof
c        jd(*)      - Profile for equations

      implicit none
      integer  ndf,nen1
      integer  i,ii, j,jj, mm, n,nad,nneq

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

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

      integer         maxf
      common /frdata/ maxf

      integer         ior,iow
      common /iofile/ ior,iow

c     Set up equation numbers

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

c     Compute column heights

      call pconsi(jd,neq,0)
      do 50 n = 1,numel
        mm  = 0
        nad = 0
        do 30 i = 1,nen
          ii = iabs(ix(i,n))
          if(ii.gt.0) then
            do 20 j = 1,ndf
              jj = id(j,ii)
              if(jj.gt.0) then
                if(mm.eq.0) mm = jj
                mm = min(mm,jj)
                nad = nad + 1
                idl(nad) = jj
              endif
20          continue
          endif
30      continue
        if(nad.gt.0) then
          do 40 i = 1,nad
            ii = idl(i)
            jj = jd(ii)
            jd(ii) = max(jj,ii-mm)
40        continue
        endif
50    continue

c     Compute diagonal pointers for profile

      nad = 0
      jd(1) = 0
      if(neq.gt.1) then
        do 60 n = 2,neq
          jd(n) = jd(n) + jd(n-1)
60      continue
        nad = jd(neq)
      endif

c     Set element search order to sequential

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

c     Equation summary

      maxf = 0
      mm   = 0
      if(neq.gt.0) mm = (nad+neq)/neq
      write(iow,2001) neq,numnp,mm,numel,nad,nummat
      if(ior.lt.0) write(*,2001) neq,numnp,mm,numel,nad,nummat

2001  format(/'   E q u a t i o n    P r o b l e m   S u m m a r y:'//
     &   5x,'Number of equations  =',i7,5x,'Number nodes      =',i5/
     &   5x,'Average col. height  =',i7,5x,'Number elements   =',i5/
     &   5x,'No. terms in profile =',i7,5x,'Number materials  =',i5/)

      end

      subroutine rsolve(b,dr,a,ipd,ipr,maxf,nv,neq,nev,engy,ifl)

c     Purpose: Resolution for profile solution

c     Inputs: 
c        b(*)      - Not used for profile solution
c        dr(neq,*) - Right hand side of equations
c        a(*)      - Profile storage for matrix
c        ipd       - Precision of "real*8" variables
c        ipr       - Precision of "real*4" variables
c        maxf      - Not used for profile solution
c        nv        - Not used for profile solution
c        neq       - Number of equations in matrix
c        nev       - Number of right hand side vectors to solve
c        ifl       - Not used for profile solution

c     Outputs:
c        dr(neq,*) - Solutions of equations
c        engy      - Energy of solution

      implicit none
      integer  ipd,ipr,maxf,nv,neq,nev,ifl
      integer  ne,n12
      real*8   a(*),b(*),dr(neq,*),engy

      real*8   dm
      real           rm
      integer              im
      common   dm(1),rm(1),im(1)

      n12 = neq*ipd - ipd - ipr + 1
      do 100 ne = 1,nev
        call dasol(a(neq+1),a(neq+1),a,dr(1,ne),im(n12),neq, engy)
100   continue

      end

      subroutine delfrt()
c     Dummy routine: not used for profile
      end
