dlinsd.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 dlinsd (NEQ, Y, T, YPRIME, CJ, P, PNRM, WT, LSOFF,
     *                   stptol, iret, res, ires, wm, iwm,
     *                   fnrm, icopt, id, r, ynew, ypnew, icnflg,
     *                   icnstr, rlx, rpar, ipar)
C
C***BEGIN PROLOGUE  DLINSD
C***REFER TO  DNSID
C***DATE WRITTEN   941025   (YYMMDD)
C***REVISION DATE  941215   (YYMMDD)
C***REVISION DATE  960129   Moved line RL = ONE to top block.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C     DLINSD uses a linesearch algorithm to calculate a new (Y,YPRIME)
C     pair (YNEW,YPNEW) such that 
C
C     f(YNEW,YPNEW) .le. (1 - 2*ALPHA*RL)*f(Y,YPRIME) ,
C
C     where 0 < RL <= 1.  Here, f(y,y') is defined as
C
C      f(y,y') = (1/2)*norm( (J-inverse)*G(t,y,y') )**2 ,
C
C     where norm() is the weighted RMS vector norm, G is the DAE
C     system residual function, and J is the system iteration matrix
C     (Jacobian).
C
C     In addition to the parameters defined elsewhere, we have
C
C     P       -- Approximate Newton step used in backtracking.
C     PNRM    -- Weighted RMS norm of P.
C     LSOFF   -- Flag showing whether the linesearch algorithm is
C                to be invoked.  0 means do the linesearch, and
C                1 means turn off linesearch.
C     STPTOL  -- Tolerance used in calculating the minimum lambda
C                value allowed.
C     ICNFLG  -- Integer scalar.  If nonzero, then constraint violations
C                in the proposed new approximate solution will be
C                checked for, and the maximum step length will be
C                adjusted accordingly.
C     ICNSTR  -- Integer array of length NEQ containing flags for
C                checking constraints.
C     RLX     -- Real scalar restricting update size in DCNSTR.
C     YNEW    -- Array of length NEQ used to hold the new Y in
C                performing the linesearch.
C     YPNEW   -- Array of length NEQ used to hold the new YPRIME in
C                performing the linesearch.
C     Y       -- Array of length NEQ containing the new Y (i.e.,=YNEW).
C     YPRIME  -- Array of length NEQ containing the new YPRIME 
C                (i.e.,=YPNEW).
C     FNRM    -- Real scalar containing SQRT(2*f(Y,YPRIME)) for the
C                current (Y,YPRIME) on input and output.
C     R       -- Work array of length NEQ, containing the scaled 
C                residual (J-inverse)*G(t,y,y') on return.
C     IRET    -- Return flag.
C                IRET=0 means that a satisfactory (Y,YPRIME) was found.
C                IRET=1 means that the routine failed to find a new
C                       (Y,YPRIME) that was sufficiently distinct from
C                       the current (Y,YPRIME) pair.
C                IRET=2 means IRES .ne. 0 from RES.
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C   DFNRMD, DYYPNW, DCOPY
C
C***END PROLOGUE  DLINSD
C
      IMPLICIT DOUBLE PRECISION(a-h,o-z)
      EXTERNAL  res
      dimension y(*), yprime(*), wt(*), r(*), id(*)
      dimension wm(*), iwm(*)
      dimension ynew(*), ypnew(*), p(*), icnstr(*)
      dimension rpar(*), ipar(*)
      CHARACTER msg*80
C
      parameter(lnre=12, lkprin=31)
C
      SAVE alpha, one, two
      DATA alpha/1.0d-4/, one/1.0d0/, two/2.0d0/
C
      kprin=iwm(lkprin)
C
      f1nrm = (fnrm*fnrm)/two
      ratio = one
      IF (kprin .GE. 2) THEN
        msg = '------ IN ROUTINE DLINSD-- PNRM = (R1) )'
        CALL xerrwd(msg, 40, 901, 0, 0, 0, 0, 1, pnrm, 0.0d0)
        ENDIF
      tau = pnrm
      ivio = 0
      rl = one
C-----------------------------------------------------------------------
C Check for violations of the constraints, if any are imposed.
C If any violations are found, the step vector P is rescaled, and the 
C constraint check is repeated, until no violations are found.
C-----------------------------------------------------------------------
      IF (icnflg .NE. 0) THEN
 10      CONTINUE
         CALL dyypnw(neq,y,yprime,cj,rl,p,icopt,id,ynew,ypnew)
         CALL dcnstr(neq, y, ynew, icnstr, tau, rlx, iret, ivar)
         IF (iret .EQ. 1) THEN
            ivio = 1
            ratio1 = tau/pnrm
            ratio = ratio*ratio1
            DO 20 i = 1,neq
 20           p(i) = p(i)*ratio1
            pnrm = tau
            IF (kprin .GE. 2) THEN
              msg = '------ CONSTRAINT VIOL., PNRM = (R1), INDEX = (I1)'
              CALL xerrwd(msg, 50, 902, 0, 1, ivar, 0, 1, pnrm, 0.0d0)
              ENDIF
            IF (pnrm .LE. stptol) THEN
              iret = 1
              RETURN
              ENDIF
            go to 10
            ENDIF
         ENDIF
C
      slpi = (-two*f1nrm)*ratio
      rlmin = stptol/pnrm
      IF (lsoff .EQ. 0 .AND. kprin .GE. 2) THEN
        msg = '------ MIN. LAMBDA = (R1)'
        CALL xerrwd(msg, 25, 903, 0, 0, 0, 0, 1, rlmin, 0.0d0)
        ENDIF
C-----------------------------------------------------------------------
C Begin iteration to find RL value satisfying alpha-condition.
C If RL becomes less than RLMIN, then terminate with IRET = 1.
C-----------------------------------------------------------------------
 100  CONTINUE
      CALL dyypnw(neq,y,yprime,cj,rl,p,icopt,id,ynew,ypnew)
      CALL dfnrmd(neq, ynew, t, ypnew, r, cj, wt, res, ires,
     *              fnrmp, wm, iwm, rpar, ipar)
      iwm(lnre) = iwm(lnre) + 1
      IF (ires .NE. 0) THEN
        iret = 2
        RETURN
        ENDIF
      IF (lsoff .EQ. 1) go to 150
C
      f1nrmp = fnrmp*fnrmp/two
      IF (kprin .GE. 2) THEN
        msg = '------ LAMBDA = (R1)'
        CALL xerrwd(msg, 20, 904, 0, 0, 0, 0, 1, rl, 0.0d0)
        msg = '------ NORM(F1) = (R1),  NORM(F1NEW) = (R2)'
        CALL xerrwd(msg, 43, 905, 0, 0, 0, 0, 2, f1nrm, f1nrmp)
        ENDIF
      IF (f1nrmp .GT. f1nrm + alpha*slpi*rl) go to 200
C-----------------------------------------------------------------------
C Alpha-condition is satisfied, or linesearch is turned off.
C Copy YNEW,YPNEW to Y,YPRIME and return.
C-----------------------------------------------------------------------
 150  iret = 0
      CALL dcopy(neq, ynew, 1, y, 1)
      CALL dcopy(neq, ypnew, 1, yprime, 1)
      fnrm = fnrmp
      IF (kprin .GE. 1) THEN
        msg = '------ LEAVING ROUTINE DLINSD, FNRM = (R1)'
        CALL xerrwd(msg, 42, 906, 0, 0, 0, 0, 1, fnrm, 0.0d0)
        ENDIF
      RETURN
C-----------------------------------------------------------------------
C Alpha-condition not satisfied.  Perform backtrack to compute new RL
C value.  If no satisfactory YNEW,YPNEW can be found sufficiently 
C distinct from Y,YPRIME, then return IRET = 1.
C-----------------------------------------------------------------------
 200  CONTINUE
      IF (rl .LT. rlmin) THEN
        iret = 1
        RETURN
        ENDIF
C
      rl = rl/two
      go to 100
C
C----------------------- END OF SUBROUTINE DLINSD ----------------------
      END