      subroutine sdas1 (x, y, yprime, neq, ldd, sdasf, info, h, wt,
     &                   idid, phi, delta, e, wm, iwork, rwork)
c Copyright (c) 2006, Math a la Carte, Inc.
c>> 2009-09-30 sdas1 Krogh  Removed unused variable sc.
c>> 2008-08-26 sdas1 Hanson add argument of leading dimension to sdasf
c>> 2006-05-18 sdas1 Hanson added checks for inconsistent constraints
c>> 2006-04-16 sdas1 Krogh  declared j.
c>> 2004-04-13 sdas1 Hanson cleared phi array at outset.
c>> 2003-03-06 sdas1 Hanson changed weights to use reciprocals
c>> 2001-12-12 sdas1 Krogh  Changed code for reverse communication
c>> 2001-11-23 sdas1 Krogh  Changed many names per library conventions.
c>> 2001-11-01 sdas1 Hanson Provide code to Math a la Carte.
c--S replaces "?": ?das1, ?dasf, ?dasj, ?daslv, ?dasco,
c--&  ?daslx, ?dasnm, ?dasdb
c***BEGIN PROLOGUE  SDAS1
c***SUBSIDIARY
c***PURPOSE  Initialization routine for SDASLX.
c***LIBRARY   SLATEC (SDASLX)
c***TYPE      DOUBLE PRECISION (SDAS1-S, SDAS1-D)
c***AUTHOR  Petzold, Linda R., (LLNL)
c***DESCRIPTION
c-----------------------------------------------------------------
c     SDAS1 TAKES ONE STEP OF SIZE H OR SMALLER WITH THE BACKWARD EULER
c     METHOD, TO FIND YPRIME.  X AND Y ARE UPDATED TO BE CONSISTENT WITH
c     THE NEW STEP.  A MODIFIED DAMPED NEWTON ITERATION IS USED TO SOLVE
c     THE CORRECTOR ITERATION.
c
c     THE INITIAL GUESS FOR YPRIME IS USED IN THE PREDICTION, AND IN
c     FORMING THE ITERATION MATRIX, BUT IS NOT INVOLVED IN THE ERROR
c     TEST. THIS MAY HAVE TROUBLE CONVERGING IF THE INITIAL GUESS IS NO
c     GOOD, OR IF G(X,Y,YPRIME) DEPENDS NONLINEARLY ON YPRIME.
c
c     THE PARAMETERS REPRESENT:
c     X --         INDEPENDENT VARIABLE
c     Y --         SOLUTION VECTOR AT X
c     YPRIME --    DERIVATIVE OF SOLUTION VECTOR
c     NEQ --       NUMBER OF EQUATIONS
c     LDD --       FIRST (ROW) DIMENSION FOR MATIX.
c     SDASF --     THE USER DEFINED ROUTINE.
c     INFO --      THE USER SUPPLIED ARRAY DEFINING OPTIONS.
c     H --         STEPSIZE. IMDER MAY USE A STEPSIZE SMALLER THAN H.
c     WT --        VECTOR OF WEIGHTS FOR ERROR CRITERION
c     IDID --      COMPLETION CODE WITH THE FOLLOWING MEANINGS
c       IDID= 1 -- YPRIME WAS FOUND SUCCESSFULLY
c       IDID=-10 -- SDAS1 FAILED TO FIND YPRIME
c       IDID=-28 -- IRES NOT RESET WHEN IRES=5 WITH USER DEFINED MATRIX.
c     PHI --       WORK SPACE FOR SDAS1
c     DELTA,E --   WORK SPACE FOR SDAS1
c     WM,IWORK --  REAL AND INTEGER ARRAYS STORING
c                  MATRIX INFORMATION
c     RWORK --     THE USUAL WORK ARRAY.
c
c-----------------------------------------------------------------
c***ROUTINES CALLED  SDASF, SDASNM, SDASLV
c***REVISION HISTORY  (YYMMDD)
c   830315  DATE WRITTEN
c   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
c   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
c   901026  Added explicit declarations for all variables and minor
c           cosmetic changes to prologue.  (FNF)
c   901030  Minor corrections to declarations.  (FNF)
c   981118  Use one external user routine, RJH.
c***END PROLOGUE  SDAS1
c
      integer    neq, ldd, idid, iwork(*), info(16)
      integer    locate
      real             x, y(*), yprime(*), h, wt(*), phi(neq,*),
     &     delta(*), e(*), wm(*), rwork(*)
      external   sdasf
