anchor.diff:          if(numslp.ne.0) call daxpy(ndof*numnp,one,dxcur,ione,dx,ione)
           if(numslp.ne.0) call daxpy(ndof*numnp,one,dxcur,ione,deldx,ione)
   c
  -c...integrate the stresses and compute the equivalent nodal loads
  -c
  -cdebug        write(6,*) "From iterate right before stresn call:"
  -cdebug        write(6,*) "gauss:",(gauss(idb),idb=1,(nsd+1)*ngauss)
  -        call stresn(x,b,d,dx,tfault,stn,deps,beta,betb,scur,st0,dbeta,
  -     &   dbetb,skew,ien,lm,lmx,lmf,dmat,mat,prop,histry,infin,gauss,
  -     &   rtimdat,stol,iddmat,nen,numel,ndof,nsd,numnp,neq,nee,
  -     &   nstr,ngauss,nppts,ngem,nskdim,nhist,nprop,numat,numfn,numslp,
  -     &   numrot,lastep,nstep,lgdefp,ibbarp,ivisc,iplas,imhist,ipstrs,
  -     &   nprestr,nddmat,ndmat,idebug,idout,kto,kw,fulout)
  +c...  integrate the stresses and compute the equivalent nodal loads,
  +c     updating the material and stiffness matrices if requested.
  +c
  +        if(reform) then
  +          call stress_mat_drv(
  +     &     alnz,ja,nnz,                                                 ! sparse
  +     &     b,neq,                                                       ! force
  +     &     x,d,iwink,wink,numnp,nwink,                                  ! global
  +     &     dx,iwinkx,winkx,numslp,numsn,nwinkx,                         ! slip
  +     &     tfault,numfn,                                                ! fault
  +     &     s,stemp,                                                     ! stiff
  +     &     state,dstate,dmat,ien,lm,lmx,lmf,infiel,iddmat,nstatesz,     ! elemnt
  +     &     ndmatsz,numelt,nconsz,                                       ! elemnt
  +     &     prop,mhist,infmat,numat,npropsz,                             ! materl
  +     &     gauss,sh,shj,infetype,                                       ! eltype
  +     &     histry,rtimdat,rgiter,ntimdat,nhist,lastep,stress_mat_cmp,   ! timdat
  +     &     skew,numrot,                                                 ! skew
  +     &     getshape,bmatrix,                                            ! bbar
  +     &     ierr,errstrng)                                               ! errcode
  +        else
  +          call stress_drv(
  +     &     b,neq,                                                       ! force
  +     &     x,d,numnp,                                                   ! global
  +     &     dx,numslp,                                                   ! slip
  +     &     tfault,numfn,                                                ! fault
  +     &     state,dstate,dmat,ien,lm,lmx,lmf,infiel,nstatesz,ndmatsz,    ! elemnt
  +     &     numelt,nconsz,                                               ! elemnt
  +     &     prop,mhist,infmat,numat,npropsz,                             ! materl
  +     &     gauss,sh,shj,infetype,                                       ! eltype
  +     &     histry,rtimdat,rgiter,ntimdat,nhist,lastep,stress_cmp,       ! timdat
  +     &     skew,numrot,                                                 ! skew
  +     &     getshape,bmatrix,                                            ! bbar
  +     &     ierr,errstrng)                                               ! errcode
  +        end if
   c
   c...compute the out-of-balance forces and convergence criteria
   c
           call residu(b,bres,btot,gvec2,gtol,gi,gprev,gcurr,id,idx,neq,
  -     &   ndof,numnp,i,itmaxp,idebug,idout,kto,kw,converge)
  +     &   numnp,iter,itmaxp,idebug,idout,kto,kw,converge)
   c
   c...if solution has converged, set equilibrium stresses and creep
   c   strains to their current values
   c
           if(converge) then
  -          call dcopy(nstr*numel*ngauss,scur,ione,stn,ione)
  -          if(ivisc.eq.1) call daxpy(nstr*numel*ngauss,one,dbeta,ione,
  -     &     beta,ione)
  -          if(iplas.eq.1) call daxpy(nstr*numel*ngauss,one,dbetb,ione,
  -     &     betb,ione)
  -          call daxpy(nstr*numel*ngauss,-one,eps,ione,deps,ione)
  -          call daxpy(nstr*numel*ngauss,one,deps,ione,eps,ione)
  +          call update_state(state,dstate,infiel,infmat,infetype,
  +     &     nstatesz,numelt,numat,ierr,errstrng)
             return
           end if
           reform=.false.
  -c*        call flush(kto)
  -c*        call flush(kw)
  -c*        call flush(kp)
         end do
         return
         end
   c
   c version
  -c $Id: iterate.f,v 1.1 2004/04/14 21:18:30 willic3 Exp $
  +c $Id: iterate.f,v 1.2 2004/07/05 20:37:17 willic3 Exp $
   c
   c Generated automatically by Fortran77Mill on Wed May 21 14:15:03 2003
   c
  
positive.code: c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  c
  c
        subroutine iterate(alnz,pcg,zcg,ja,                               ! sparse
       & b,btot,bres,pvec,gvec1,gvec2,                                    ! force
       & x,d,dx,deld,deldx,dprev,dcur,dxcur,id,idx,skew,histry,           ! global
       & ien,infin,mat,lm,lmx,lmf,prop,gauss,                             ! elemnt
       & dmat,stn,scur,st0,eps,deps,beta,dbeta,betb,dbetb,iddmat,         ! stress
       & ielno,iside,ihistry,pres,pdir,                                   ! tractn
       & nfault,dfault,tfault,                                            ! split
       & idslp,ipslp,                                                     ! slip
       & iwink,wink,iwinkx,winkx,                                         ! wink
       & s,stemp,                                                         ! local
       & gcurr,gi,gprev,grav,gtol,ncodat,ndimens,                         ! info
       & npar,nprint,nsiter,nsysdat,ntimdat,nunits,nvisdat,rgiter,        ! info
       & rmin,rmult,rtimdat)                                              ! info
  c
  c...subroutine to loop over iterations within each time step
  c
  
negative.code: c
          call fill(dcur,zero,ndof*numnp)
          call fill(dxcur,zero,ndof*numnp)
          call disp(gvec2,dcur,id,numnp,neq)
          if(numslp.ne.0) call disp(gvec2,dxcur,idx,numnp,neq)
  c
  c...rotate skewed coordinates to global system
  c
          if(numrot.ne.0) then
            call rdisp(dcur,skew,numnp)
            if(numslp.ne.0) call rdisp(dxcur,skew,numnp)
          end if
  c
  c...compute contribution to global load vector from winkler boundary
  