dnsid.f

Go to the documentation of this file.

C Work performed under the auspices of the U.S. Department of Energy
C by Lawrence Livermore National Laboratory under contract number 
C W-7405-Eng-48.
C
      SUBROUTINE dnsid(X,Y,YPRIME,NEQ,ICOPT,ID,RES,WT,RPAR,IPAR,
     *   delta,r,yic,ypic,wm,iwm,cj,epcon,ratemx,maxit,stptol,
     *   icnflg,icnstr,iernew)
C
C***BEGIN PROLOGUE  DNSID
C***REFER TO  DDASPK
C***DATE WRITTEN   940701   (YYMMDD)
C***REVISION DATE  950713   (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C     DNSID solves a nonlinear system of algebraic equations of the
C     form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME
C     in the initial conditions.
C
C     The method used is a modified Newton scheme.
C
C     The parameters represent
C
C     X         -- Independent variable.
C     Y         -- Solution vector.
C     YPRIME    -- Derivative of solution vector.
C     NEQ       -- Number of unknowns.
C     ICOPT     -- Initial condition option chosen (1 or 2).
C     ID        -- Array of dimension NEQ, which must be initialized
C                  if ICOPT = 1.  See DDASIC.
C     RES       -- External user-supplied subroutine to evaluate the
C                  residual.  See RES description in DDASPK prologue.
C     WT        -- Vector of weights for error criterion.
C     RPAR,IPAR -- Real and integer arrays used for communication
C                  between the calling program and external user
C                  routines.  They are not altered within DASPK.
C     DELTA     -- Residual vector on entry, and work vector of
C                  length NEQ for DNSID.
C     WM,IWM    -- Real and integer arrays storing matrix information
C                  such as the matrix of partial derivatives,
C                  permutation vector, and various other information.
C     CJ        -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).
C     R         -- Array of length NEQ used as workspace by the 
C                  linesearch routine DLINSD.
C     YIC,YPIC  -- Work vectors for DLINSD, each of length NEQ.
C     EPCON     -- Tolerance to test for convergence of the Newton
C                  iteration.
C     RATEMX    -- Maximum convergence rate for which Newton iteration
C                  is considered converging.
C     MAXIT     -- Maximum allowed number of Newton iterations.
C     STPTOL    -- Tolerance used in calculating the minimum lambda
C                  value allowed.
C     ICNFLG    -- Integer scalar.  If nonzero, then constraint
C                  violations in the proposed new approximate solution
C                  will be checked for, and the maximum step length 
C                  will be adjusted accordingly.
C     ICNSTR    -- Integer array of length NEQ containing flags for
C                  checking constraints.
C     IERNEW    -- Error flag for Newton iteration.
C                   0  ==> Newton iteration converged.
C                   1  ==> failed to converge, but RATE .le. RATEMX.
C                   2  ==> failed to converge, RATE .gt. RATEMX.
C                   3  ==> other recoverable error (IRES = -1, or
C                          linesearch failed).
C                  -1  ==> unrecoverable error (IRES = -2).
C
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C   DSLVD, DDWNRM, DLINSD, DCOPY
C
C***END PROLOGUE  DNSID
C
C
      IMPLICIT DOUBLE PRECISION(a-h,o-z)
      dimension y(*),yprime(*),wt(*),r(*)
      dimension id(*),delta(*), yic(*), ypic(*)
      dimension wm(*),iwm(*), rpar(*),ipar(*)
      dimension icnstr(*)
      EXTERNAL  res
C
      parameter(lnni=19, llsoff=35)
C
C
C     Initializations.  M is the Newton iteration counter.
C
      lsoff = iwm(llsoff)
      m = 0
      rate = 1.0d0
      rlx = 0.4d0
C
C     Compute a new step vector DELTA by back-substitution.
C
      CALL dslvd(neq, delta, wm, iwm)
C
C     Get norm of DELTA.  Return now if norm(DELTA) .le. EPCON.
C
      delnrm = ddwnrm(neq,delta,wt,rpar,ipar)
      fnrm = delnrm
      IF (fnrm .LE. epcon) RETURN
C
C     Newton iteration loop.
C
 300  CONTINUE
      iwm(lnni) = iwm(lnni) + 1
C
C     Call linesearch routine for global strategy and set RATE
C
      oldfnm = fnrm
C
      CALL dlinsd(neq, y, x, yprime, cj, delta, delnrm, wt, lsoff,
     *             stptol, iret, res, ires, wm, iwm, fnrm, icopt, id,
     *             r, yic, ypic, icnflg, icnstr, rlx, rpar, ipar)
C
      rate = fnrm/oldfnm
C
C     Check for error condition from linesearch.
      IF (iret .NE. 0) go to 390
C
C     Test for convergence of the iteration, and return or loop.
C
      IF (fnrm .LE. epcon) RETURN
C
C     The iteration has not yet converged.  Update M.
C     Test whether the maximum number of iterations have been tried.
C
      m = m + 1
      IF (m .GE. maxit) go to 380
C
C     Copy the residual to DELTA and its norm to DELNRM, and loop for
C     another iteration.
C
      CALL dcopy(neq, r, 1, delta, 1)
      delnrm = fnrm      
      go to 300
C
C     The maximum number of iterations was done.  Set IERNEW and return.
C
 380  IF (rate .LE. ratemx) THEN
         iernew = 1
      ELSE
         iernew = 2
      ENDIF
      RETURN
C
 390  IF (ires .LE. -2) THEN
         iernew = -1
      ELSE
         iernew = 3
      ENDIF
      RETURN
C
C
C------END OF SUBROUTINE DNSID------------------------------------------
      END