c
      external   sdasj, sdasnm, sdaslv
      real             sdasnm
c
      integer    i, ires, j, jcalc, m, maxit, mjac, ncf,
     &           nef, nsf
      real             damp, delnrm, err, oldnrm, r, rate, s,
     &           xold, ynorm
      logical    convgd
c
c     POINTERS INTO IWORK
      integer  lml, lmu, lires, ldelt, lwm, lmxord, lk, lkold, lmat,
     &     lcnstr, lns, lnstl, lnst, lnre, lnje, letf, lctf, lnpd,
     &     ljcalc, lphase, revloc, mxstep, le, lwt, lphi, ntemp, lipvt
      parameter (lml=1, lmu=2, lires=3, ldelt=4, lwm=5, lmxord=6, lk=7,
     &     lkold=8, lmat=9, lcnstr=10, lns=11, lnstl=12, lnst=13,
     &     lnre=14, lnje=15, letf=16, lctf=17, lnpd=18, ljcalc=19,
     &     lphase=20, revloc=21, mxstep=22, le=23, lwt=24, lphi=25,
     &     ntemp=26, lipvt=31)
c
c     POINTERS INTO RWORK
      integer  lcj, ltstop, lhmax, lh, ltn, lcjold, lhold, lnjac,
     &     lround, lhmin, lalpha, lbeta, lgamma, lpsi, lsigma, ldelta
      parameter  (lcj=1, ltstop=2, lhmax=3, lh=4, ltn=5, lcjold=6,
     &     lhold=7, lnjac=8, lround=9, lhmin=10, lalpha=11, lbeta=17,
     &     lgamma=23, lpsi=29, lsigma=35, ldelta=46)
c
c     POINTERS INTO INFO
      integer  itol, iout, istop, imat, idb, imaxh, ih0, iord, icnstr,
     &     inityp, ixstep
      parameter  (itol=2, iout=3,  istop=4, imat=5, idb=6, imaxh=7,
     &     ih0=8, iord=9, icnstr=10, inityp=11, ixstep= 12)
 
      data maxit/10/, mjac/5/
      data damp/0.75e0/
c
      save
c
c---------------------------------------------------
c     BLOCK 1.
c     INITIALIZATIONS.
c---------------------------------------------------
c
c***FIRST EXECUTABLE STATEMENT  SDAS1
      if (iwork(revloc) .ne. 0) then
         ires = iwork(lires)
         locate = mod(iwork(revloc), 8)
         iwork(revloc) = iwork(revloc) / 8
c                1   2    3    4   5    6
         go to (50, 70, 100, 320, 75, 360), locate
      end if
c     Control drops through here on first call:
      idid   = 1
      nef    = 0
      ncf    = 0
      nsf    = 0
      xold   = x
      ynorm  = sdasnm(neq,y,wt,rwork,iwork)
c
c     SAVE Y AND YPRIME IN PHI
      do 10  i=1, neq
         phi(i,1) = y(i)
         phi(i,2) = yprime(i)
   10 continue
c
c  Clear rest of phi(*,*) array to define differences.
      do 15 j=3, iwork(lmxord)+1
        do 15 i=1,neq
           phi(i,j) = 0.e0
   15 continue
 
c
c
c----------------------------------------------------
c     BLOCK 2.
c     DO ONE BACKWARD EULER STEP.
c----------------------------------------------------
c
c     SET UP FOR START OF CORRECTOR ITERATION
   20 rwork(lcj) = 1.0e0/h
      x = x + h
c
c     PREDICT SOLUTION AND DERIVATIVE
      do 30  i=1, neq
         y(i) = y(i) + h*yprime(i)
   30 continue
c
      jcalc  = -1
      m      = 0
      convgd = .true.
