c***************************** file: mg1.f *****************************
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine pltmg(vx,vy,xm,ym,itnode,ibndry,ja,a,ip,rp,sp,w,
     +      a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          itnode(5,*),ibndry(6,*),ip(100),ja(*)
            double precision
     +          vx(*),vy(*),xm(*),ym(*),w(*),rp(100),a(*)
            character*80
     +          sp(100)
            external a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy
c
c       user specified ip variables
c
        if(ip(6).lt.0.or.ip(6).gt.3) ip(6)=1
        if(ip(7).lt.-5.or.ip(7).gt.6) ip(7)=1
        if(ip(12).ne.1) ip(12)=0
        if(ip(8).ne.1) ip(8)=0
        ip(10)=max0(1,ip(10))
        ip(11)=max0(1,ip(11))
        rp(3)=dmax1(rp(3),0.0d0)
        ip(25)=0
        if(ip(6).ne.0) ip(24)=0
c
c       storage allocation
c
        if(ip(6).ne.0) then
            call stor(ip)
            if(ip(25).ne.0) go to 20
        endif
c
c       error flags
c
        if(itnode(3,1).eq.0) then
            ip(25)=25
            go to 20
        endif
c
c       array pointers...in the order that they
c       occur in the w array
c
        iuu=ip(90)
        itdof=ip(91)
        jtime=ip(92)
        jhist=ip(93)
        jpath=ip(94)
        ka=ip(95)
        jstat=ip(96)
        iee=ip(97)
        ipath=ip(98)
        iz=ip(99)
c
        ntf=ip(1)
        nvf=ip(2)
        nbf=ip(4)
        ispd=ip(8)
        iprob=ip(7)
        itask=ip(9)
        iord=ip(26)
        ndof=(iord+1)*(iord+2)/2
        lenw=ip(82)
        mpisw=ip(48)
        nproc=ip(49)
        irgn=ip(50)
        ndf=ip(5)
        maxd=ip(89)
        ibegin=iz
        iend=lenw
c
c       check for mpi status
c
        if(iprob.lt.0) then
            if(mpisw.ne.1) then
                ip(25)=48
                go to 20
            endif
            call exflag(ip(24))
            if(ip(24).ne.0) then
                ip(25)=24
                go to 20
            endif
        endif
c
        if(ip(6).ne.0) then
            call timer(w(jtime),-2)
            call hist2(w(jhist),rp,0,0)
            call updpth(w(jpath),1,1,rp)
            call pstat1(ntf,nproc,w(jstat),itnode,w(iee),0)
            call dschek(vx,vy,xm,ym,itnode,ibndry,ip,rp,sp,w)
            if(ip(25).ne.0) return
            rp(21)=rp(1)
            rp(31)=rp(1)
            rp(33)=1.0d0
            rp(45)=0.0d0
            rp(53)=1.0d0
            rp(59)=0.0d0
            rp(60)=0.0d0
            rp(64)=1.0d-3
            if(ip(7).eq.3.and.ip(9).lt.3) ip(9)=3
            if(ip(7).eq.4) ip(9)=8
            ip(6)=0
            ip(70)=0
c
c       setup itdof
c
            call memptr(isv,0,'mark',ibegin,iend,iflag)
            call memptr(iequv,nvf,'head',ibegin,iend,iflag)
            call memptr(itedge,3*ntf,'head',ibegin,iend,iflag)
            call memptr(ibedge,2*nbf,'head',ibegin,iend,iflag)
            ll=3*ntf+nbf+nvf
            call memptr(mark,ll,'head',ibegin,iend,iflag)
            call mkdof(ip,itnode,ibndry,w(itedge),w(ibedge),
     +          w(iequv),w(mark),ndof,w(itdof))
            call memptr(isv,0,'free',ibegin,iend,iflag)
c
            ndf=ip(5)
            maxd=ip(89)
c
            call gfinit(ip,maxd,w(iuu),w(iee))
        else
            call timer(w(jtime),-1)
        endif
c
        call gfptr(iprob,itask,maxd,iuu,iu0,iudot,iu0dot,
     +          ievr,ievl,ivx0,ivy0,ium,iuc,ngf,nef)
c
        call memptr(ibedge,2*nbf,'head',ibegin,iend,iflag)
c
        ndd=max0(0,ip(33))
        nvdd=max0(0,ip(71))
        jbegin=ibegin
        call lsptr(iprob,ispd,iord,ndf,ndd,nvdd,jbegin,jend,
     +      ihh,isu,ism,ia0,ih0,ig0,isu0,ism0,ibb,idd,
     1      ird,ipp,idl,ibdlwr,ibdupr,idu,idum,iduc)
        ll=jend-jbegin
        call memptr(ils,ll,'head',ibegin,iend,iflag)
        if(iprob.lt.0) then
            maxja0=9*iord*nvdd/2
            call memptr(iudd,11*nvdd,'head',ibegin,iend,iflag)
            call memptr(jeq,2*nvdd,'head',ibegin,iend,iflag)
            call memptr(ja0,maxja0,'head',ibegin,iend,iflag)
        endif
c
        maxn=2*ndf
        nvdd=0
        if(iprob.lt.0) nvdd=ip(71)
        if(ispd.eq.1) then
            ii=max0(13*ndf,7*ndf+2*maxn,12*ndf+nvdd)
        else
            ii=max0(13*ndf,12*ndf+2*maxn,12*ndf+nvdd)
        endif
        if(iabs(iprob).eq.4) then
            ii=max0(ii,6*ndf+10*nvdd)
        else if(iabs(iprob).eq.5) then
            ii=max0(ii,6*ndf+22*nvdd,20*ndf)
        else
            ii=max0(ii,6*ndf+4*nvdd)
        endif
        call memptr(izz,ii,'head',ibegin,iend,iflag)
        maxja=ip(87)
        if(maxja.lt.10*ndf) iflag=82
        maxa=ip(88)
        if(iabs(iprob).eq.5) then
            if(ispd.eq.1) then
                maxa=maxa/2
            else
                maxa=2*maxa/3
            endif
            igg=maxa+1
        else
            igg=1
        endif
        if(ispd.eq.1) then
            if(maxa.lt.8*ndf) iflag=82
        else
            if(maxa.lt.14*ndf) iflag=82
        endif
        if(iprob.lt.0) call exflag(iflag)
        if(iflag.ne.0) then
            ip(25)=82
            go to 20
        endif
c
c
        call cedge3(nvf,ntf,nbf,itnode,ibndry,w(ibedge),
     +      w(izz),iflag)
        if(iprob.lt.0) call exflag(iflag)
        if(iflag.ne.0) then
            ip(25)=iflag
            go to 20
        endif
c
c       continuation options
c
        if(iprob.eq.3) then
c
            call pltmgc(ip,rp,vx,vy,xm,ym,itnode,ibndry,w(ibedge),
     +          w(iuu),w(iu0),w(iudot),w(iu0dot),w(ievr),
     1          w(ievl),w(ium),w(iuc),w(ivx0),w(ivy0),ndof,w(itdof),
     2          w(ka),ja,a,w(ihh),a(igg),w(isu),w(ism),w(ibb),w(idd),
     3          w(ird),w(ipp),w(idl),w(ibdlwr),w(ibdupr),w(idu),
     4          w(idum),w(iduc),w(jtime),w(jhist),w(jpath),w(izz),
     5          a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
c       time dependent options
c
        else if(iprob.eq.6) then
            call pltmgp(ip,rp,vx,vy,xm,ym,itnode,ibndry,w(ibedge),
     +          w(iuu),w(iu0),w(iudot),w(iu0dot),w(ievr),
     1          w(ievl),w(ium),w(iuc),w(ivx0),w(ivy0),ndof,w(itdof),
     2          w(ka),ja,a,w(ihh),a(igg),w(isu),w(ism),w(ibb),w(idd),
     3          w(ird),w(ipp),w(idl),w(ibdlwr),w(ibdupr),w(idu),
     4          w(idum),w(iduc),w(jtime),w(jhist),w(jpath),w(izz),
     5          a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
c       obstacle problem
c
        else if(iprob.eq.1.or.iprob.eq.2) then
            call pltmgo(ip,rp,vx,vy,xm,ym,itnode,ibndry,w(ibedge),
     +          w(iuu),w(iu0),w(iudot),w(iu0dot),w(ievr),
     1          w(ievl),w(ium),w(iuc),w(ivx0),w(ivy0),ndof,w(itdof),
     2          w(ka),ja,a,w(ihh),a(igg),w(isu),w(ism),w(ibb),w(idd),
     3          w(ird),w(ipp),w(idl),w(ibdlwr),w(ibdupr),w(idu),
     4          w(idum),w(iduc),w(jtime),w(jhist),w(jpath),w(izz),
     5          a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
c       parameter identification problem
c
        else if(iprob.eq.4.or.iprob.eq.5) then
            call pltmgi(ip,rp,vx,vy,xm,ym,itnode,ibndry,w(ibedge),
     +          w(iuu),w(iu0),w(iudot),w(iu0dot),w(ievr),
     1          w(ievl),w(ium),w(iuc),w(ivx0),w(ivy0),ndof,w(itdof),
     2          w(ka),ja,a,w(ihh),a(igg),w(isu),w(ism),w(ibb),w(idd),
     3          w(ird),w(ipp),w(idl),w(ibdlwr),w(ibdupr),w(idu),
     4          w(idum),w(iduc),w(jtime),w(jhist),w(jpath),w(izz),
     5          a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
c       domain decomposition solve
c
        else if(iprob.lt.0) then
            call pltmgd(ip,rp,vx,vy,xm,ym,itnode,ibndry,w(ibedge),
     +          w(iuu),w(iu0),w(iudot),w(iu0dot),w(ievr),
     1          w(ievl),w(ium),w(iuc),w(ivx0),w(ivy0),ndof,w(itdof),
     2          w(ka),ja,a,w(ihh),a(igg),w(isu),w(ism),w(ibb),w(idd),
     3          w(ird),w(ipp),w(idl),w(ibdlwr),w(ibdupr),w(idu),
     4          w(idum),w(iduc),w(ipath),w(jeq),w(ja0),w(ia0),w(ih0),
     5          w(ig0),w(isu0),w(ism0),nvdd,w(iudd),w(jtime),w(jhist),
     6          w(jpath),w(izz),a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
        else
            ip(25)=7
        endif
c
        call timer(w(jtime),33)
c
   20   iflag=ip(25)
c
c       successful return
c
        if(iflag.eq.0) then
             if(ip(7).lt.0) then
                write(unit=sp(11),fmt='(a17,i2,a8,i2,a6,i8,a1)')
     +              'pltmg: ok (iprob=',ip(7),', itask=',ip(9),
     1             ', ndg=',ip(40),')'
             else
                write(unit=sp(11),fmt='(a17,i2,a8,i2,a6,i6,a1)')
     +              'pltmg: ok (iprob=',ip(7),', itask=',ip(9),
     1             ', ndf=',ip(5),')'
             endif
c
c       insufficient storage errors, wrong input data structure
c
        else if(iflag.ge.82.and.iflag.le.89) then
            write(unit=sp(11),fmt='(a11,i3,a22)')
     +          'pltmg error',iflag,': insufficient storage'
            if(nproc.gt.1) ip(24)=irgn
        else if(iflag.eq.25) then
            write(unit=sp(11),fmt='(a11,i3,a28)')
     +          'pltmg error',iflag,': wrong input data structure'
c
c       convergence errors
c
        else if(iflag.eq.1) then
            write(unit=sp(11),fmt='(a11,i2,a29)')
     +          'pltmg error',iflag,': zero pivot in factorization'
            if(nproc.gt.1) ip(24)=irgn
        else if(iflag.eq.2) then
            write(unit=sp(11),fmt='(a11,i2,a27)')
     +          'pltmg error',iflag,': newton line search failed'
            if(nproc.gt.1) ip(24)=irgn
        else if(iflag.eq.7) then
            write(unit=sp(11),fmt='(a11,i2,a22)')
     +          'pltmg error',iflag,': illegal problem type'
        else if(iflag.eq.9) then
            write(unit=sp(11),fmt='(a11,i2,a31)')
     +          'pltmg error',iflag,': continuation procedure failed'
        else if(iflag.eq.10) then
            write(unit=sp(11),fmt='(a11,i3,a29)')
     +          'pltmg error',iflag,': multigraph iteration failed'
            if(nproc.gt.1) ip(24)=irgn
        else if(iflag.eq.11) then
            if(ip(7).lt.0) then
                write(unit=sp(11),fmt='(a11,i3,a28)')
     +          'pltmg error',iflag,': newton/dd iteration failed'
            else
                write(unit=sp(11),fmt='(a11,i3,a25)')
     +          'pltmg error',iflag,': newton iteration failed'
            endif
            if(nproc.gt.1) ip(24)=irgn
        else if(iflag.eq.24) then
            write(unit=sp(11),fmt='(a11,i3,a8,i4)')
     +          'pltmg error',iflag,': region',ip(24)
        else if(iflag.eq.48) then
            write(unit=sp(11),fmt='(a11,i3,a12)')
     +          'pltmg error',iflag,': mpi is off'
        else if(iflag.eq.71.or.iflag.eq.72) then
            write(unit=sp(11),fmt='(a11,i3,a27)')
     +          'pltmg error',iflag,': dd solver not initialized'
            if(nproc.gt.1) ip(24)=irgn
        else
            write(unit=sp(11),fmt='(a11,i3,a15)')
     +          'pltmg error',iflag,': unknown error'
            if(nproc.gt.1) ip(24)=irgn
        endif
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine pltmgd(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +      u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1      ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2      ipath,jequv,ja0,a0,h0,g0,su0,sm0,nn,gf,
     3      time,hist,path,z,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100),itnode(5,*),ibndry(6,*),ibedge(2,*),
     1          ka(*),ja(*),ipath(6,*),jequv(*),ja0(*),
     2          itdof(ndof,*)
            double precision
     +          rp(100),vx(*),vy(*),xm(*),ym(*),u(*),u0(*),udot(*),
     1          u0dot(*),evr(*),evl(*),um(*),uc(*),vx0(*),vy0(*),
     2          a(*),h(*),g(*),su(*),sm(*),b(*),p(*),
     3          d(*),rd(*),dl(*),bdlwr(*),bdupr(*),du(*),dum(*),
     4          duc(*),time(3,*),hist(22,*),path(101,*),z(*),
     5          a0(*),h0(*),g0(*),su0(*),sm0(*),gf(nn,*),t(20)
            external a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy
            data ibit/0/
c
c       make sure the system is solved on each domain
c
        iprob=iabs(ip(7))
        ip(7)=iprob
        jflag=0
c
        if(iprob.eq.3)  then
            if(ip(9).lt.5.or.ip(9).gt.7) ip(9)=7
            call ctheta(ip,rp,jflag)
            if(jflag.ne.0) then
                ip(25)=9
                return
            endif
        else if(iprob.eq.1) then
            if(ip(9).ne.9) ip(9)=0
        else
            ip(9)=0
        endif
c
        call nwtt(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +      u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1      ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2      time,hist,z,-1,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
        irgn=ip(50)
        ip(7)=-ip(7)
        call exflag(ip(25))
        if(ip(25).ne.0) return
c
c       initialize global parameters on all processors
c
        nproc=ip(49)
        irgn=ip(50)
        eps=ceps(ibit)*1.0d2
        epsmg=dmax1(1.0d-4,eps)
        ndf=ip(5)
        if(iprob.eq.4) then
            t(1)=rp(21)
            call pl2ip(t,1)
            rl=t(1)/dfloat(nproc)
            rmu=rp(3)
            rllwr=rp(4)
            rlupr=rp(5)
            tol=dmax1(1.0d-2*rmu,eps)
c
            if(rlupr.ne.0.0d0) then
                rup=dabs(rlupr)*tol
            else
                rup=tol
            endif
            if(rllwr.ne.0.0d0) then
                rlw=dabs(rllwr)*tol
            else
                rlw=tol
            endif
            if(rllwr+rlw.le.rlupr-rup) then
                rl=dmax1(rl,rllwr+rlw)
                rl=dmin1(rl,rlupr-rup)
            else
                rr=tol*(rlupr-rllwr)
                rl=dmax1(rl,rllwr+rr)
                rl=dmin1(rl,rlupr-rr)
            endif
c
            rp(21)=rl
        else if(iprob.eq.3) then
            do k=1,ndf
                u0(k)=u(k)
                u0dot(k)=udot(k)
            enddo
            do k=1,7
               t(k)=rp(20+k)
            enddo
            t(8)=rp(68)
            call pl2ip(t,8)
            do k=1,7
               rp(20+k)=t(k)/dfloat(nproc)
            enddo
            rp(68)=t(8)/dfloat(nproc)
            do k=1,5
                rp(30+k)=rp(20+k)
            enddo
        endif
c
c       initialize for domain decomposition solve
c
        ntf=ip(1)
        ndf=ip(5)
        ising=ip(12)
        itask=ip(9)
        newndf=ip(30)
        ndd=ip(33)
        ndi=ip(36)
        iord=ip(26)
c
        nvdd=ip(71)
        lipath=ip(72)
        if(nvdd.le.0) then
            ip(25)=71
        else if(lipath.le.0) then
            ip(25)=72
        else if(ipath(2,nproc+2).lt.ipath(1,nproc+2)) then
            ip(25)=72
        else
            ip(25)=0
        endif
        call exflag(ip(25))
        if(ip(25).ne.0) return
c
c       initialize for domain decomposition
c
        rp(52)=1.0d0
        rp(56)=1.0d0
        rp(57)=1.0d0
        rp(54)=0.0d0
c*      if(iprob.eq.2) rp(63)=rp(3)
c*      if(iprob.eq.4) rp(63)=rp(3)
c*      if(iprob.eq.5) rp(63)=rp(3)
c
        mxdamp=10
        ievals=1
        itype=0
        iconv=0
        ispd=ip(8)
        jspd=1
        if(ispd.ne.1) jspd=-1
        mxcg=ip(10)
        mxnwtt=ip(11)
        jnwtt=mxnwtt
        if(iprob.eq.3) jnwtt=mxnwtt+1
c
        ir=1
        img=ir+ndf
        igm=img+ndf
        iadu=igm+ndf
        ihdu=iadu+ndf
        iadm=ihdu+ndf
        ismdm=iadm+ndf
        ismdc=ismdm+ndf
        isudu=ismdc+ndf
        isudc=isudu+ndf
        igdc=isudc+ndf
        isv=igdc+ndf
        iin=isv+ndf
        imk=iadu
        id1=ihdu
        id2=iadm
        mm1=1
        mm2=mm1+3*ndf
        mm3=mm2+3*ndf
        mm4=mm3+3*ndf
        mm5=mm4+3*ndf
        mm6=mm5+3*ndf
c
        maxja0=9*iord*nvdd/2
        call cequv2(irgn,nproc,iord,newndf,ndi,ipath,jequv,ja0)
        call setgr2(irgn,nproc,ntf,ndd,newndf,ndi,itnode,
     +      ndof,itdof,ipath,jequv,ja0,z(ir),maxja0,kflag)
        if(kflag.ne.0) stop 9011
c
        iqptr=ja(ndf+1)
c
c       linear system
c
        call timer(time,33)
        call rgnsys(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +      u,u0,udot,um,uc,z(id1),z(id2),vx0,vy0,ndof,itdof,
     1      ja,ja(iqptr),a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,z(imk),
     2      z(imk),jequv,ipath,ja0,a0,h0,g0,su0,sm0,nn,gf,
     3      z(igm),z(iin),1,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
        call timer(time,32)
c
c       solve equations
c
        call hist3(hist(1,11),-1,1.0d0,1.0d0)
        do itnum=1,jnwtt
c
            if(itnum.gt.0) then
                call timer(time,33)
                call mgilu(ja,a,ka,z(img))
                call timer(time,22)
            endif
c
c       multi-level solution of newton equations
c
            call timer(time,33)
            if(iprob.eq.3) then
                call blk3(ip,rp,vx,vy,ndof,itdof,itnode,du,dum,
     +              ka,ja,a,b,rd,p,udot,u0dot,epsmg,z(ir),
     1              z(img),hist,jflag,1)
                call timer(time,26)
                if(iconv.eq.1) go to 170
                if(itnum.gt.mxnwtt) go to 100
            else if(iprob.eq.4) then
                call blk4(ip,rp,du,dum,ka,ja,a,h,b,p,dl,
     +              rd,udot,epsmg,z(ir),z(img),hist,jflag,1)
                call timer(time,27)
            else if(iprob.eq.5) then
                call mgilu(ja,g,ka(101),z(img))
                call timer(time,22)
                call blk5(ip,epsmg,ka,ja,a,h,g,su,sm,
     +              du,dum,duc,p,b,dl,z(mm1),z(mm2),z(mm3),
     1              z(mm4),z(mm5),hist(1,7),z(mm6),reler5,jflag)
c               call blk5x(ip,du,dum,duc,ka,ja,a,h,g,su,sm,
c    +              b,p,dl,epsmg,z(ir),z(img),hist,jflag)
                call timer(time,28)
            else
                do i=1,ndf
                    z(ir+i-1)=b(i)
                enddo
                call mg(ispd,mxcg,epsmg,ja,a,du,z(ir),
     +              ka,ising,reler1,jflag,z(img),hist(1,7))
                if(iprob.eq.1.and.itask.eq.9) then
                    do i=1,ndf
                        z(ir+i-1)=p(i)
                    enddo
                    call mg(jspd,mxcg,epsmg,ja,a,dum,z(ir),
     +                  ka,ising,reler2,jflag,z(img),hist(1,8))
                endif
                call timer(time,21)
            endif
c
c    line search loop
c
            isw=0
            call timer(time,33)
            call tpickd(ip,rp,vx,vy,itnode,u,um,uc,z(isv),
     +          z(ir),z(img),ja,a,h,g,su,sm,b,d,p,dl,bdlwr,bdupr,
     1          du,dum,duc,ipath,jequv,ja0,a0,h0,g0,su0,sm0,nn,gf,
     2          z(iadu),z(iadm),z(ihdu),z(igdc),z(isudu),z(isudc),
     3          z(ismdm),z(ismdc),z(igm),z(iin),isw,itnum,ndof,itdof)
            call timer(time,31)
c
c       initializization for itnum=1
c
            dnew=rp(58)
            if(dnew.gt.0.0d0) then
                call hist3(hist(1,11),itnum,rp(56),rp(54))
                iconv=icvtst(itnum,-iprob,itask,itype,rp)
c****           iconv=jcvtst(itnum,-iprob,itask,itype,rp)
                if(iconv.eq.1) go to 170
                ip(25)=2
                if(jflag.ne.0) ip(25)=11
                go to 130
            endif
            iter=0
   70       iter=iter+1
c
c       linear system
c
            call timer(time,33)
            call rgnsys(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +          u,u0,udot,um,uc,z(id1),z(id2),vx0,vy0,ndof,itdof,
     1          ja,ja(iqptr),a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,z(imk),
     2          z(imk),jequv,ipath,ja0,a0,h0,g0,su0,sm0,nn,gf,
     3          z(igm),z(iin),0,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
            ievals=ievals+1
            call timer(time,32)
c
            call tpickd(ip,rp,vx,vy,itnode,u,um,uc,z(isv),
     +          z(ir),z(img),ja,a,h,g,su,sm,b,d,p,dl,bdlwr,bdupr,
     1          du,dum,duc,ipath,jequv,ja0,a0,h0,g0,su0,sm0,nn,gf,
     2          z(iadu),z(iadm),z(ihdu),z(igdc),z(isudu),z(isudc),
     3          z(ismdm),z(ismdc),z(igm),z(iin),isw,itnum,ndof,itdof)
            call timer(time,31)
            if(isw.ge.0) then
                if(iter.lt.mxdamp) go to 70
                ip(25)=2
                return
            endif
c
c       convergence check
c
            call hist3(hist(1,11),itnum,rp(56),rp(54))
            iconv=icvtst(itnum,-iprob,itask,itype,rp)
c****       iconv=jcvtst(itnum,-iprob,itask,itype,rp)
            if(iprob.ne.3.and.iconv.eq.1) go to 170
        enddo
  100   itnum=jnwtt
        if(iconv.eq.-1) go to 170
        ip(25)=11
c
c       newton iteration failed to converge...reset u, udot, and rp
c
  130   ip(79)=ievals
        ip(80)=itnum
        return
c
c       newton iteration was successful
c
  170   ip(25)=0
        ip(79)=ievals
        ip(80)=itnum
        if(iprob.eq.3) then
            ip(80)=itnum-1
            call updpth(path,-1,5,rp)
        else if(iprob.ne.1) then
            call updip(path,-1,4,rp,ip)
        endif
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine pltmgp(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +      u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1      ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2      time,hist,path,z,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100),itnode(5,*),ibndry(6,*),ibedge(2,*),
     1          ka(*),ja(*),itdof(ndof,*)
            double precision
     +          rp(100),vx(*),vy(*),xm(*),ym(*),u(*),u0(*),udot(*),
     1          u0dot(*),evr(*),evl(*),um(*),uc(*),vx0(*),vy0(*),
     2          a(*),h(*),g(*),su(*),sm(*),b(*),p(*),
     3          d(*),rd(*),dl(*),bdlwr(*),bdupr(*),du(*),dum(*),
     4          duc(*),time(3,*),hist(22,*),path(101,*),z(*)
            character*80
     +          iostr,msg
            external a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy
            save msg
            data msg/'pltmg:      tcur       deltat        utnorm'/
c
c       main routine for parabolic problems
c
        itask=ip(9)
        ntf=ip(1)
        nvf=ip(2)
        ndf=ip(5)
        iord=ip(26)
c
        call filutl(msg,0)
c
        call uinit(ip,rp,itnode,ibndry,ibedge,vx,vy,xm,ym,
     +      u,um,uc,z,ndof,itdof,z(ndf+1),gdxy)
c
        if(itask.eq.10) then
            tstart=rp(42)
            tend=rp(43)
            if(tstart.ge.tend) return
            mxstep=max0(1,ip(15))
            mxfail=5
            rp(46)=tstart
            rp(49)=tend-tstart
            rp(48)=rp(49)/dfloat(mxstep)
            tnew=rp(46)
            ifirst=1
c
c       compute time step
c
   60       call dtpick(ntf,ndf,itnode,vx,vy,u,u0,rp,z,itflag,ifirst,
     +          z(ndf+1),ndof,itdof)
c
c       update solution
c
            if(itflag.ne.-1.and.ifirst.ne.-1) then
                rp(46)=tnew
                do i=1,ndf
                    u0(i)=u(i)
                enddo
                do i=1,nvf
                    vx0(i)=vx(i)
                    vy0(i)=vy(i)
                enddo
                idsp=0
            endif
            if(ifirst.eq.-1) then
                rp(46)=tnew
                rp(42)=tnew
                rp(43)=tnew
            endif
c
c       save time history
c
            if(ifirst.eq.1) then
                if(itflag.le.-3)  then
                    call updtm(path,1,itflag,rp)
                else
                    call updtm(path,0,itflag,rp)
                endif
            else
                if(itflag.eq.-1)  then
                    call updtm(path,0,itflag,rp)
                else
                    call updtm(path,-1,itflag,rp)
                endif
            endif
            write(unit=iostr,fmt='(2i3,3(1x,e12.5))')
     +          ip(25),itflag,rp(46),rp(47),rp(50)
            call filutl(iostr,0)
            if(ifirst.eq.-1) return
            ifirst=0
c
c       solve equations
c
  220       idsp=idsp+1
            tcur=rp(46)
            deltat=dmax1(rp(47),rp(48))
            rp(21)=tcur+deltat
            if(deltat.gt.0) then
                rp(45)=1.0d0/deltat
            else
                rp(45)=0.0d0
            endif
            call nwtt(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +          u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1          ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2          time,hist,z,0,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
            if(ip(25).ne.0) then
                if(idsp.lt.mxfail) then
                    rp(47)=rp(47)/2.0d0
                    go to 220
                else
                    do i=1,ndf
                        u(i)=u0(i)
                    enddo
                    write(unit=iostr,fmt='(2i3,3(1x,e12.5))')
     +                  ip(25),itflag,rp(46),rp(47),rp(50)
                    call filutl(iostr,0)
                    return
                endif
            else
                tnew=rp(46)+rp(47)
                if(itflag.eq.2.and.idsp.eq.1) ifirst=-1
                if(itflag.eq.-4.and.idsp.eq.1) ifirst=-1
                go to 60
            endif
        else
            tcur=rp(46)
            deltat=dmax1(rp(47),rp(48))
            rp(21)=tcur
            if(deltat.gt.0) then
                rp(45)=1.0d0/deltat
            else
                rp(45)=0.0d0
            endif
            call nwtt(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +          u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1          ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2          time,hist,z,0,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
            itflag=3
            write(unit=iostr,fmt='(2i3,3(1x,e12.5))')
     +          ip(25),itflag,rp(46),rp(47),rp(50)
            call filutl(iostr,0)
            call updtm(path,0,3,rp)
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine pltmgc(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +      u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1      ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2      time,hist,path,z,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100),itnode(5,*),ibndry(6,*),ibedge(2,*),
     1          ka(*),ja(*),itdof(ndof,*)
            double precision
     +          rp(100),vx(*),vy(*),xm(*),ym(*),u(*),u0(*),udot(*),
     1          u0dot(*),evr(*),evl(*),um(*),uc(*),vx0(*),vy0(*),
     2          a(*),h(*),g(*),su(*),sm(*),b(*),p(*),
     3          d(*),rd(*),dl(*),bdlwr(*),bdupr(*),du(*),dum(*),
     4          duc(*),time(3,*),hist(22,*),path(101,*),z(*)
            character*80
     +          iostr,msg(7)
            external a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy
            save msg
            data msg(1)/'pltmg:     lambda        rho       lambda dot
     +  rho dot     eigenvalue'/
            data msg(2)/'pltmg: find limit / bifurcation point'/
            data msg(3)/'pltmg: probable limit point'/
            data msg(4)/'pltmg: probable regular point'/
            data msg(5)/'pltmg: probable bifurcation point'/
            data ibit/0/
c
c       continuation
c
        itask=ip(9)
        ispd=ip(8)
        ntf=ip(1)
        nbf=ip(4)
        ndf=ip(5)
        eps=1.0d2*ceps(ibit)
        iord=ip(26)
c
c
        call filutl(msg(1),0)
c
        istep=0
        idsp=0
        mxbis=10
        mxfail=10
        mxstep=10
c
        igm=1
        ipp=igm+ndf
        iz1=ipp+ndf
        iz2=iz1+ndf
        iz3=iz2+ndf
        iz4=iz3+ndf
c
c       restore solution
c
        call uinit(ip,rp,itnode,ibndry,ibedge,vx,vy,xm,ym,
     +      u,um,uc,z(igm),ndof,itdof,z(ipp),gdxy)
        do i=1,ndf
            u(i)=u0(i)
            udot(i)=u0dot(i)
        enddo
        do i=1,5
            rp(20+i)=rp(30+i)
        enddo
        rltrgt=rp(1)
        rtrgt=rp(2)
        rp(26)=rp(31)
        rp(27)=rp(32)
c
c       change itask if things look inconsistant
c
        dd=dabs(rltrgt-rp(21))+dabs(rtrgt-rp(22))
        if(dd.eq.0.0d0.and.itask.le.1) itask=7
        if(dd.ne.0.0d0.and.itask.ge.5) itask=0
        ip(9)=itask
c
c       switch branches at bifurcation point
c
        if(itask.eq.2) then
            call timer(time,33)
            call swbrch(ndf,ntf,nbf,itnode,ibndry,ndof,itdof,vx,vy,
     +          xm,ym,evl,evr,udot,u,u0dot,z(ipp),z(iz1),z(iz2),
     1          z(iz3),z(igm),z(iz4),rp,ibedge,ispd,iord,
     2          a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy,0)
            call timer(time,29)
            call updpth(path,0,6,rp)
            do i=1,ndf
                u0dot(i)=udot(i)
            enddo
            rp(33)=rp(23)
            rp(34)=rp(24)
            ip(9)=0
            ip(80)=0
            write(unit=iostr,fmt='(2i3,5(1x,e12.5))')
     +          ip(25),ip(80),(rp(k),k=21,25)
            call filutl(iostr,0)
            return
        endif
c
c       switch functional and/or parameters
c
        if(itask.ge.3) then
            call ctheta(ip,rp,iflag)
            if(iflag.ne.0) then
                ip(25)=9
                return
            endif
            call nwtt(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +          u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1          ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2          time,hist,z,0,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
            write(unit=iostr,fmt='(2i3,5(1x,e12.5))')
     +          ip(25),ip(80),(rp(k),k=21,25)
            call filutl(iostr,0)
            if(ip(25).ne.0) then
                ip(9)=itask
                rp(1)=rp(21)
                rp(2)=rp(22)
                return
            else
                if(itask.le.4) then
                    call updpth(path,1,1,rp)
                    ip(9)=0
                else
                    call updpth(path,-1,3,rp)
                endif
                go to 40
            endif
        endif
c
c       get set for an arc length continuation step
c
   10   idsp=0
        istep=istep+1
        if(istep.gt.mxstep) then
            ip(25)=9
            ip(9)=itask
            rp(1)=rp(21)
            rp(2)=rp(22)
            return
        endif
c
c       step picker
c
        call timer(time,33)
        call predct(ip,itnode,ibndry,vx,vy,xm,ym,z(ipp),z(iz1),
     +      z(igm),u0,u0dot,rp,ibedge,idsp,mxfail,ndof,itdof,
     1      z(iz2),z(iz3),z(iz3),a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
        call timer(time,25)
        if(idsp.gt.mxfail) then
            ip(25)=9
            ip(9)=itask
            rp(1)=rp(21)
            rp(2)=rp(22)
            return
        endif
c
c       solve nonlinear equations
c
        call nwtt(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +      u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1      ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2      time,hist,z,0,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
        write(unit=iostr,fmt='(2i3,5(1x,e12.5))')
     +      ip(25),ip(80),(rp(k),k=21,25)
        call filutl(iostr,0)
        if(ip(25).ne.0) then
            ip(9)=itask
            rp(1)=rp(21)
            rp(2)=rp(22)
            return
        endif
        sval=rp(25)
        sval0=rp(35)
        if(istep.eq.1) then
            call updpth(path,-1,4,rp)
        else
            call updpth(path,0,4,rp)
        endif
        if(sval0*sval.ge.0.0d0.or.itask.eq.0) go to 40
c
c       change in sign in determinent
c
        call filutl(msg(2),0)
c
c       information for testing type of singular point
c
        rqmx=dmax1(dabs(sval),dabs(sval0))
        rlsign=rp(23)*rp(33)
        idsp=0
        isw=0
        call hist3(hist(1,15),-2,sval,sval0)
c
        do istep=1,mxbis
c
c       bisection/secant step
c
            call bisect(rp,isw,rqup,rqlow)
            call hist3(hist(1,15),istep,rqup,rqlow)
            if(isw.eq.-1) go to 30
            sigma=rp(71)
   20       call nwtt(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +          u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1          ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2          time,hist,z,1,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
            write(unit=iostr,fmt='(2i3,5(1x,e12.5))')
     +          ip(25),ip(80),(rp(k),k=21,25)
            call filutl(iostr,0)
            if(ip(25).ne.0) then
                if(dabs(sigma).eq.dabs(rp(71))) then
                    rp(71)=sigma*(1.0d0-eps)
                    go to 20
                else if(dabs(sigma).lt.dabs(rp(71))) then
                    rp(71)=sigma*(1.0d0+eps)
                    go to 20
                else
                    rp(1)=rp(31)
                    rp(2)=rp(32)
                    return
                endif
            endif
        enddo
   30   ip(9)=0
c
c       fixup tangent for the case of bifurcation
c
        dnorm=rl2nrm(ndf,udot)*rl2nrm(ndf,evr)
        if(dnorm.gt.0.0d0) dnorm=1.0d0/dnorm
        udr=rl2ip(ndf,evr,udot)*dnorm
        if(dabs(udr).gt.1.0d-1.and.rlsign.lt.0.0d0) then
            call filutl(msg(3),0)
            call updpth(path,0,2,rp)
        else if(dabs(rp(25)).gt.rqmx*1.0d-2) then
            call filutl(msg(4),0)
            call updpth(path,0,4,rp)
        else
            call filutl(msg(5),0)
            call updpth(path,0,6,rp)
            call timer(time,33)
            call swbrch(ndf,ntf,nbf,itnode,ibndry,ndof,itdof,vx,vy,
     +          xm,ym,evl,evr,udot,u,u0dot,z(ipp),z(iz1),z(iz2),
     1          z(iz3),z(igm),z(iz4),rp,ibedge,ispd,iord,
     2          a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy,1)
            call timer(time,29)
            ip(80)=0
            write(unit=iostr,fmt='(2i3,5(1x,e12.5))')
     +          ip(25),ip(80),(rp(k),k=21,25)
            call filutl(iostr,0)
        endif
c
c       successful continuation
c
   40   do i=1,5
            rp(30+i)=rp(20+i)
        enddo
        do i=1,ndf
            u0(i)=u(i)
            u0dot(i)=udot(i)
        enddo
        if(idsp.ne.0) go to 10
        rp(1)=rp(31)
        rp(2)=rp(32)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine pltmgo(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +      u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1      ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2      time,hist,path,z,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100),itnode(5,*),ibndry(6,*),ibedge(2,*),
     1          ka(*),ja(*),itdof(ndof,*)
            double precision
     +          rp(100),vx(*),vy(*),xm(*),ym(*),u(*),u0(*),udot(*),
     1          u0dot(*),evr(*),evl(*),um(*),uc(*),vx0(*),vy0(*),
     2          a(*),h(*),g(*),su(*),sm(*),b(*),p(*),
     3          d(*),rd(*),dl(*),bdlwr(*),bdupr(*),du(*),dum(*),
     4          duc(*),time(3,*),hist(22,*),path(101,*),z(*)
            character*80
     +          iostr
            external a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy
            save isw
            data isw/1/
c
c       solve equations
c
        rp(21)=rp(1)
        rp(22)=rp(2)
        iprob=ip(7)
        if(iprob.eq.2) then
            rp(63)=rp(3)
        else
            if(ip(9).ne.9) ip(9)=0
        endif
        call nwtt(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +      u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1      ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2      time,hist,z,0,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
        if(iprob.eq.2.and.ip(25).eq.0) then
            if(isw.eq.1) then
                call updip(path,1,1,rp,ip)
                isw=0
            else
                call updip(path,-1,2,rp,ip)
            endif
            write(unit=iostr,fmt='(a11,e12.5,3x,a3,e12.5)')
     +          'pltmg: rho=',rp(22),'mu=',rp(63)
            call filutl(iostr,0)
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine pltmgi(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +      u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1      ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2      time,hist,path,z,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100),itnode(5,*),ibndry(6,*),ibedge(2,*),
     1          ka(*),ja(*),itdof(ndof,*)
            double precision
     +          rp(100),vx(*),vy(*),xm(*),ym(*),u(*),u0(*),udot(*),
     1          u0dot(*),evr(*),evl(*),um(*),uc(*),vx0(*),vy0(*),
     2          a(*),h(*),g(*),su(*),sm(*),b(*),p(*),
     3          d(*),rd(*),dl(*),bdlwr(*),bdupr(*),du(*),dum(*),
     4          duc(*),time(3,*),hist(22,*),path(101,*),z(*)
            character*80
     +          iostr
            external a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy
            save isw
            data isw/1/
            data ibit/0/
c
c       solve equations
c
        iprob=ip(7)
        itask=ip(9)
        if(iprob.eq.4) then
            if(itask.eq.8) then
                rp(21)=rp(1)
cc              rp(21)=(rp(4)+rp(5))/2.0e0
                ip(9)=0
            endif
            rmu=rp(3)
            rllwr=rp(4)
            rlupr=rp(5)
            eps=100.0d0*ceps(ibit)
            tol=dmax1(1.0d-2*rmu,eps)
            rl=rp(21)
c
            if(rlupr.ne.0.0d0) then
                rup=dabs(rlupr)*tol
            else
                rup=tol
            endif
            if(rllwr.ne.0.0d0) then
                rlw=dabs(rllwr)*tol
            else
                rlw=tol
            endif
            if(rllwr+rlw.le.rlupr-rup) then
                rl=dmax1(rl,rllwr+rlw)
                rl=dmin1(rl,rlupr-rup)
            else
                rr=tol*(rlupr-rllwr)
                rl=dmax1(rl,rllwr+rr)
                rl=dmin1(rl,rlupr-rr)
            endif
c
            rp(21)=rl
        endif
        rp(63)=rp(3)
        call nwtt(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +      u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1      ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2      time,hist,z,0,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
        if(ip(25).eq.0) then
            if(isw.eq.1) then
                call updip(path,1,1,rp,ip)
                isw=0
            else if(iprob.eq.4.and.itask.eq.8) then
                call updip(path,-1,3,rp,ip)
            else
                call updip(path,-1,2,rp,ip)
            endif
        endif
        if(iprob.eq.4) then
            write(unit=iostr,fmt='(a11,e12.5,3x,a7,e12.5,3x,a3,e12.5)')
     +          'pltmg: rho=',rp(22),'lambda=',rp(21),'mu=',rp(63)
        else
            write(unit=iostr,fmt='(a11,e12.5,3x,a3,e12.5)')
     +          'pltmg: rho=',rp(22),'mu=',rp(63)
        endif
        call filutl(iostr,0)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine nwtt(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +      u,u0,udot,u0dot,evr,evl,um,uc,vx0,vy0,ndof,itdof,
     1      ka,ja,a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,du,dum,duc,
     2      time,hist,z,itype,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100),itnode(5,*),ibndry(6,*),ibedge(2,*),
     1          ka(*),ja(*),itdof(ndof,*)
            double precision
     +          rp(100),vx(*),vy(*),xm(*),ym(*),u(*),u0(*),udot(*),
     1          u0dot(*),evr(*),evl(*),um(*),uc(*),vx0(*),vy0(*),
     2          a(*),h(*),g(*),su(*),sm(*),b(*),d(*),
     3          rd(*),p(*),dl(*),bdlwr(*),bdupr(*),du(*),dum(*),
     4          duc(*),time(3,*),hist(22,*),z(*),rpsv(100)
            external a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy
            data ibit/0/
c
c       approximate newton method
c
        ntf=ip(1)
        ndf=ip(5)
        itask=ip(9)
        iprob=ip(7)
        ising=ip(12)
        iord=ip(26)
        eps=1.0d2*ceps(ibit)
c
        rp(52)=1.0d0
        rp(56)=1.0d0
        rp(57)=1.0d0
        if(itype.eq.0) then
            epsmg=dmax1(1.0d-3,eps)
        else
            epsmg=dmax1(1.0d-4,eps)
        endif
        epsmg0=epsmg
        mxdamp=20
        iconv=0
        jflag=0
c
        ncfact=4
        maxlvl=20
        ispd=ip(8)
        jspd=1
        if(ispd.ne.1) jspd=-1
        maxja=ip(87)
        maxa=ip(88)
        if(iprob.eq.5) then
            if(ispd.eq.1) then
                maxa=maxa/2
            else
                maxa=2*maxa/3
            endif
            maxg=ip(88)-maxa
        endif
        lenz=30*ndf
        mxcg=ip(10)
        mxnwtt=ip(11)
        jnwtt=mxnwtt
        if(iprob.eq.3) jnwtt=mxnwtt+1
        dtol=rp(6)
c
c       get pointers  (len = 13 ispd=0, 16/17 ispd=0,iprob eq/ne 3
c
        isav=1
        m1=isav+ndf
        m2=m1+ndf
        m3=m2+ndf
        m4=m3+ndf
        m5=m4+ndf
        m6=m5+ndf
        img=m1
        if(ispd.eq.1) then
            ns1=m1
            ns2=m2
            ns3=m3
            kmg=m4
        else
            ns1=m4
            ns2=m5
            ns3=m6
            kmg=m6+ndf
        endif
        mm1=1
        mm2=mm1+3*ndf
        mm3=mm2+3*ndf
        mm4=mm3+3*ndf
        mm5=mm4+3*ndf
        mm6=mm5+3*ndf
c
c       save rp
c
        do i=1,100
            rpsv(i)=rp(i)
        enddo
c
        maxlnk=ndf*4
        if(iord.eq.2) maxlnk=ndf*25/4
        if(iord.eq.3) maxlnk=ndf*81/9
        call setgr1(ntf,ndf,ndof,itdof,ja,a,maxlnk,kflag)
        call mdinit(ndf,ispd,ja,ka,lenz,z,kflag)
        if(kflag.ne.0) then
            ip(25)=kflag
            return
        endif
c
        iqptr=ja(ndf+1)
        call uinit(ip,rp,itnode,ibndry,ibedge,vx,vy,xm,ym,
     +      u,um,uc,z(m1),ndof,itdof,z(m2),gdxy)
        if(iprob.eq.3) call evinit(ip,evl,evr,ndof,itdof,ibndry,
     +      ibedge,ja(iqptr),z(m1))
        if(iprob.eq.2) call bdinit(ip,rp,u,vx,vy,xm,ym,ndof,itdof,
     +      itnode,ibndry,ibedge,bdlwr,bdupr,z(m1),z(m2),gdxy)
        if(iprob.eq.5) call bdinit(ip,rp,uc,vx,vy,xm,ym,ndof,itdof,
     +      itnode,ibndry,ibedge,bdlwr,bdupr,z(m1),z(m2),gdxy)
c
        if(iprob.eq.3.and.itask.le.1) then
            seqdot=rp(74)
            sigma=rp(71)
            rl0dot=rp(33)
            if(seqdot.ne.0.0d0) then
                ss=dsqrt(eps)*rl0dot*sigma/seqdot
            else
                ss=dsqrt(eps)*rl0dot
            endif
            do j=1,ndf
                u(j)=u(j)+ss*u0dot(j)
            enddo
        endif
c
c       first matrix and right hand side
c
        call timer(time,33)
        call linsys(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +      u,u0,udot,um,uc,z(m3),z(m4),vx0,vy0,ndof,itdof,
     1      ja,ja(iqptr),a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,z(m5),
     2      z(m6),1,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
        call timer(time,24)
c
c       print out matrix and return (for mutigraph testing)
c
c       call mkmtx(ndf,ispd,ja,a,b,ja(iqptr),z(m1))
c       ip(25)=82
c       if(ip(25).ne.0) return
c
        ievals=1
c
c       compute ordering symbolic factorization
c
        call timer(time,33)
        call mgilu0(maxja,ja,maxa,a,ncfact,maxlvl,ka,dtol,z,kflag)
        call timer(time,20)
        if(kflag.ne.0) then
            ip(25)=kflag
            return
        endif
        ip(73)=ka(17)
        ip(74)=ka(18)
        ip(75)=ka(2)
c
c       ka array for matrix g
c
        if(iprob.eq.5) then
            do i=1,100
                ka(i+100)=ka(i)
            enddo
            ka(106)=ka(6)+ka(10)
            ka(107)=ka(3)
            ka(114)=1
            call timer(time,33)
            call mgilu0(maxja,ja,maxg,g,ncfact,maxlvl,ka(101),
     +          dtol,z,kflag)
            call timer(time,20)
            if(kflag.ne.0) then
                ip(25)=kflag
                return
            endif
        endif
c
c       the main loop
c
        call hist3(hist(1,11),0,1.0d0,1.0d0)
        do itnum=1,jnwtt
c
c       compute approximate factorization
c
            if(itnum.gt.1) then
                call timer(time,33)
                call mgilu(ja,a,ka,z(img))
                call timer(time,22)
                if(itype.eq.0) then
                    epsmg=dmax1(epsmg0,rp(57))
                    epsmg=dmin1(1.0d-1,rp(57))
                endif
            endif
c
c       compute singular vectors
c
            if(iprob.eq.3) then
                call timer(time,33)
                call cev(ip,rp,ja,a,ka,evl,evr,z(m1),z(ns1),
     +              z(m2),z(ns2),z(m3),z(ns3),z(kmg),hist)
                call timer(time,23)
            endif
c
c       multi-level solution of newton equations
c
            call timer(time,33)
            if(iprob.eq.3) then
                call blk3(ip,rp,vx,vy,ndof,itdof,itnode,du,dum,
     +              ka,ja,a,b,rd,p,udot,u0dot,epsmg,z(isav),
     1              z(img),hist,jflag,0)
                call timer(time,26)
                if(iconv.eq.1) go to 170
                if(itnum.gt.mxnwtt) go to 100
            else if(iprob.eq.4) then
                call blk4(ip,rp,du,dum,ka,ja,a,h,b,p,dl,
     +              rd,udot,epsmg,z(isav),z(img),hist,jflag,0)
                call timer(time,27)
            else if(iprob.eq.5) then
                if(itnum.gt.1) then
                    call mgilu(ja,g,ka(101),z(img))
                    call timer(time,22)
                endif
                call blk5(ip,epsmg,ka,ja,a,h,g,su,sm,
     +              du,dum,duc,p,b,dl,z(mm1),z(mm2),z(mm3),
     1              z(mm4),z(mm5),hist(1,7),z(mm6),reler5,jflag)
c               call blk5x(ip,du,dum,duc,ka,ja,a,h,g,su,sm,
c    +              b,p,dl,epsmg,z(isav),z(img),hist,jflag)
                call timer(time,28)
            else
                do j=1,ndf
                    z(isav+j-1)=b(j)
                enddo
                call mg(ispd,mxcg,epsmg,ja,a,du,z(isav),
     +              ka,ising,reler1,jflag,z(img),hist(1,7))
                if(iprob.eq.1.and.itask.eq.9) then
                    do j=1,ndf
                        z(isav+j-1)=p(j)
                    enddo
                    call mg(jspd,mxcg,epsmg,ja,a,dum,z(isav),
     +                  ka,ising,reler2,jflag,z(img),hist(1,8))
                endif
                call timer(time,21)
            endif
c
c       line search and sufficient decrease loop
c
            isw=0
            call timer(time,33)
            call tpick(ip,rp,vx,vy,itnode,ndof,itdof,u,um,uc,
     +          z(isav),z(m1),z(m2),ja,a,h,g,su,sm,b,d,p,dl,bdlwr,
     1          bdupr,du,dum,duc,z(m3),z(m4),z(m5),z(m6),isw,itnum)
            call timer(time,30)
            dnew=rp(58)
cc          write(6,*) itnum,dnew,rp(52)
            if(dnew.gt.0.0d0) then
                call hist3(hist(1,11),itnum,rp(56),rp(54))
                iconv=icvtst(itnum,iprob,itask,itype,rp)
                if(iconv.eq.1) go to 170
                ip(25)=2
                if(jflag.ne.0) ip(25)=11
                go to 130
            endif
            iter=0
   70       iter=iter+1
c
            call timer(time,33)
            call linsys(ip,rp,vx,vy,xm,ym,itnode,ibndry,ibedge,
     +          u,u0,udot,um,uc,z(m3),z(m4),vx0,vy0,ndof,itdof,
     1          ja,ja(iqptr),a,h,g,su,sm,b,d,rd,p,dl,bdlwr,bdupr,z(m5),
     2          z(m6),0,a1xy,a2xy,fxy,gnxy,gdxy,p1xy,p2xy)
            ievals=ievals+1
            call timer(time,24)
            call tpick(ip,rp,vx,vy,itnode,ndof,itdof,u,um,uc,
     +          z(isav),z(m1),z(m2),ja,a,h,g,su,sm,b,d,p,dl,bdlwr,
     1          bdupr,du,dum,duc,z(m3),z(m4),z(m5),z(m6),isw,itnum)
            call timer(time,30)
c
c       test for sufficient decrease
c
            if(isw.ge.0) then
                if(iter.lt.mxdamp) go to 70
                ip(25)=2
                if(jflag.ne.0) ip(25)=11
                go to 130
            endif
c
c       convergence test
c
            call hist3(hist(1,11),itnum,rp(56),rp(54))
            iconv=icvtst(itnum,iprob,itask,itype,rp)
            if(iprob.ne.3.and.iconv.eq.1) go to 170
c
        enddo
  100   itnum=jnwtt
        if(iconv.eq.-1) go to 170
        ip(25)=11
c
c       newton iteration failed to converge...reset u, udot, and rp
c
  130   ip(79)=ievals
        ip(80)=itnum
        do i=1,100
            rp(i)=rpsv(i)
        enddo
        if(iprob.eq.3) then
            do j=1,ndf
                udot(j)=u0dot(j)
                u(j)=u0(j)
            enddo
        else if(iprob.eq.6) then
            do j=1,ndf
                u(j)=u0(j)
            enddo
        endif
        return
c
c       newton iteration was successful
c
  170   ip(25)=0
        ip(79)=ievals
        ip(80)=itnum
        if(iprob.eq.3) ip(80)=itnum-1
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        integer function icvtst(itnum,iprob,itask,itype,rp)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            double precision
     +          rp(100)
            save tola,tolb,eps,erf,egf,tole,tolr,isw,trf
            data ibit/0/
c
c       convergence test for outer newton loop
c
c       icvtst = -1  making progress
c       icvtst =  0  not converged
c       icvtst =  1  converged
c
        ii=0
        if(iabs(iprob).ne.3.or.itask.ge.5) ii=1
        if(iprob.lt.0) ii=1
        if(itype.lt.0) ii=2
c
        if(itnum.le.1) then
            isw=0
            eps=1.0d2*ceps(ibit)
            tola=eps
            if(itype.eq.1) tola=dsqrt(tola)
            tolb=tola
            trf=0.5d0
            erf=1.0d0-eps
            egf=0.1d0
            if(ii.eq.1) then
                tole=1.0d-1
                tolr=1.0d-2
            else if(ii.eq.2) then
                tole=1.0d-2
                tolr=1.0d-4
            else
                tole=1.0d-2
                tolr=1.0d-4
            endif
        endif
c
        reler0=rp(53)
        relerr=rp(54)
        relres=rp(56)
        ratio=rp(57)
c
c       revise tol if indicated
c
        if(isw.eq.0.and.ii.ge.1.and.relerr.lt.trf) then
            isw=1
            tola=dmax1(relerr*tole,tola)
            if(relres.le.1.0d0) tolb=dmax1(relres*tole,tolb)
            if(reler0.lt.0.5d0) tola=dmax1(reler0*tole,tola)
        endif
c
c       convergence test
c
        icvtst=0
        if(relerr.lt.tole.or.relres.lt.tolr) icvtst=-1
        if(relerr.le.tola.and.ratio.le.erf) icvtst=-1
        if(relerr.le.tola.and.relres.le.tolb) icvtst=1
        if(relres.lt.eps.and.ratio.ge.egf) icvtst=1
        if(relerr.lt.eps) icvtst=1
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        integer function jcvtst(itnum,iprob,itask,itype,rp)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            double precision
     +          rp(100)
            save tola,tolb,eps,erf,egf,tole,tolr,isw,trf
            data ibit/0/
c
c       debug convergence test for outer parallel solve
c
c       jcvtst = -1  making progress
c       jcvtst =  0  not converged
c       jcvtst =  1  converged
c
        ii=0
        if(iabs(iprob).ne.3.or.itask.ge.5) ii=1
        if(iprob.lt.0) ii=1
        if(itype.lt.0) ii=2
c
        if(itnum.le.1) then
            isw=0
            eps=1.0d2*ceps(ibit)
            tola=eps
            if(itype.eq.1) tola=dsqrt(tola)
            tolb=tola
            trf=0.5d0
            erf=1.0d0-eps
            egf=0.1d0
            if(ii.eq.1) then
                tole=1.0d-4
                tolr=1.0d-6
            else if(ii.eq.2) then
                tole=1.0d-2
                tolr=1.0d-4
            else
                tole=1.0d-2
                tolr=1.0d-4
            endif
        endif
c
        reler0=rp(53)
        relerr=rp(54)
        relres=rp(56)
        ratio=rp(57)
c
c       revise tol if indicated
c
        if(isw.eq.0.and.ii.ge.1.and.relerr.lt.trf) then
            isw=1
            tola=dmax1(relerr*tole,tola)
            if(relres.le.1.0d0) tolb=dmax1(relres*tole,tolb)
            if(reler0.lt.0.5d0) tola=dmax1(reler0*tole,tola)
        endif
c
c       convergence test
c
        jcvtst=0
        if(relerr.lt.tole.or.relres.lt.tolr) jcvtst=-1
        if(relerr.le.tola.and.ratio.le.erf) jcvtst=-1
        if(relerr.le.tola.and.relres.le.tolb) jcvtst=1
        if(relres.lt.eps.and.ratio.ge.egf) jcvtst=1
        if(relerr.lt.eps) jcvtst=1
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine cev(ip,rp,ja,a,ka,evl,evr,br,bl,devr,devl,
     +      evr0,evl0,z,hist)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ja(*),ip(100),ka(*)
            double precision
     +          rp(100),a(*),evl(*),evr(*),bl(*),br(*),evr0(*),
     1          evl0(*),devr(*),devl(*),z(*),hist(22,*)
            save ibit
            data ibit/0/
c
c       compute approximate left and right singular vectors
c       this routine should work when evl is allocated and ispd=1
c       and also when evl=evr but ispd=0
c
        ndf=ip(5)
        ispd=ip(8)
        mxcg=ip(10)
        anorm=rp(55)
        jspd=1
        if(ispd.ne.1) jspd=-1
        i1=1
        i2=i1+ndf
        i3=i2+ndf
        iqptr=ka(3)
        japtr=ka(4)
        iaptr=ka(5)
        juptr=ka(6)
        iuptr=ka(7)
        nsum=ka(8)
c
        tol=1.0d-3
        eps=ceps(ibit)*1.0d2
        tol1=dmax1(tol*1.0d-1,dsqrt(eps))
        itmax=max0(10,mxcg)
c
c       check for null vectors
c
        sval=1.0d0
        evrn=rl2nrm(ndf,evr)
        evln=rl2nrm(ndf,evl)
        if(evrn.eq.0.0d0.or.evln.eq.0.0d0) then
            rp(25)=sval
            return
        endif
c
c       normalize initial vectors
c
        dp=rl2ip(ndf,evl,evr)
        if(dp.lt.0.0d0) evln=-evln
        do i=1,ndf
            ee=evr(i)/evrn
            evl(i)=evl(i)/evln
            evr(i)=ee
        enddo
        do i=1,ndf
            ii=ja(iqptr+i-1)
            z(i1+ii-1)=evr(i)
            z(i2+ii-1)=evl(i)
        enddo
        do i=1,ndf
            evr(i)=z(i1+i-1)
            evl(i)=z(i2+i-1)
            devr(i)=0.0d0
            devl(i)=0.0d0
            evr0(i)=0.0d0
            evl0(i)=0.0d0
        enddo
c
c       inverse iteration loop
c
        call hist1(hist(1,14),0,1.0d0)
        do itnum=1,itmax
c
c       a evr = sval * evl
c       a(transpose)  evl = sval * evr
c
            call mtxmlt(ndf,ja,a,evr,z,ispd)
            sval=rl2ip(ndf,evl,z)
            do i=1,ndf
                br(i)=sval*evl(i)-z(i)
                bl(i)=br(i)
            enddo
            brnorm=rl2nrm(ndf,br)
            dr=brnorm/(dabs(sval)+tol*anorm)
            if(ispd.ne.1) then
                call mtxmlt(ndf,ja,a,evl,z,jspd)
                svll=rl2ip(ndf,evr,z)
                do i=1,ndf
                    bl(i)=svll*evr(i)-z(i)
                enddo
                blnorm=rl2nrm(ndf,bl)
                dl=blnorm/(dabs(svll)+tol*anorm)
                dr=dmax1(dr,dl)
            endif
            call hist1(hist(1,14),itnum,dr)
            if(dr.lt.tol.and.itnum.gt.1) go to 100
c
c       add small perturbation
c
            ee=tol1*brnorm
            do i=1,ndf
                br(i)=br(i)+ee
                ee=-ee
            enddo
c
            call cycle(ispd,ka,ja(japtr),a(iaptr),
     +          ja(juptr),a(iuptr),devr,br,z,z(nsum+1))
            call csv(ndf,ja,a,evr,z(i1),devr,z(i2),evr0,z(i3),ispd)
c
            if(ispd.ne.1) then
                ee=tol1*blnorm
                do i=1,ndf
                    bl(i)=bl(i)+ee
                    ee=-ee
                enddo
c
                call cycle(jspd,ka,ja(japtr),a(iaptr),
     +              ja(juptr),a(iuptr),devl,bl,z,z(nsum+1))
                call csv(ndf,ja,a,evl,z(i1),devl,z(i2),evl0,z(i3),jspd)
            else
                do i=1,ndf
                    evl(i)=evr(i)
                    evl0(i)=evr0(i)
                    devl(i)=devr(i)
                enddo
            endif
c
        enddo
        itnum=itmax
c
c       final computation of singular value
c       sign determined such that evl * evr is positive
c
  100   dp=rl2ip(ndf,evr,evl)
        if(dp.lt.0.0d0) then
            sval=-sval
            do i=1,ndf
                evl(i)=-evl(i)
            enddo
        endif
        do i=1,ndf
             ii=ja(iqptr+i-1)
             z(i1+i-1)=evr(ii)
             z(i2+i-1)=evl(ii)
        enddo
        do i=1,ndf
             evr(i)=z(i1+i-1)
             evl(i)=z(i2+i-1)
        enddo
        rp(25)=sval
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine csv(n,ja,a,ev,aev,dev,adev,ev0,aev0,ispd)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ja(*)
            double precision
     +          a(*),ev(*),aev(*),dev(*),adev(*),ev0(*),aev0(*),
     1          aa(3,3),r(3),q(3,3)
c
c       orthogonalize
c
        call orthog(n,ev,dev,ev0,irank)
c
        call mtxmlt(n,ja,a,ev,aev,ispd)
        call mtxmlt(n,ja,a,dev,adev,ispd)
        call mtxmlt(n,ja,a,ev0,aev0,ispd)
c
c       compute inner products for quadratic equation
c
        aa(1,1)=rl2ip(n,aev,aev)
        aa(1,2)=rl2ip(n,aev,adev)
        aa(1,3)=rl2ip(n,aev,aev0)
        aa(2,1)=aa(1,2)
        aa(2,2)=rl2ip(n,adev,adev)
        aa(2,3)=rl2ip(n,adev,aev0)
        aa(3,1)=aa(1,3)
        aa(3,2)=aa(2,3)
        aa(3,3)=rl2ip(n,aev0,aev0)
        call ev3x3(aa,r,q,irank)
c
c       reset ev
c
        do i=1,n
            s=q(2,1)*dev(i)+q(3,1)*ev0(i)
            ev(i)=q(1,1)*ev(i)+s
            ev0(i)=s
        enddo
        evnorm=rl2nrm(n,ev)
        if(evnorm.gt.0.0d0) evnorm=1.0d0/evnorm
        ev0nrm=rl2nrm(n,ev0)
        if(ev0nrm.gt.0.0d0) ev0nrm=1.0d0/ev0nrm
        do i=1,n
            ev(i)=ev(i)*evnorm
            ev0(i)=ev0(i)*ev0nrm
        enddo
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine c3x3(a,b,num)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          index(3,3)
            double precision
     +          a(3,3),b(3,3)
            save index
            data index/1,2,3,2,3,1,3,1,2/
c
c       this routine solves 3 x 3 linear systems
c
c       scale rows so that largest element is one
c
        do j=1,3
            rmax=dmax1(dabs(a(j,1)),dabs(a(j,2)),dabs(a(j,3)))
            if(rmax.ne.0.0d0) rmax=1.0d0/rmax
            do k=1,3
                a(j,k)=a(j,k)*rmax
            enddo
            do k=1,num
                b(j,k)=b(j,k)*rmax
            enddo
        enddo
c
        j1=1
        if(dabs(a(1,1)).lt.dabs(a(2,1))) j1=2
        if(dabs(a(j1,1)).lt.dabs(a(3,1))) j1=3
        j2=index(2,j1)
        j3=index(3,j1)
c
        if(a(j1,1).ne.0.0d0) a(j1,1)=1.0d0/a(j1,1)
        q2=a(j2,1)*a(j1,1)
        q3=a(j3,1)*a(j1,1)
        do k=1,num
            b(j2,k)=b(j2,k)-b(j1,k)*q2
            b(j3,k)=b(j3,k)-b(j1,k)*q3
        enddo
c
        a(j2,2)=a(j2,2)-a(j1,2)*q2
        a(j2,3)=a(j2,3)-a(j1,3)*q2
        a(j3,2)=a(j3,2)-a(j1,2)*q3
        a(j3,3)=a(j3,3)-a(j1,3)*q3
c
        if(dabs(a(j2,2)).lt.dabs(a(j3,2))) j2=j3
        j3=6-j1-j2
        if(a(j2,2).ne.0.0d0) a(j2,2)=1.0d0/a(j2,2)
        q3=a(j3,2)*a(j2,2)
        do k=1,num
            b(j3,k)=b(j3,k)-b(j2,k)*q3
        enddo
        a(j3,3)=a(j3,3)-a(j2,3)*q3
c
        if(a(j3,3).ne.0.0d0) a(j3,3)=1.0d0/a(j3,3)
        do k=1,num
            x3=b(j3,k)*a(j3,3)
            x2=(b(j2,k)-a(j2,3)*x3)*a(j2,2)
            b(1,k)=(b(j1,k)-a(j1,2)*x2-a(j1,3)*x3)*a(j1,1)
            b(2,k)=x2
            b(3,k)=x3
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine ev3x3(a,r,q,irank)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            double precision
     +          a(3,3),r(3),root(3),q(3,3)
c
c       solve 3 x 3 eigenvalue problem for all special cases
c
        do i=1,3
            r(i)=0.0d0
            do j=1,3
                q(i,j)=0.0d0
            enddo
            q(i,i)=1.0d0
        enddo
        as=dmax1(a(1,1),a(2,2),a(3,3))
        if(as.eq.0.0d0) return
        a11=a(1,1)/as
        a22=a(2,2)/as
        a33=a(3,3)/as
        a12=a(1,2)/as
        a13=a(1,3)/as
        a23=a(2,3)/as
c
c       irank=1
c
        if(irank.eq.1) then
            r(1)=a(1,1)
            return
        endif
        if(irank.eq.2) then
c
c       coefficients of quadratic
c
            d12=a12**2
            b=(a11+a22)/2.0d0
            c0=dmax1(a11*a22-d12,0.0d0)
            d=(a11-a22)/2.0d0
            d=dsqrt(d*d+d12)
            r1=c0/(b+d)
            r(1)=r1*as
            r(2)=(b+d)*as
            r(3)=0.0d0
            d1=a11-r1
            d2=a22-r1
            if(dmax1(dabs(d1),dabs(d2)).eq.0.0d0) then
                c=1.0d0
                s=0.0d0
            else if(dabs(d1).gt.dabs(d2)) then
                dd=1.0d0/dsqrt(d1**2+d12)
                c=-a12*dd
                s=d1*dd
            else
                dd=1.0d0/dsqrt(d2**2+d12)
                s=-a12*dd
                c=d2*dd
            endif
            q(1,1)=c
            q(2,1)=s
            q(1,2)=-s
            q(2,2)=c
            return
        endif
c
c       coefficients of cubic polynomial
c
        tol=1.0d-3
        d12=a12**2
        d13=a13**2
        d23=a23**2
        p=-(a11+a22+a33)/3.0d0
        qq=a11*a22+a22*a33+a33*a11-d12-d13-d23
        s=a11*d23+a22*d13+a33*d12
     +      -a11*a22*a33-2.0d0*a12*a23*a13
c
c       solve cubic equation (all roots should be real and non-neg.)
c
        aa=qq/3.0d0-p**2
        bb=p**3-(p*qq-s)/2.0d0
        if(bb**2+aa**3.ge.0.0d0) then
c
c       case of two equal roots (assume b*b+a*a*a=0)
c
            sgn=2.0d0
            if(bb.gt.0.0d0) sgn=-2.0d0
            bb=sgn*(dabs(bb)**(1.0d0/3.0d0))
            r(1)=bb-p
            r(2)=-bb/2.0d0-p
            r(3)=r(2)
        else
c
c       three distinct roots
c
            d=dsqrt(-aa)*2.0d0
            theta=2.0d0*bb/(aa*d)
            theta=dmin1(1.0d0,theta)
            theta=dmax1(-1.0d0,theta)
            theta=dacos(theta)/3.0d0
            pi=3.141592653589793d0/3.0d0
            r(1)=d*dcos(theta)-p
            r(2)=d*dcos(theta+2.0d0*pi)-p
            r(3)=d*dcos(theta+4.0d0*pi)-p
        endif
c
c       order
c
        ic1=1
        if(r(2).lt.r(1)) ic1=2
        if(r(3).lt.r(ic1)) ic1=3
        ic2=(5-ic1)/2
        ic3=6-ic1-ic2
        if(r(ic3).lt.r(ic2)) ic2=ic3
        ic3=6-ic1-ic2
        root(1)=r(ic1)
        root(2)=r(ic2)
        root(3)=r(ic3)
        r(1)=root(1)*as
        r(2)=root(2)*as
        r(3)=root(3)*as
c
c       now get eigenvectors
c
        if(r(3)-r(1).lt.tol*r(3)) then
            return
        else if(dmin1(r(2)-r(1),r(3)-r(2)).le.tol*r(2)) then
            a1=a11-root(2)
            a2=a22-root(2)
            a3=a33-root(2)
            s1=a1**2+d12+d13
            s2=a2**2+d12+d23
            s3=a3**2+d13+d23
            if(s1.gt.dmax1(s2,s3)) then
                qq=1.0d0/dsqrt(s1)
                v1=qq*a1
                v2=qq*a12
                v3=qq*a13
            else if(s2.gt.s3) then
                qq=1.0d0/dsqrt(s2)
                v1=qq*a12
                v2=qq*a2
                v3=qq*a23
            else
                qq=1.0d0/dsqrt(s3)
                v1=qq*a13
                v2=qq*a23
                v3=qq*a3
            endif
            if(v1.eq.0.0d0) then
                w1=1.0d0
                w2=0.0d0
                w3=0.0d0
            else if(v2.eq.0.0d0) then
                w1=0.0d0
                w2=1.0d0
                w3=0.0d0
            else
                qq=1.0d0/dsqrt(v1**2+v2**2)
                w1=-v2*qq
                w2=v1*qq
                w3=0.0d0
            endif
            z1=v2*w3-v3*w2
            z2=v3*w1-v1*w3
            z3=v1*w2-v2*w1
            if(r(2)-r(1).le.tol*r(2)) then
                dd=dsqrt((z1-w2)**2+(z2+w1)**2)
                c=(z2+w1)/dd
                s=(z1-w2)/dd
                q(1,1)=c*w1+s*z1
                q(2,1)=c*w2+s*z2
                q(3,1)=c*w3+s*z3
                q(1,2)=c*z1-s*w1
                q(2,2)=c*z2-s*w2
                q(3,2)=c*z3-s*w3
                q(1,3)=v1
                q(2,3)=v2
                q(3,3)=v3
            else
                dd=dsqrt((z2-w3)**2+(z3+w2)**2)
                c=(z3+w2)/dd
                s=(z2-w3)/dd
                q(1,1)=v1
                q(2,1)=v2
                q(3,1)=v3
                q(1,2)=c*w1+s*z1
                q(2,2)=c*w2+s*z2
                q(3,2)=c*w3+s*z3
                q(1,3)=c*z1-s*w1
                q(2,3)=c*z2-s*w2
                q(3,3)=c*z3-s*w3
            endif
        else
c
c       the general case
c
c
            if(r(2)-r(1).gt.(r(3)-r(2))*1.d-2) then
                js=1
                jf=2
            else
                js=2
                jf=3
            endif
            do i=js,jf
                a1=a11-root(i)
                a2=a22-root(i)
                a3=a33-root(i)
                v1=a2*a3-d23
                v2=a13*a23-a12*a3
                v3=a12*a23-a13*a2
                vv=v1**2+v2**2+v3**2
                w1=v2
                w2=a1*a3-d13
                w3=a13*a12-a23*a1
                ww=w1**2+w2**2+w3**2
                z1=v3
                z2=w3
                z3=a1*a2-d12
                zz=z1**2+z2**2+z3**2
                if(vv.gt.dmax1(ww,zz)) then
                    qq=1.0d0/dsqrt(vv)
                    q(1,i)=qq*v1
                    q(2,i)=qq*v2
                    q(3,i)=qq*v3
                else if(ww.gt.zz) then
                    qq=1.0d0/dsqrt(ww)
                    q(1,i)=qq*w1
                    q(2,i)=qq*w2
                    q(3,i)=qq*w3
                else
                    qq=1.0d0/dsqrt(zz)
                    q(1,i)=qq*z1
                    q(2,i)=qq*z2
                    q(3,i)=qq*z3
                endif
            enddo
            ic=6-js-jf
            q(1,ic)=q(2,js)*q(3,jf)-q(3,js)*q(2,jf)
            q(2,ic)=q(3,js)*q(1,jf)-q(1,js)*q(3,jf)
            q(3,ic)=q(1,js)*q(2,jf)-q(2,js)*q(1,jf)
        endif
        do i=1,3
            if(q(i,i).lt.0.0d0) then
                q(1,i)=-q(1,i)
                q(2,i)=-q(2,i)
                q(3,i)=-q(3,i)
            endif
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine orthog(n,v1,v2,v3,irank)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            double precision
     +          v1(*),v2(*),v3(*),r(3)
c
c       orthogonalize, normalize, determine rank
c
        tol=1.d-1
        a11=0.0d0
        a22=0.0d0
        a33=0.0d0
        do i=1,n
            a11=a11+v1(i)**2
            a22=a22+v2(i)**2
            a33=a33+v3(i)**2
        enddo
        if(a11.gt.0.0d0) a11=1.0d0/dsqrt(a11)
        if(a22.gt.0.0d0) a22=1.0d0/dsqrt(a22)
        if(a33.gt.0.0d0) a33=1.0d0/dsqrt(a33)
        d12=0.0d0
        d13=0.0d0
        do i=1,n
            v1(i)=v1(i)*a11
            v2(i)=v2(i)*a22
            v3(i)=v3(i)*a33
            d12=d12+v1(i)*v2(i)
            d13=d13+v1(i)*v3(i)
        enddo
        a22=0.0d0
        a33=0.0d0
        do i=1,n
            v2(i)=v2(i)-d12*v1(i)
            v3(i)=v3(i)-d13*v1(i)
            a22=a22+v2(i)**2
            a33=a33+v3(i)**3
        enddo
        if(a22.gt.0.0d0) a22=1.0d0/dsqrt(a22)
        if(a33.gt.0.0d0) a33=1.0d0/dsqrt(a33)
        d23=0.0d0
        do i=1,n
            v2(i)=v2(i)*a22
            v3(i)=v3(i)*a33
            d23=d23+v2(i)*v3(i)
        enddo
        a33=0.0d0
        do i=1,n
            v3(i)=v3(i)-d23*v2(i)
            a33=a33+v3(i)**2
        enddo
        if(a33.gt.0.0d0) a33=1.0d0/dsqrt(a33)
        a12=0.0d0
        a13=0.0d0
        a23=0.0d0
        do i=1,n
            v3(i)=v3(i)*a33
            a12=a12+v1(i)*v2(i)
            a13=a13+v1(i)*v3(i)
            a23=a23+v2(i)*v3(i)
        enddo
c
c       coefficients of cubic polynomial
c
        if(a11.gt.0.0d0) a11=1.0d0
        if(a22.gt.0.0d0) a22=1.0d0
        if(a33.gt.0.0d0) a33=1.0d0
        d12=a12**2
        d13=a13**2
        d23=a23**2
        p=-(a11+a22+a33)/3.0d0
        qq=a11*a22+a22*a33+a33*a11-d12-d13-d23
        s=a11*d23+a22*d13+a33*d12
     +      -a11*a22*a33-2.0d0*a12*a23*a13
c
c       solve cubic equation (all roots should be real and non-neg.)
c
        aa=qq/3.0d0-p**2
        bb=p**3-(p*qq-s)/2.0d0
        if(bb**2+aa**3.ge.0.0d0) then
c
c       case of two equal roots (assume b*b+a*a*a=0)
c
            sgn=2.0d0
            if(bb.gt.0.0d0) sgn=-2.0d0
            bb=sgn*(dabs(bb)**(1.0d0/3.0d0))
            r(1)=bb-p
            r(2)=-bb/2.0d0-p
            r(3)=r(2)
        else
c
c       three distinct roots
c
            d=dsqrt(-aa)*2.0d0
            theta=2.0d0*bb/(aa*d)
            theta=dmin1(1.0d0,theta)
            theta=dmax1(-1.0d0,theta)
            theta=dacos(theta)/3.0d0
            pi=3.141592653589793d0/3.0d0
            r(1)=d*dcos(theta)-p
            r(2)=d*dcos(theta+2.0d0*pi)-p
            r(3)=d*dcos(theta+4.0d0*pi)-p
        endif
c
c       order
c
        ic1=1
        if(r(2).lt.r(1)) ic1=2
        if(r(3).lt.r(ic1)) ic1=3
        ic2=(5-ic1)/2
        ic3=6-ic1-ic2
        if(r(ic3).lt.r(ic2)) ic2=ic3
        ic3=6-ic1-ic2
c
        irank=1
        if(r(ic2).gt.tol) irank=2
        if(r(ic1).gt.tol) irank=3
c
        if(irank.eq.1) then
            do i=1,n
                v2(i)=0.0d0
                v3(i)=0.0d0
            enddo
        else if(irank.eq.2.and.a33.gt.0.0d0) then
            if(a22.le.0.0d0) then
                do i=1,n
                    v2(i)=v3(i)
                enddo
            else if(dabs(a13).lt.dabs(a12)) then
                do i=1,n
                    v2(i)=v3(i)
                enddo
            endif
            do i=1,n
                v3(i)=0.0d0
            enddo
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine grad(ux,uy,vx,vy,u,iv,isw)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          iv(3)
            double precision
     +          vx(*),vy(*),u(*)
c
c       compute the gradient of u in element defined by iv
c
        iv1=iv(1)
        iv2=iv(2)
        iv3=iv(3)
        x2=vx(iv2)-vx(iv1)
        x3=vx(iv3)-vx(iv1)
        y2=vy(iv2)-vy(iv1)
        y3=vy(iv3)-vy(iv1)
        if(isw.eq.1) then
            u2=u(2)-u(1)
            u3=u(3)-u(1)
        else
            u2=u(iv2)-u(iv1)
            u3=u(iv3)-u(iv1)
        endif
        det=x2*y3-x3*y2
        ux=(u2*y3-u3*y2)/det
        uy=(x2*u3-x3*u2)/det
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        double precision function carea(ntf,itnode,itedge,
     +      ibndry,vx,vy,xm,ym)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          itnode(5,*),ibndry(6,*),itedge(3,*),index(3,3)
            double precision
     +          vx(*),vy(*),xm(*),ym(*),c(3)
            save index
            data index/1,2,3,2,3,1,3,1,2/
c
c       compute area of domain
c
        carea=0.0d0
        pi=3.141592653589793d0
        do i=1,ntf
            x2=vx(itnode(2,i))-vx(itnode(1,i))
            x3=vx(itnode(3,i))-vx(itnode(1,i))
            y2=vy(itnode(2,i))-vy(itnode(1,i))
            y3=vy(itnode(3,i))-vy(itnode(1,i))
            det=dabs(x2*y3-x3*y2)/2.0d0
c
c       curved edges
c
            do 5 j=1,3
                if(itedge(j,i).ge.0) go to 5
                k=-itedge(j,i)
                if(ibndry(3,k).le.0) go to 5
                kt=ibndry(3,k)
                iv1=itnode(index(2,j),i)
                iv2=itnode(index(3,j),i)
                call arc(vx(iv1),vy(iv1),vx(iv2),vy(iv2),
     +              xm(kt),ym(kt),theta1,theta2,rad,alen)
                call bari(xm(kt),ym(kt),vx,vy,itnode(1,i),c)
                theta=dabs(theta2-theta1)*pi
                aa=(rad**2/2.0d0)*(theta-dsin(theta))
                if(c(j).lt.0.0d0) aa=-aa
                det=det+aa
    5       enddo
            carea=carea+det
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine bari(x,y,vx,vy,iv,c)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          iv(3)
            double precision
     +          vx(*),vy(*),c(3)
c
c       compute the barycentric coordinates of the point (x,y)
c
        iv1=iv(1)
        iv2=iv(2)
        iv3=iv(3)
        x2=vx(iv2)-vx(iv1)
        y2=vy(iv2)-vy(iv1)
        x3=vx(iv3)-vx(iv1)
        y3=vy(iv3)-vy(iv1)
        xr=x-vx(iv1)
        yr=y-vy(iv1)
        det=x2*y3-x3*y2
        c(2)=(xr*y3-x3*yr)/det
        c(3)=(x2*yr-xr*y2)/det
        c(1)=1.0d0-c(2)-c(3)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        double precision function rl2nrm(n,b)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            double precision
     +          b(*)
c
c       compute norm of b and update history
c
        bnorm=0.0d0
        bmax=0.0d0
        do i=1,n
            if(dabs(b(i)).lt.bmax) then
                bnorm=bnorm+(b(i)/bmax)**2
            else if(b(i).ne.0.0d0) then
                bnorm=1.0d0+bnorm*(bmax/b(i))**2
                bmax=dabs(b(i))
            endif
        enddo
        rl2nrm=dsqrt(bnorm)*bmax
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        double precision function rl2ip(n,x,y)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            double precision
     +          x(*),y(*)
c
c       compute dot product
c
        rl2ip=0.0d0
        spmax=0.0d0
        snmax=0.0d0
        sp=0.0d0
        sn=0.0d0
        do i=1,n
            t=x(i)*y(i)
            if(t.ge.0.0d0) then
                if(t.lt.spmax) then
                    sp=sp+t/spmax
                else if(t.ne.0.0d0) then
                    sp=1.0d0+sp*(spmax/t)
                    spmax=t
                endif
            else
                if(-t.lt.snmax) then
                    sn=sn+t/snmax
                else
                    sn=-(1.0d0+sn*(snmax/t))
                    snmax=-t
                endif
            endif
        enddo
        rl2ip=sp*spmax+sn*snmax
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine pnorm(ip,rp,itnode,vx,vy,lenb,bump,maxd,nef,
     +      u,mark,ndof,itdof,hist)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100),itnode(5,*),mark(*),idof(10),itdof(ndof,*)
            double precision
     +          bump(lenb,*),rp(100),u(maxd,*),vx(*),vy(*),t(6),
     1          hist(22,*)
c
c       compute global errors estimate on distributed mesh
c
        irgn=ip(50)
        ntf=ip(1)
        ndf=ip(5)
        iord=ip(26)
c
        do i=1,ndf
            mark(i)=0
        enddo
        enorm1=0.0d0
        enorm2=0.0d0
        unorm1=0.0d0
        unorm2=0.0d0
        do i=1,ntf
            if(itnode(4,i).eq.irgn) then
                call l2gmap(i,idof,ndof,itdof)
                do j=1,ndof
                    mark(idof(j))=1
                enddo
                e1=tqual(i,itnode,vx,vy,lenb,bump,iord)
                e2=tqual2(i,itnode,vx,vy,lenb,bump,iord)
                enorm1=enorm1+e1
                enorm2=enorm2+e2
c
                do ifn=1,nef
                    unorm1=unorm1+eh1nrm(i,itnode,vx,vy,
     +                  u(1,ifn),idof,iord)
                    unorm2=unorm2+el2nrm(i,itnode,vx,vy,
     +                  u(1,ifn),idof,iord)
                 enddo
            endif
        enddo
        ndg=0
        do i=1,ndf
            if(mark(i).eq.1) ndg=ndg+1
        enddo
        t(1)=unorm2
        t(2)=unorm1
        t(3)=enorm2
        t(4)=enorm1
        t(5)=dfloat(ndg)
c
        call pl2ip(t,5)
c
        enorm1=dsqrt(t(4))
        rp(37)=enorm1
        unorm1=dsqrt(t(2))
        rp(38)=unorm1
        rp(39)=dsqrt(t(3))
        rp(40)=dsqrt(t(1))
        ndg=idint(t(5))
        relerr=1.0d0
        if(unorm1.ne.0.0d0) relerr=enorm1/unorm1
        if(unorm1+enorm1.le.0.0d0) relerr=0.0d0
        rp(53)=relerr
c
        call hist2(hist,rp,-2,ndg)
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine l2mtx1(nvf,ntf,vx,vy,itnode,ja,a)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          itnode(5,*),ja(*)
            double precision
     +          vx(*),vy(*),a(*)
c
c       mass matrix for linear elements with node numbering
c
        do i=1,ja(nvf+1)-1
            a(i)=0.0d0
        enddo
        ishift=0
        do i=1,ntf
            iv1=itnode(1,i)
            iv2=itnode(2,i)
            iv3=itnode(3,i)
            x2=vx(iv2)-vx(iv1)
            x3=vx(iv3)-vx(iv1)
            y2=vy(iv2)-vy(iv1)
            y3=vy(iv3)-vy(iv1)
            det=dabs(x2*y3-x3*y2)/24.0d0
            do k=1,3
                ivk=itnode(k,i)
                a(ivk)=a(ivk)+2.0d0*det
                do j=k+1,3
                    call jamap(ivk,itnode(j,i),kj,jk,ja,ishift)
                    a(jk)=a(jk)+det
                enddo
            enddo
        enddo
        anorm=0.0d0
        do i=1,nvf
            anorm=dmax1(anorm,dabs(a(i)))
        enddo
        do i=1,nvf
            if(a(i).eq.0.0d0) a(i)=anorm
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine h1mtx1(nvf,ntf,vx,vy,itnode,ja,a)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          itnode(5,*),ja(*)
            double precision
     +          vx(*),vy(*),a(*),tx(3),ty(3),x(3),y(3)
c
c       stiffness matrix for linear elements with node numbering
c
        do i=1,ja(nvf+1)-1
            a(i)=0.0d0
        enddo
        ishift=0
        do i=1,ntf
            call afmap(i,itnode,vx,vy,tx,ty,x,y,det)
            det=dabs(det)/2.0d0
            do k=1,3
                ivk=itnode(k,i)
                a(ivk)=a(ivk)+det*(x(k)**2+y(k)**2)
                do j=k+1,3
                    call jamap(ivk,itnode(j,i),kj,jk,ja,ishift)
                    a(jk)=a(jk)+det*(x(k)*x(j)+y(k)*y(j))
                enddo
            enddo
        enddo
        anorm=0.0d0
        do i=1,nvf
            anorm=dmax1(anorm,dabs(a(i)))
        enddo
        do i=1,nvf
            if(a(i).eq.0.0d0) a(i)=anorm
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine l2mtx(n,ntf,vx,vy,itnode,ja,a,iord,ndof,itdof)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          itnode(5,*),itdof(ndof,*),ja(*),idof(10)
            double precision
     +          vx(*),vy(*),a(*),ea(10,10)
c
c
        do i=1,ja(n+1)-1
            a(i)=0.0d0
        enddo
c
c       compute mass matrix
c
        ishift=0
        do i=1,ntf
            call l2gmap(i,idof,ndof,itdof)
            call elel2(i,itnode,vx,vy,ea,iord)
            do k=1,ndof
                ivk=idof(k)
                a(ivk)=a(ivk)+ea(k,k)
                do j=k+1,ndof
                    call jacmap(ivk,idof(j),kj,jk,ja,ishift)
                    a(jk)=a(jk)+ea(j,k)
                enddo
            enddo
        enddo
        anorm=0.0d0
        do i=1,n
            anorm=dmax1(anorm,dabs(a(i)))
        enddo
        do i=1,n
            if(a(i).eq.0.0d0) a(i)=anorm
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine h1mtx(n,ntf,vx,vy,itnode,ja,a,iord,ndof,itdof)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          itnode(5,*),itdof(ndof,*),ja(*),idof(10)
            double precision
     +          vx(*),vy(*),a(*),ea(10,10)
c
c
        do i=1,ja(n+1)-1
            a(i)=0.0d0
        enddo
c
c       compute stiffness matrix for laplacian
c
        ishift=0
        do i=1,ntf
            call l2gmap(i,idof,ndof,itdof)
            call eleh1(i,itnode,vx,vy,ea,iord)
            do k=1,ndof
                ivk=idof(k)
                a(ivk)=a(ivk)+ea(k,k)
                do j=k+1,ndof
                    call jacmap(ivk,idof(j),kj,jk,ja,ishift)
                    a(jk)=a(jk)+ea(j,k)
                enddo
            enddo
        enddo
        anorm=0.0d0
        do i=1,n
            anorm=dmax1(anorm,dabs(a(i)))
        enddo
        do i=1,n
            if(a(i).eq.0.0d0) a(i)=anorm
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine sgscg(n,ja,a,dr,b,mxcg,eps,p,ap,z,hist,iflag)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ja(*)
            double precision
     +          a(*),hist(*),dr(*),b(*),p(*),ap(*),z(*)
c
c       initialize
c
        iflag=0
        ispd=1
        epsmin=0.5d0
        relerr=0.0d0
c
c       compute initial residual
c
        call mtxmlt(n,ja,a,dr,ap,ispd)
        do i=1,n
            b(i)=b(i)-ap(i)
        enddo
c
c       compute initial norm of b
c
        bnorm=rl2nrm(n,b)
        call hist1(hist,0,bnorm)
        if(bnorm.le.0.0d0) return
        rnorm=bnorm
c
c       compute initial p and ap
c
        call sgs(n,ja,a,p,b,ispd)
        call mtxmlt(n,ja,a,p,ap,ispd)
        bp=rl2ip(n,p,b)
c
c       the main loop
c
        do itnum=1,mxcg
c
c       compute alpha and precondition
c
            pap=rl2ip(n,p,ap)
            alpha=bp/pap
            do i=1,n
                dr(i)=dr(i)+alpha*p(i)
                b(i)=b(i)-alpha*ap(i)
            enddo
            call sgs(n,ja,a,z,b,ispd)
c
c       compute coefficients
c
            bz=rl2ip(n,z,b)
            beta=bz/bp
            bp=bz
            do i=1,n
                p(i)=z(i)+beta*p(i)
            enddo
            call mtxmlt(n,ja,a,p,ap,ispd)
c
c       convergence test
c
            rnorm=rl2nrm(n,b)
            call hist1(hist,itnum,rnorm)
            relerr=rnorm/bnorm
            if(relerr.le.eps) return
            call mtxmlt(n,ja,a,p,ap,ispd)
        enddo
        if(relerr.gt.epsmin) iflag=10
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine sgs(n,ja,a,x,b,ispd)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ja(*),lmtx,umtx
            double precision
     +          a(*),x(*),b(*)
c
c       ispd = 1   symmetric
c            = 0   non-symmetric
c            =-1   non-symmetric for a-transpose
c
        lmtx=0
        umtx=0
        if(ispd.eq.0) lmtx=ja(n+1)-ja(1)
        if(ispd.eq.-1) umtx=ja(n+1)-ja(1)
c
c
c       solve   sgs * x = b
c
        do i=1,n
            x(i)=b(i)
        enddo
c
c       the lower triangular system
c
        do i=1,n
            s=x(i)/a(i)
            if(ja(i).lt.ja(i+1)) then
                do jj=ja(i),ja(i+1)-1
                    j=ja(jj)
                    x(j)=x(j)-a(jj+lmtx)*s
                enddo
            endif
        enddo
c
c       the upper triangular system
c
        do i=n,1,-1
            s=0.0d0
            if(ja(i).lt.ja(i+1)) then
                do jj=ja(i),ja(i+1)-1
                    j=ja(jj)
                    s=s+a(jj+umtx)*x(j)
                enddo
            endif
            x(i)=(x(i)-s)/a(i)
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine sgscg1(n,ja,a,x,r,mxcg,ap,p,z,eps)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ja(*)
            double precision
     +          a(*),ap(*),p(*),z(*),x(*),r(*)
c
c       sgs-cg using just one matrix multiply per iteration
c
c       initialize
c
        zdz=0.0d0
        relerr=1.0d0
        do i=1,n
            p(i)=0.0d0
            ap(i)=0.0d0
            sum=a(i)*x(i)
            do j=ja(i),ja(i+1)-1
                sum=sum+x(ja(j))*a(j)
            enddo
            r(i)=r(i)-sum
            z(i)=r(i)
        enddo
c
c       the main loop
c
        do itnum=1,mxcg
c
c       forward sweep
c
            sum=0.0d0
            do i=1,n
                t=z(i)/a(i)
                sum=sum+t*z(i)
                do j=ja(i),ja(i+1)-1
                    z(ja(j))=z(ja(j))-(t+x(i))*a(j)
                enddo
            enddo
c
c       test for convergence
c
            if(itnum.gt.1) then
                if(zdz.eq.0.0d0) return
                beta=sum/zdz
                relerr=relerr*beta
                if(dsqrt(relerr).lt.eps) return
            else
                beta=0.0d0
            endif
            zdz=sum
c
c       backward sweep
c
            pap=0.0d0
            do i=n,1,-1
                ap(i)=z(i)+beta*ap(i)
                sum=0.0d0
                do j=ja(i),ja(i+1)-1
                    sum=sum+z(ja(j))*a(j)
                enddo
                z(i)=(z(i)-sum)/a(i)
                p(i)=z(i)+beta*p(i)
                pap=pap+p(i)*(2.0d0*ap(i)-p(i)*a(i))
            enddo
            if(pap.eq.0.0d0) return
            alpha=zdz/pap
c
c       update x,r
c
            do i=1,n
                x(i)=x(i)+alpha*p(i)
                r(i)=r(i)-alpha*ap(i)
                z(i)=r(i)
            enddo
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine jcg(n,ja,a,x,r,mxcg,ap,p,z,eps)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ja(*)
            double precision
     +          a(*),x(*),r(*),p(*),ap(*),z(*)
c
c       cg with identity preconditioner
c
c       initialize
c
        ispd=1
        zdz=0.0d0
        relerr=1.0d0
        call mtxmlt(n,ja,a,x,ap,ispd)
        do i=1,n
            r(i)=r(i)-ap(i)
            p(i)=0.0d0
            ap(i)=0.0d0
        enddo
c
c       the main loop
c
        do itnum=1,mxcg
c
c       compute alpha and pecondition
c
            do i=1,n
cc              z(i)=r(i)/a(i)
                z(i)=r(i)
            enddo
            sum=rl2ip(n,z,r)
            if(itnum.gt.1) then
                if(zdz.eq.0.0d0) return
                beta=sum/zdz
                relerr=relerr*beta
                if(dsqrt(relerr).lt.eps) return
            else
                beta=0.0d0
            endif
            zdz=sum
            do i=1,n
                p(i)=z(i)+beta*p(i)
            enddo
            call mtxmlt(n,ja,a,p,ap,ispd)
            pap=rl2ip(n,p,ap)
            if(pap.eq.0.0d0) return
            alpha=zdz/pap
            do i=1,n
                x(i)=x(i)+alpha*p(i)
                r(i)=r(i)-alpha*ap(i)
            enddo
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        double precision function dl2nrm(n,b,d,isw)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            double precision
     +          b(*),d(*)
c
c       compute norm of b and update history
c
        bnorm=0.0d0
        bmax=0.0d0
        do i=1,n
            dd=0.0d0
            if(isw.ge.0) then
                dd=d(i)
            else
                if(d(i).ne.0.0d0) dd=1.0d0/d(i)
            endif
            if(dabs(b(i)).lt.bmax) then
                bnorm=bnorm+dd*(b(i)/bmax)**2
            else if(b(i).ne.0.0d0) then
                bnorm=dd+bnorm*(bmax/b(i))**2
                bmax=dabs(b(i))
            endif
        enddo
        dl2nrm=dsqrt(bnorm)*bmax
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        double precision function dl2ip(n,x,y,d,isw)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            double precision
     +          x(*),y(*),d(*)
c
c       compute dot product
c
        dl2ip=0.0d0
        spmax=0.0d0
        snmax=0.0d0
        sp=0.0d0
        sn=0.0d0
        do i=1,n
            t=0.0d0
            if(isw.ge.0) then
                t=x(i)*y(i)*d(i)
            else
                if(d(i).ne.0.0d0) t=x(i)*y(i)/d(i)
            endif
            if(t.ge.0.0d0) then
                if(t.lt.spmax) then
                    sp=sp+t/spmax
                else if(t.ne.0.0d0) then
                    sp=1.0d0+sp*(spmax/t)
                    spmax=t
                endif
            else
                if(-t.lt.snmax) then
                    sn=sn+t/snmax
                else
                    sn=-(1.0d0+sn*(snmax/t))
                    snmax=-t
                endif
            endif
        enddo
        dl2ip=sp*spmax+sn*snmax
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine mkgm(ndf,ntf,vx,vy,gm,itnode,ndof,itdof)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          itnode(5,*),itdof(ndof,*),idof(10)
            double precision
     +          vx(*),vy(*),gm(*),bw(10)
            save bw
            data bw/1.0d0,1.0d0,1.0d0,3.0d0,3.0d0,3.0d0,3.0d0,3.0d0,
     +          3.0d0,6.0d0/
c
c       compute diag of gram matrix
c
        do i=1,ndf
            gm(i)=0.0d0
        enddo
        do i=1,ntf
            x2=vx(itnode(2,i))-vx(itnode(1,i))
            y2=vy(itnode(2,i))-vy(itnode(1,i))
            x3=vx(itnode(3,i))-vx(itnode(1,i))
            y3=vy(itnode(3,i))-vy(itnode(1,i))
            det=dabs(x2*y3-x3*y2)/12.0d0
c
            call l2gmap(i,idof,ndof,itdof)
            do j=1,ndof
                gm(idof(j))=gm(idof(j))+det*bw(j)
            enddo
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine centre(x1,y1,x2,y2,x3,y3,xc,yc)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
c
c       compute the center of the circle which passes
c       through (x1,y1), (x2,y2), and (x3,y3)
c
        z1=x1-x3
        z2=x2-x3
        w1=y1-y3
        w2=y2-y3
        det=(z1*w2-z2*w1)*2.0d0
        if(det.ne.0.0d0) then
            r1=(z1*(x1+x3)+w1*(y1+y3))/det
            r2=(z2*(x2+x3)+w2*(y2+y3))/det
            xc=r1*w2-r2*w1
            yc=z1*r2-z2*r1
        else
            xc=x1
            yc=y1
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine midpt(x1,y1,x2,y2,xc,yc,x,y)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
c
c       compute the midpoint of the circle with center (xc,yc)
c       which passes through the points (x1,y1),(x2,y2).
c       the midpoint (x,y) is relative to the shorter
c       of the two arcs.
c
        x=(x1+x2)/2.0d0
        y=(y1+y2)/2.0d0
        c=(x+x1-2.0d0*xc)*(x1-x)+(y+y1-2.0d0*yc)*(y1-y)
        if(c.le.0.0d0) return
        dy=y1-y2
        dx=x1-x2
        b=(x-xc)*dy-(y-yc)*dx
        a=dx*dx+dy*dy
        a=c/(dabs(b)+dsqrt(b*b+a*c))
        if(b.lt.0.0d0)  a=-a
        x=x+a*dy
        y=y-a*dx
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        double precision function ceps(ibit)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          eptst
            save isw,sveps,jbit
            data isw,jbit,sveps/1,0,0.0d0/
c
c       compute machine epsilon
c
        if(isw.eq.0) then
            ceps=sveps
            ibit=jbit
            return
        else
            ibit=-4
            eps=1.0d0
    3       eps1=1.0d0+eps
            if(eptst(eps1).eq.1) then
                ceps=2.0d0**(-ibit)
                sveps=ceps
                jbit=ibit
                isw=0
                return
            else
                eps=eps/2.0d0
                ibit=ibit+1
                go to 3
            endif
        endif
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        integer function eptst(x)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
c
c       this is to force a store of eps1 to memory
c
        if(x.eq.1.0d0) then
            eptst=1
        else
            eptst=0
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine gfptr(itype,itask,maxd,iuu,iu0,iudot,iu0dot,
     +          ievr,ievl,ivx0,ivy0,ium,iuc,ngf,nef)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
c
        iprob=iabs(itype)
c
c       set grid function pointers
c
        iuu=1
c
        iu0=1
        iudot=1
        iu0dot=1
        ievr=1
        ievl=1
c
        ivx0=1
        ivy0=1
c
        ium=1
        iuc=1
c
        ngf=1
        nef=1
c
        if(iprob.eq.1) then
            ium=iuu+maxd
            ngf=2
            if(itask.eq.9) nef=2
        else if(iprob.eq.3) then
            iu0=iuu+maxd
            iudot=iu0+maxd
            iu0dot=iudot+maxd
            ievr=iu0dot+maxd
            ievl=ievr+maxd
            ngf=6
c
c       ngf=6 in case the user changes ispd=1 to ispd=0 after init.
c
        else if(iprob.eq.4) then
            ium=iuu+maxd
            iudot=ium+maxd
            ngf=3
            nef=2
        else if(iprob.eq.5) then
            ium=iuu+maxd
            iuc=ium+maxd
            ngf=3
            nef=3
        else if(iprob.eq.6) then
            iu0=iuu+maxd
            ivx0=iu0+maxd
            ivy0=ivx0+maxd
            ngf=4
            nef=2
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine lsptr(iprob,ispd,iord,ndf,ndd,nvdd,ibeg,iend,
     +      ihh,isu,ism,ia0,ih0,ig0,isu0,ism0,ibb,idd,
     1      ird,ipp,idl,ibdlwr,ibdupr,idu,idum,iduc)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
c
c       matrices
c
        ihh=ibeg
        isu=ibeg
        ism=ibeg
        maxjas=ndf*4
        if(iord.eq.2) maxjas=ndf*25/4
        if(iord.eq.3) maxjas=ndf*81/9
        maxjan=2*maxjas-ndf
c
c       interface matrices
c
        ia0=ibeg
        ih0=ibeg
        ig0=ibeg
        isu0=ibeg
        ism0=ibeg
c
c       vectors
c
        ibb=ibeg
        idd=ibeg
        ird=ibeg
        ipp=ibeg
        idl=ibeg
        ibdlwr=ibeg
        ibdupr=ibeg
c
c       increments
c
        idu=ibeg
        idum=ibeg
        iduc=ibeg
c
c
c
        if(iabs(iprob).eq.1) then
            ibb=ibeg
            ipp=ibb+ndf
            idu=ipp+ndf
            idum=idu+ndf
            iend=idum+ndf
        else if(iabs(iprob).eq.2) then
            ibb=ibeg
            idu=ibb+ndf
            ibdlwr=idu+ndf
            ibdupr=ibdlwr+ndf
            iend=ibdupr+ndf
        else if(iabs(iprob).eq.3) then
            ibb=ibeg
            ipp=ibb+ndf
            if(iprob.eq.3) then
                idd=ipp+ndf
                idu=idd+ndf
            else
                idd=ipp+ndf+ndd
                idu=idd+ndf+ndd
            endif
            ird=idu+ndf
            idum=ird+ndf
            iend=idum+ndf
        else if(iabs(iprob).eq.4) then
            ibb=ibeg
            idd=ibb+ndf
            ird=idd+ndf
            if(iprob.eq.4) then
                idl=ird+ndf
                ipp=idl+ndf
            else
                idl=ird+ndf+ndd
                ipp=idl+ndf+ndd
            endif
            idu=ipp+ndf
            idum=idu+ndf
            ihh=idum+ndf
            iend=ihh+maxjas
        else if(iabs(iprob).eq.5) then
            ibb=ibeg
            ipp=ibb+ndf
            idl=ipp+ndf
            ibdlwr=idl+ndf
            ibdupr=ibdlwr+ndf
            idu=ibdupr+ndf
            idum=idu+ndf
            iduc=idum+ndf
            ihh=iduc+ndf
            isu=ihh+maxjas
            ism=isu+maxjan
            iend=ism+maxjan
        else if(iabs(iprob).eq.6) then
            ibb=ibeg
            idu=ibb+ndf
            iend=idu+ndf
        endif
c
        if(iprob.gt.0) return
c
        maxja0=9*iord*nvdd/2
        maxa0s=maxja0
        maxa0n=2*maxja0-nvdd
        ia0=iend
        if(ispd.eq.1) then
            iend=ia0+maxa0s
        else
            iend=ia0+maxa0n
        endif
c
        if(iabs(iprob).eq.4) then
            ih0=iend
            iend=ih0+maxa0s
        else if(iabs(iprob).eq.5) then
            ih0=iend
            ig0=ih0+maxa0s
            isu0=ig0+maxa0s
            ism0=isu0+maxa0n
            iend=ism0+maxa0n
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine stor(ip)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100)
c
c       set up /val*/ common blocks
c
        call setval
c
c       storage allocation
c
        iprob=iabs(ip(7))
        itask=ip(9)
        lenw=ip(82)
        maxt=ip(83)
        maxv=ip(84)
        nproc=ip(49)
        ifirst=ip(6)
        iord=ifirst
        ip(26)=iord
        ndof=(iord+1)*(iord+2)/2
        maxd=maxv*iord**2
        ip(89)=maxd
        if(nproc.gt.1) then
            ss=dsqrt(dfloat(maxv))
            maxpth=idint(7.5d0*ss)*nproc
        else
            maxpth=0
        endif
        ip(81)=maxpth
c
c       determine grid functions
c
        call gfptr(iprob,itask,maxd,iuu,iu0,iudot,iu0dot,
     +          ievr,ievl,ivx0,ivy0,ium,iuc,ngf,nef)
c
c       pointers
c
        iuu=1
        if(nproc.gt.1) ngf=ngf+1
c
        itdof=iuu+ngf*maxd
        jtime=itdof+maxt*ndof
        jhist=jtime+150
        jpath=jhist+660
        ka=jpath+606
        jstat=ka+1000
        iee=jstat+10*nproc
        ipath=iee+maxt
        iz=ipath+6*maxpth
c
        ip(76)=nef
        ip(77)=ngf
c
        ip(90)=iuu
        ip(91)=itdof
        ip(92)=jtime
        ip(93)=jhist
        ip(94)=jpath
        ip(95)=ka
        ip(96)=jstat
        ip(97)=iee
        ip(98)=ipath
        ip(99)=iz
        if(iz.gt.lenw) then
            ip(25)=82
        else
            ip(25)=0
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine memptr(newptr,length,type,ibegin,iend,iflag)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          list(2,100),link(100)
            character*4 type
            save len,ifirst,link,list,level
            data len,ifirst/100,1/
c
c       this is a very crude memory manager, mainly for arrays
c       it assumes generally higest priority stuff is allocated first, low
c       priority stuff last (or from the tail, and that freeing goes in
c       reverse order from allocating.
c
        iflag=0
c
c       allocate from the head of available space
c
        if(type.eq.'head') then
            newptr=ibegin
            ibegin=ibegin+length
            if(ibegin.gt.iend+1) iflag=82
c
c       allocate from the tail of the available space
c
        elseif(type.eq.'tail') then
            iend=iend-length
            newptr=iend+1
            if(ibegin.gt.iend+1) iflag=82
c
c       save the current state of allocation (to allow a massive free)
c
        elseif(type.eq.'mark') then
            if(ifirst.eq.1) then
                ifirst=0
                level=1
                do i=1,len
                    link(i)=i+1
                enddo
                link(len)=0
            endif
            newptr=level
            if(level.gt.0) then
                level=link(level)
                list(1,newptr)=ibegin
                list(2,newptr)=iend
            else
                iflag=82
            endif
c
c       restore to state newptr
c
        elseif(type.eq.'free') then
            link(newptr)=level
            level=newptr
            ibegin=list(1,level)
            iend=list(2,level)
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine dtpick(ntf,ndf,itnode,vx,vy,u,u0,rp,z,iflag,isw,
     +      gm,ndof,itdof)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          itnode(5,*),itdof(ndof,*)
            double precision
     +          vx(*),vy(*),u(*),u0(*),rp(100),z(*),gm(*)
c
c       compute time step
c
c       iflag = -5  initialize, deltat=dtmin, next to last step
c       iflag = -4  initialize, deltat=dtmin, last step
c       iflag = -3  initialize, deltat=dtmin
c       iflag = -2  step failed, accept step (dt=dtmin)
c       iflag = -1  step failed, retake step
c       iflag =  0  normal step accepted
c       iflag =  1  next to last step
c       iflag =  2  last step
c       iflag =  3  just computed utnorm
c
c
c       initialize
c
        deltat=rp(47)
        if(isw.eq.1) then
            tcur=rp(46)
        else
            tcur=rp(46)+deltat
        endif
        dtmin=rp(48)
        dtmax=rp(49)
        utnorm=rp(50)
        tend=rp(43)
        tmtol=rp(44)
        ratio=10.0d0
        fudge=0.9d0
        iflag=3
c
c       the main loop
c
        if(isw.eq.1) go to 30
        call mkgm(ndf,ntf,vx,vy,gm,itnode,ndof,itdof)
        do i=1,ndf
            z(i)=u(i)-u0(i)
        enddo
        unorm=dl2nrm(ndf,u,gm,1)
        utnorm=dl2nrm(ndf,z,gm,1)
        if(unorm.gt.0.0d0) utnorm=utnorm/unorm
        rp(50)=utnorm
        if(isw.eq.-1) return
c
c       compute a new tentative time step
c
   30   if(utnorm.gt.tmtol) then
c
c       cut step back
c
            if(deltat.le.dtmin) then
                iflag=-2
                deltat=dtmin
            else
                deltat=dmax1(dtmin,deltat/ratio,
     +              deltat*tmtol*fudge/utnorm)
                iflag=-1
            endif
        else
            if(utnorm.gt.0.0d0) then
c
c       increase step (slight cutback if utnorm > tmtol*fudge)
c
                deltat=dmin1(dtmax,deltat*ratio,
     +              deltat*tmtol*fudge/utnorm)
                deltat=dmax1(dtmin,deltat)
                iflag=0
            else
                iflag=-3
                deltat=dtmin
            endif
        endif
c
c       check for end of interval
c
        if(tcur+deltat.ge.tend) then
            deltat=tend-tcur
            if(iflag.ne.-3) then
                iflag=2
            else
                iflag=-4
            endif
        else if(tcur+2.0d0*deltat.ge.tend) then
            if(tend-tcur-deltat.le.2.0d0*deltat/ratio)
     +          deltat=tend-tcur-2.0d0*deltat/ratio
            if(iflag.ne.-3) then
                iflag=1
            else
                iflag=-5
            endif
        endif
        rp(47)=deltat
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine bisect(rp,isw,rqup0,rqlow0)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            double precision
     +          rp(100)
            save tol,sigup,siglow,signew,sigold,
     +          rqup,rqlow,rqnew,rqold,rqmx
            data ibit/0/
c
c       this routine carries out a bisection
c       or secant iteration
c
c       isw = 0   initialize
c           > 0   update
c           < 0   converged
c
        if(isw.eq.0) then
            tol=dmax1(1.0d-6,ceps(ibit)*1.0d2)
            sigup=rp(71)
            siglow=0.0d0
            signew=sigup
            sigold=siglow
            rqup=rp(25)
            rqlow=rp(35)
            rqnew=rqup
            rqold=rqlow
            rqmx=dmax1(dabs(rqup),dabs(rqlow))
            isw=1
        else
            sigold=signew
            signew=rp(71)
            rqold=rqnew
            rqnew=rp(25)
            if(rqnew*rqlow.lt.0.0d0) then
                sigup=signew
                rqup=rqnew
            else
                siglow=signew
                rqlow=rqnew
            endif
        endif
c
c       return rqup, rqlow just for the history file
c
        rqup0=rqup
        rqlow0=rqlow
        sigma=(sigup+siglow)/2.0d0
        ds=dabs(sigup-siglow)
c
c       convergence test
c
        if(sigma.eq.signew.or.ds.lt.tol*dabs(sigma).or.
     +      dabs(rqnew).lt.tol*rqmx) then
            isw=-1
            return
        endif
c
        if(rqnew-rqold.ne.0.0d0) then
            qq=signew-rqnew*(signew-sigold)/(rqnew-rqold)
            qlow=dabs(qq-siglow)
            qup=dabs(qq-sigup)
            if(dmax1(qlow,qup).le.ds*(1.0d0-tol)) then
                sigma=qq
            else if(qlow.le.ds*tol) then
                sigma=siglow+(sigup-siglow)*tol
            else if(qup.le.ds*tol) then
                sigma=sigup+(siglow-sigup)*tol
            endif
        endif
        rp(71)=sigma
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine tpick(ip,rp,vx,vy,itnode,ndof,itdof,u,um,uc,
     +      usv,umsv,ucsv,ja,a,h,g,su,sm,b,d,p,dl,bdlwr,bdupr,
     1      du,dum,duc,gm,adu,hdu,z,isw,itnum)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100),itnode(5,*),ja(*),itdof(ndof,*)
            double precision
     +          rp(100),vx(*),vy(*),u(*),um(*),uc(*),usv(*),umsv(*),
     1          ucsv(*),a(*),h(*),g(*),su(*),sm(*),b(*),d(*),p(*),
     2          dl(*),bdlwr(*),bdupr(*),du(*),dum(*),duc(*),gm(*),
     3          adu(*),hdu(*),z(*)
            save rlsv,step0
c
c       this routine carries out a bisection
c       or secant iteration
c
c       isw = 0   initialize
c           > 0   update
c           < 0   converged
c
        ntf=ip(1)
        ndf=ip(5)
        iprob=ip(7)
        itask=ip(9)
c
c       compute norms
c
        call mkgm(ndf,ntf,vx,vy,gm,itnode,ndof,itdof)
        if(iprob.eq.1.and.itask.eq.9) then
            call norm1(ip,rp,isw,itnum,u,du,um,dum,
     +          ja,a,b,p,adu,hdu,gm,z)
            if(isw.le.0) then
                do i=1,ndf
                    usv(i)=u(i)
                    umsv(i)=um(i)
                enddo
                step0=1.0d0
            endif
        else if(iprob.eq.2) then
            call norm2(ip,rp,isw,itnum,u,du,ja,a,b,adu,gm,z)
            if(isw.le.0) then
                do i=1,ndf
                    usv(i)=u(i)
                enddo
                step0=stepmx(ndf,u,du,bdlwr,bdupr)
            endif
        else if(iprob.eq.3) then
            call norm3(ip,rp,isw,itnum,u,du,ja,a,b,p,d,adu,gm,z)
            if(isw.le.0) then
                do i=1,ndf
                    usv(i)=u(i)
                enddo
                rlsv=rp(21)
                step0=1.0d0
            endif
        else if(iprob.eq.4) then
            call norm4(ip,rp,isw,itnum,u,um,du,dum,ja,a,h,b,p,
     +          d,dl,adu,hdu,gm,z)
            if(isw.le.0) then
                do i=1,ndf
                    usv(i)=u(i)
                    umsv(i)=um(i)
                enddo
                rlsv=rp(21)
                rllwr=rp(4)
                rlupr=rp(5)
                delta=rp(72)
                if(delta.lt.0.0d0) then
                    step0=dmin1((rllwr-rlsv)/delta,1.0d0)
                else if(delta.gt.0.0d0) then
                    step0=dmin1((rlupr-rlsv)/delta,1.0d0)
                else
                    step0=1.0d0
                endif
            endif
        else if(iprob.eq.5) then
            call norm5(ip,rp,isw,itnum,u,um,uc,du,dum,duc,
     +          ja,a,h,g,su,sm,b,p,dl,adu,hdu,gm,z)
            if(isw.le.0) then
                do i=1,ndf
                    usv(i)=u(i)
                    umsv(i)=um(i)
                    ucsv(i)=uc(i)
                enddo
                step0=stepmx(ndf,uc,duc,bdlwr,bdupr)
            endif
        else
            call norm6(ip,rp,isw,itnum,u,du,ja,a,b,adu,gm,z)
            if(isw.le.0) then
                do i=1,ndf
                    usv(i)=u(i)
                enddo
                step0=1.0d0
            endif
        endif
c
c       compute new step
c
        call cstep(rp,0,isw,step0)
        if(isw.eq.-1) return
c
c    update solution with current step
c
        step=rp(52)
        delta=rp(72)
        if(iprob.eq.1.and.itask.eq.9) then
            do i=1,ndf
                um(i)=umsv(i)+step*dum(i)
            enddo
        else if(iprob.eq.3) then
            rp(21)=rlsv+step*delta
        else if(iprob.eq.4) then
            rp(21)=rlsv+step*delta
            do i=1,ndf
                um(i)=umsv(i)+step*dum(i)
            enddo
        else if(iprob.eq.5) then
            do i=1,ndf
                um(i)=umsv(i)+step*dum(i)
                uc(i)=ucsv(i)+step*duc(i)
            enddo
        endif
        do i=1,ndf
            u(i)=usv(i)+step*du(i)
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine cstep(rp,iexsw,isw,step0)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            double precision
     +          rp(100)
            save ksw,ibit,tol,eps,snew,sold,sleft,sright,
     +          dnew,dold,fnew,fold
            data ibit/0/
c
c       this routine carries out a bisection
c       or secant iteration
c
c       isw = 0   initialize
c           > 0   update
c           < 0   converged
c
c       initialization
c
        if(isw.le.0) then
            eps=1.0d2*ceps(ibit)
            tol=1.d-2
            snew=0.0d0
            sleft=0.0d0
            sright=0.0d0
            dnew=rp(58)
            fnew=rp(56)**2/2.0d0
            step=rp(52)
            ratio=rp(57)
            step=step/(step+(1.0d0-step)*ratio/100.0d0)
            if(step0.lt.1.0d0) then
                frac=dmax1(0.75d0,0.98d0-rp(63))
                step=dmin1(step,frac*step0)
            endif
            if(iexsw.eq.1) call exstep(step)
            isw=1
            ksw=0
            rp(52)=step
            return
        endif
c
c       the case isw > 0
c
        isw=isw+1
        sold=snew
        snew=rp(52)
        dold=dnew
        dnew=rp(58)
        fold=fnew
        fnew=rp(56)**2/2.0d0
        relres=rp(56)
        ratio=rp(57)
        relerr=rp(54)
c
        if(sright.le.0.0d0.or.dnew.gt.0.0d0.or.ksw.eq.1) then
            sright=snew
            if(dnew.le.0.0d0) then
                ksw=1
            else
                ksw=0
            endif
        else
            sleft=snew
        endif
c
c       sufficient decrease
c
        ds=sright-sleft
        if(ds.le.tol.and.dnew.le.0.0d0) isw=-1
        if(ratio.le.1.0d0-eps*snew.and.dnew.le.0.0d0) isw=-1
        if(dmin1(relerr,relres).le.eps) isw=-1
        if(isw.eq.-1) return
c
c       bisection step
c
        rp(52)=(sleft+sright)/2.0d0
        if(ksw.eq.0) then
c
c       secant step
c
            if(dold.eq.dnew) return
            step=snew-dnew*(snew-sold)/(dnew-dold)
        else
c
c       cubic interpolation step
c
            ff=-(fold-fnew)*6.0d0/(sold-snew)
            gg=(dold+dnew)
            a=ff+gg*3.0d0
            b=-(ff+2.0d0*(gg+dnew))
            c=dnew
            if(snew.gt.sold) then
                a=-a
                b=-b
                c=-c
            endif
            rr=dmax1(dabs(a),dabs(b),dabs(c))*eps
c
c       quadratic case
c
            if(dabs(a).lt.rr) then
c
c       b > 0 for min
c
                if(b.le.rr) return
                step=snew-(c/b)*(sold-snew)
            else
c
c       cubic case
c
                b=b/(2.0d0*a)
                c=c/a
                discr=b**2-c
                if(discr.le.0.0d0) return
                d=dsqrt(discr)
                if(b.lt.0.0d0) then
c
c       the min occurs for 2*a r + b > 0 (not b/2a above)
c
                    if(a.gt.0.0d0) then
                        r=-(b-d)
                    else
                        r=-c/(b-d)
                    endif
                else
                    if(a.lt.0.0d0) then
                        r=-(b+d)
                    else
                        r=-c/(b+d)
                    endif
                endif
                step=snew+r*(sold-snew)
            endif
        endif
c
c       choose alternative
c
        dl=dabs(step-sleft)
        dr=dabs(step-sright)
        if(dmax1(dl,dr).le.ds*(1.0d0-tol)) then
            rp(52)=step
        else if(dl.le.ds*tol) then
            rp(52)=sleft+ds*tol
        else if(dr.le.ds*tol) then
            rp(52)=sright-ds*tol
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        double precision function stepmx(n,u,du,bdlwr,bdupr)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            double precision
     +          u(*),du(*),bdlwr(*),bdupr(*)
c
c       compute maximum step for interior point
c
        stepmx=1.0d0
        do i=1,n
            if(du(i).lt.0.0d0) then
                stepmx=dmin1((bdlwr(i)-u(i))/du(i),stepmx)
            else if(du(i).gt.0.0d0) then
                stepmx=dmin1((bdupr(i)-u(i))/du(i),stepmx)
            endif
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine norm1(ip,rp,isw,itnum,u,du,um,dum,
     +      ja,a,b,p,adu,adum,gm,z)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100),ja(*)
            double precision
     +          u(*),du(*),a(*),b(*),adu(*),gm(*),z(*),um(*),
     1          dum(*),p(*),rp(100),adum(*)
            save bnorm0,bmnrm0,blast,bmlast,ibit,eps
            data bnorm0,bmnrm0,blast,bmlast,ibit/
     +          0.0d0,0.0d0,0.0d0,0.0d0,0/
c
c       compute norms -- iprob=1
c
        ndf=ip(5)
        ispd=ip(8)
        jspd=1
        if(ispd.ne.1) jspd=-1
c
        call mtxml0(ndf,ja,a,du,adu,z,ispd)
        bnorm=dl2nrm(ndf,b,gm,-1)
        dgamma=dl2ip(ndf,b,adu,gm,-1)
c
        call mtxml0(ndf,ja,a,dum,adum,z,jspd)
        bmnorm=dl2nrm(ndf,p,gm,-1)
        gammam=dl2ip(ndf,p,adum,gm,-1)
c
        if(isw.le.0) then
            eps=1.0d2*ceps(ibit)
c
            enorm=dl2nrm(ndf,du,gm,1)
            unorm=dl2nrm(ndf,u,gm,1)
            relerr=1.0d0
            if(unorm.gt.enorm) relerr=enorm/unorm
            if(unorm+enorm.le.0.0d0) relerr=0.0d0
            emnorm=dl2nrm(ndf,dum,gm,1)
            umnorm=dl2nrm(ndf,um,gm,1)
            relerm=1.0d0
            if(umnorm.gt.emnorm) relerm=emnorm/umnorm
            if(umnorm+emnorm.le.0.0d0) relerm=0.0d0
            rp(54)=relerr+relerm
            rp(54)=relerr
c
            if(bnorm.le.0.0d0) bnorm=eps
            if(bmnorm.le.0.0d0) bmnorm=eps
            if(itnum.eq.1) then
                bnorm0=dmax1(bnorm,rp(59))
                rp(59)=bnorm0
                bmnrm0=dmax1(bmnorm,rp(60))
                rp(60)=bmnrm0
            endif
        else
            rp(56)=bnorm/bnorm0+bmnorm/bmnrm0
            rp(57)=bnorm/blast+bmnorm/bmlast
            rp(56)=bnorm/bnorm0
            rp(57)=bnorm/blast
        endif
c
        ddnew=-dgamma/bnorm0**2
        dmdnew=-gammam/bmnrm0**2
        rp(58)=ddnew+dmdnew
        rp(58)=ddnew
        blast=bnorm
        bmlast=bmnorm
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine norm2(ip,rp,isw,itnum,u,du,ja,a,b,adu,gm,z)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100),ja(*)
            double precision
     +          rp(100),u(*),du(*),a(*),b(*),adu(*),gm(*),z(*)
            save bnorm0,blast,ibit,eps
            data bnorm0,blast,ibit/0.0d0,0.0d0,0/
c
c       compute norms -- iprob=2
c
        ndf=ip(5)
        ispd=ip(8)
c
        call mtxml0(ndf,ja,a,du,adu,z,ispd)
        bnorm=dl2nrm(ndf,b,gm,-1)
        dgamma=dl2ip(ndf,b,adu,gm,-1)
c
        if(isw.eq.0) then
            eps=1.0d2*ceps(ibit)
c
            enorm=dl2nrm(ndf,du,gm,1)
            unorm=dl2nrm(ndf,u,gm,1)
            relerr=1.0d0
            if(unorm.gt.enorm) relerr=enorm/unorm
            if(unorm+enorm.le.0.0d0) relerr=0.0d0
            rp(54)=relerr
c
            if(bnorm.le.0.0d0) bnorm=eps
            if(itnum.eq.1) then
                bnorm0=dmax1(bnorm,rp(59))
                rp(59)=bnorm0
            endif
        else
            rp(56)=bnorm/bnorm0
            rp(57)=bnorm/blast
        endif
c
        ddnew=-dgamma/bnorm0**2
        rp(58)=ddnew
        blast=bnorm
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine norm3(ip,rp,isw,itnum,u,du,ja,a,b,p,d,adu,gm,z)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100),ja(*)
            double precision
     +          rp(100),u(*),du(*),a(*),b(*),p(*),d(*),
     1          adu(*),gm(*),z(*)
            save bnorm0,blast,ibit,eps
            data bnorm0,blast,ibit/0.0d0,0.0d0,0/
c
c       compute norms -- iprob=3
c
        ndf=ip(5)
        ispd=ip(8)
        rl=rp(21)
        scale=dsqrt(rp(68))
        scleqn=rp(67)*scale
        thetal=rp(69)*scale
        thetar=rp(70)*scale
        delta=rp(72)
        drdrl=rp(73)
c
c    compute adu
c
        call mtxml0(ndf,ja,a,du,adu,z,ispd)
        ss=thetar*(rl2ip(ndf,p,du)+drdrl*delta)+thetal*delta
        bnorm=dsqrt(dl2nrm(ndf,b,gm,-1)**2+scleqn**2)
        dgamma=dl2ip(ndf,b,adu,gm,-1)
        bd=dl2ip(ndf,b,d,gm,-1)
c
        if(isw.le.0) then
            eps=1.0d2*ceps(ibit)
c
c       compute relerr
c
            enorm=dl2nrm(ndf,du,gm,1)
            unorm=dl2nrm(ndf,u,gm,1)
            relerr=1.0d0
            if(unorm.gt.enorm) relerr=enorm/unorm
            if(unorm+enorm.le.0.0d0) relerr=0.0d0
            rlerr=1.0d0
            if(dabs(rl).gt.dabs(delta)) rlerr=dabs(delta)/dabs(rl)
            if(dabs(rl)+dabs(delta).eq.0.0d0) rlerr=0.0d0
            rp(54)=relerr+rlerr
c
            if(bnorm.le.0.0d0) bnorm=eps
            if(itnum.eq.1) then
                bnorm0=dmax1(bnorm,rp(59))
                rp(59)=bnorm0
            endif
        else
            rp(56)=bnorm/bnorm0
            rp(57)=bnorm/blast
        endif
c
        ddnew=(-dgamma+ss*scleqn+bd*delta)/bnorm0**2
        rp(58)=ddnew
        blast=bnorm
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 10.0 - - - september, 2007
c
c-----------------------------------------------------------------------
        subroutine norm4(ip,rp,isw,itnum,u,um,du,dum,ja,a,h,b,p,
     +      d,dl,adu,hdu,gm,z)
c
            implicit double precision (a-h,o-z)
            implicit integer (i-n)
            integer
     +          ip(100),ja(*)
            double precision
     +          rp(100),u(*),um(*),du(*),dum(*),a(*),h(*),b(*),p(*),
     1          d(*),dl(*),adu(*),hdu(*),gm(*),z(*)
            save bnorm0,blast,ibit,eps
            data bnorm0,blast,ibit/0.0d0,0.0d0,0/
c
c       compute norms -- iprob=4
c
        ndf=ip(5)
        ispd=ip(8)
        jspd=1
        if(ispd.ne.1) jspd=-1
        scleqn=rp(67)
        seqdot=rp(74)
        delta=rp(72)
        rl=rp(21)
c
c       matrix multiplies
c
        call mtxml0(ndf,ja,h,du,hdu,z,1)
        call mtxml0(ndf,ja,a,dum,adu,z,jspd)
        do i=1,ndf
            hdu(i)=hdu(i)+adu(i)-delta*dl(i)
        enddo
        call mtxml0(ndf,ja,a,du,adu,z,ispd)
        do i=1,ndf
            adu(i)=adu(i)-delta*d(i)
        enddo
        bnorm=dl2nrm(ndf,b,gm,-1)
        dgamma=dl2ip(ndf,b,adu,gm,-1) 
        pnorm=dl2nrm(ndf,p,gm,-1)
        pgamma=dl2ip(ndf,p,hdu,gm,-1)
        bnorm=dsqrt(scleqn**2+bnorm**2+pnorm**2)
        c=-rl2ip(ndf,du,dl)-rl2ip(ndf,dum,d)-seqdot*delta
c
        if(isw.le.0) then
            eps=1.0d2*ceps(ibit)
c
            uunorm=dl2nrm(ndf,u,gm,1)
            umnorm=dl2nrm(ndf,um,gm,1)
            eunorm=dl2nrm(ndf,du,gm,1)
            emnorm=dl2nrm(ndf,dum,gm,1)
c
c       compute relerr
c
            rulerr=1.0d0
            if(uunorm.gt.eunorm) rulerr=eunorm/uunorm
            if(uunorm+eunorm.le.0.0d0) rulerr=0.0d0
            rmlerr=1.0d0
            if(umnorm.gt.emnorm) rmlerr=emnorm/umnorm
            if(umnorm+emnorm.le.0.0d0) rmlerr=0.0d0
            rlerr=1.0d0
            if(dabs(rl).gt.dabs(delta)) rlerr=dabs(delta)/dabs(rl)
            if(dabs(rl)+dabs(delta).eq.0.0d0) rlerr=0.0d0
            rp(54)=rulerr+rmlerr+rlerr
c
            if(bnorm.le.0.0d0) bnorm=eps
            if(it