dlinsk.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 dlinsk (NEQ, Y, T, YPRIME, SAVR, CJ, P, PNRM, WT,
     *   SQRTN, RSQRTN, LSOFF, STPTOL, IRET, RES, IRES, PSOL, WM, IWM,
     *   RHOK, FNRM, ICOPT, ID, WP, IWP, R, EPLIN, YNEW, YPNEW, PWK,
     *   ICNFLG, ICNSTR, RLX, RPAR, IPAR)
C
C***BEGIN PROLOGUE  DLINSK
C***REFER TO  DNSIK
C***DATE WRITTEN   940830   (YYMMDD)
C***REVISION DATE  951006   (Arguments SQRTN, RSQRTN added.)
C***REVISION DATE  960129   Moved line RL = ONE to top block.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C     DLINSK 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                          ALPHA*RL*RHOK*RHOK ,
C
C     where 0 < RL <= 1, and RHOK is the scaled preconditioned norm of
C     the final residual vector in the Krylov iteration.
C     Here, f(y,y') is defined as
C
C      f(y,y') = (1/2)*norm( (P-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 P is the preconditioner used
C     in the Krylov iteration.
C
C     In addition to the parameters defined elsewhere, we have
C
C     SAVR    -- Work array of length NEQ, containing the residual
C                vector G(t,y,y') on return.
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,
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     RHOK    -- Weighted norm of preconditioned Krylov residual.
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     PWK     -- Work vector of length NEQ for use in PSOL.
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 space length NEQ for residual vector.
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 a failure in RES or PSOL.
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C   DFNRMK, DYYPNW, DCOPY
C
C***END PROLOGUE  DLINSK
C
      IMPLICIT DOUBLE PRECISION(a-h,o-z)
      EXTERNAL  res, psol
      dimension y(*), yprime(*), p(*), wt(*), savr(*), r(*), id(*)
      dimension wm(*), iwm(*), ynew(*), ypnew(*), pwk(*), icnstr(*)
      dimension wp(*), iwp(*), rpar(*), ipar(*)
      CHARACTER MSG*80
C
      parameter(lnre=12, lnps=21, lkprin=31)
C
      SAVE alpha, one, two
      DATA alpha/1.0d-4/, one/1.0d0/, two/2.0d0/
C
      kprin=iwm(lkprin)
      f1nrm = (fnrm*fnrm)/two
      ratio = one
C
      IF (kprin .GE. 2) THEN
        msg = '------ IN ROUTINE DLINSK-- PNRM = (R1) )'
        CALL xerrwd(msg, 40, 921, 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, 922, 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 + rhok*rhok)*ratio
      rlmin = stptol/pnrm
      IF (lsoff .EQ. 0 .AND. kprin .GE. 2) THEN
        msg = '------ MIN. LAMBDA = (R1)'
        CALL xerrwd(msg, 25, 923, 0, 0, 0, 0, 1, rlmin, 0.0d0)
        ENDIF
C-----------------------------------------------------------------------
C Begin iteration to find RL value satisfying alpha-condition.
C Update YNEW and YPNEW, then compute norm of new scaled residual and
C perform alpha condition test.
C-----------------------------------------------------------------------
 100  CONTINUE
      CALL dyypnw (neq,y,yprime,cj,rl,p,icopt,id,ynew,ypnew)
      CALL dfnrmk (neq, ynew, t, ypnew, savr, r, cj, wt, sqrtn, rsqrtn,
     *  res, ires, psol, 0, ier, fnrmp, eplin, wp, iwp, pwk, rpar, ipar)
      iwm(lnre) = iwm(lnre) + 1
      IF (ires .GE. 0) iwm(lnps) = iwm(lnps) + 1
      IF (ires .NE. 0 .OR. ier .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, 924, 0, 0, 0, 0, 1, rl, 0.0d0)
        msg = '------ NORM(F1) = (R1),  NORM(F1NEW) = (R2)'
        CALL xerrwd(msg, 43, 925, 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 DLINSK, FNRM = (R1)'
        CALL xerrwd(msg, 42, 926, 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 RL is less than RLMIN, i.e. no satisfactory YNEW,YPNEW can
C be found sufficiently 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 DLINSK ----------------------
      END