c
c
c     CORRECTOR LOOP.
   40 iwork(lnre) = iwork(lnre) + 1
      ires = 1
      if (info(idb) .ne. 0) call sdasdb(2, neq, x, y, yprime, info,
     &     rwork, iwork, ires,rwork, rwork)
      if (iwork(lmat) .ge. 0) then
         call sdasf (x, y, yprime, delta, wm, LDD, rwork(lcj), ires,
     &   rwork, iwork)
      else
         iwork(revloc) = 8 * iwork(revloc) + 1
         iwork(lires) = ires
         return
      end if
c REVERSE ENTRY 1:
   50 continue
c     Test signal for invalid request.
      if (info(idb) .ne. 0) call sdasdb(3, neq, x, y, yprime,
     &    info, rwork, iwork, ires, rwork, rwork)
      if (ires .lt. 0) then
         if (ires .eq. -1) go to 190
         if (ires .eq. -2) go to 240
         info(idb) = -ires
         ires = 0
      end if
c
c
c     EVALUATE THE ITERATION MATRIX
      if (jcalc .ne. -1) go to 80
      iwork(lnje) = iwork(lnje) + 1
      jcalc     = 0
c REVERSE ENTRY 2:
   70 continue
      call sdasj (neq, ldd, x, y, yprime, delta, h, wt, e, wm, iwork,
     &             rwork, sdasf, info, ires)
 
c     See if reverse communication needed:
      if (iwork(revloc) .ne. 0) then
         if (iwork(revloc) .lt. 0) then
            iwork(revloc) = 5
         else
            iwork(revloc) = 8 * iwork(revloc) + 2
         end if
         return
      end if
c
c REVERSE ENTRY 5:  After user computes and factors the matrix.
   75 continue
      s = 1000000.e0
      if (ires .lt. 0) then
         if (ires .eq. -1) go to 190
         if (ires .eq. -2) go to 240
         info(idb) = -ires
         ires = 0
      end if
      if (ires .ne. 0) go to 190
      nsf = 0
c
c
c
c     MULTIPLY RESIDUAL BY DAMPING FACTOR
   80 continue
      do 90  i=1, neq
         delta(i) = delta(i)*damp
   90 continue
c
c     COMPUTE A NEW ITERATE (BACK SUBSTITUTION)
c     STORE THE CORRECTION IN DELTA
c
      call sdaslv (neq, ldd, x, y, yprime, delta, sdasf, info,
     &             iwork, rwork)
      if (iwork(revloc) .ne. 0) then
         iwork(revloc) = 3
         return
      end if
c REVERSE ENTRY 3:
  100 continue
c
c     UPDATE Y AND YPRIME
      do 110  i=1, neq
         y(i) = y(i) - delta(i)
         yprime(i) = yprime(i) - rwork(lcj)*delta(i)
  110 continue
 
c
c     TEST FOR CONVERGENCE OF THE ITERATION.
c
      delnrm = sdasnm(neq,delta,wt,rwork,iwork)
      if (delnrm .le. 100.e0 * rwork(lround) * ynorm) go to 200
c
      if (m .gt. 0) go to 120
      oldnrm = delnrm
      go to 130
c
  120 rate = (delnrm/oldnrm)**(1.0e0/m)
      if (rate .gt. 0.90e0) go to 190
      s = rate/(1.0e0-rate)
c
  130 if (s*delnrm .le. 0.33e0) go to 200
c
c
c     THE CORRECTOR HAS NOT YET CONVERGED. UPDATE
c     M AND AND TEST WHETHER THE MAXIMUM
c     NUMBER OF ITERATIONS HAVE BEEN TRIED.
c     EVERY MJAC ITERATIONS, GET A NEW
c     ITERATION MATRIX.
c
      m = m + 1
      if (m .ge. maxit) go to 190
c
      if ((m/mjac)*mjac .eq. m) jcalc = -1
 
      go to 40
c
c
c     THE ITERATION HAS CONVERGED.
 
c     EXITS FROM CORRECTOR LOOP.
  190 convgd = .false.
      go to 220
  200 continue
c
c
c
c-----------------------------------------------------
c     BLOCK 3.
c     THE CORRECTOR ITERATION CONVERGED.
c     DO ERROR TEST.
c-----------------------------------------------------
c
      do 210  i=1, neq
         e(i) = y(i) - phi(i,1)
  210 continue
      err = sdasnm(neq,e,wt,rwork,iwork)
c
      if (err .le. 1.0e0) go to 280
c
c
c
c--------------------------------------------------------
c     BLOCK 4.
c     THE BACKWARD EULER STEP FAILED. RESTORE X, Y
c     AND YPRIME TO THEIR ORIGINAL VALUES.
c     REDUCE STEPSIZE AND TRY AGAIN, IF
c     POSSIBLE.
c---------------------------------------------------------
c
  220 continue
      x = xold
      do 230  i=1, neq
         y(i) = phi(i,1)
         yprime(i) = phi(i,2)
  230 continue
c
      if (convgd) go to 260
      if (ires .ne. 0) then
         nsf = nsf + 1
         h   = h*0.25e0
         if ((nsf.lt.3) .and. (abs(h) .ge. rwork(lhmin))) go to 270
         idid = -10
         go to 280
      end if
      if (ires .ne. -2) go to 250
  240 idid = -10
      go to 280
  250 ncf = ncf + 1
      h = h*0.25e0
      if ((ncf.lt.10) .and. (abs(h) .ge. rwork(lhmin))) go to 270
      idid = -10
      go to 280
c
  260 nef = nef + 1
      r = 0.90e0/(2.0e0*err+0.0001e0)
      r = max(0.1e0,min(0.5e0,r))
      h = h*r
      if ((abs(h).ge.rwork(lhmin)) .and. (nef.lt.10)) go to 270
      idid = -10
      go to 280
  270 go to 20
c
  280 continue
c     CHECK for CONSTRAINTS
      if (iwork(lcnstr) .eq. 0 ) go to 400
      ires = 5
 
      if (iwork(lcnstr) .ge. 0) then
         if (info(idb) .ne. 0) call sdasdb(2, neq, x, y, yprime,
     &    info, rwork, iwork, ires, rwork, rwork)
         call sdasf (x, y, yprime, delta, wm, LDD, rwork(lcj),
     &   ires, rwork, iwork)
      else
         iwork(revloc) = 4
         iwork(lires) = ires
         return
      end if
c REVERSE ENTRY 4:
  320 continue
      if (info(idb) .ne. 0) call sdasdb(3, neq, x, y, yprime,
     &    info, rwork, iwork, ires, rwork, rwork)
      if (ires .ne. 5) then
         if (ires .lt. 0) then
            if (ires .eq. -1) go to 190
            if (ires .eq. -2) go to 240
            info(idb) = -ires
            ires = 0
         end if
      else
         i = abs(iwork(lmat))
         if ((i .eq. 5) .or. (i .eq. 6)) then
            idid = -28
            go to 400
         end if
c     If IRES=5, compute a weighted least-distance Newton step
c     back to the constraints. Put move into delta(:).
         call sdasco (wm(neq+1), ldd, neq, wt, delta)
 
      end if
C Check that the moves back to the constraints are smaller than
c sizes of the error tolerances.  If they are then apply the
c moves.  Otherwise evaluate the constraints and return an error flag.
      do 330 i=1, neq
C Recall that tol_i = 1/ wt_i.
        if(abs(delta(i)*wt(i)) .gt. 1) go to 350
  330 continue
C If moving onto the constraints was reasonable, make the move.
      do 340  i=1, neq
        y(i) = y(i) - delta(i)
        yprime(i) = yprime(i) - rwork(lcj)*delta(i)
 340  continue
      go to 400
 350  continue
C Here the projection back to the constraints was
C .gt. the requested tolerance in some component.
c  Evaluate the constraints at the current values.
      if (iwork(lmat) .ge. 0) then
         call sdasf (x, y, yprime, delta, wm, LDD, rwork(lcj), ires,
     &   rwork, iwork)
      else
         iwork(revloc) = 6
         iwork(lires) = ires
         return
      end if
c REVERSE ENTRY 6:
 360  continue
C This error flag will result in printing a message
c and residuals on the constraints.
      IDID = -29
c-------------END OF SUBROUTINE SDAS1----------------------
 400  return
      end
