dcnstr.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 dcnstr (NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR)
C
C***BEGIN PROLOGUE  DCNSTR
C***DATE WRITTEN   950808   (YYMMDD)
C***REVISION DATE  950814   (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This subroutine checks for constraint violations in the proposed 
C new approximate solution YNEW.
C If a constraint violation occurs, then a new step length, TAU,
C is calculated, and this value is to be given to the linesearch routine
C to calculate a new approximate solution YNEW.
C
C On entry:
C
C   NEQ    -- size of the nonlinear system, and the length of arrays
C             Y, YNEW and ICNSTR.
C
C   Y      -- real array containing the current approximate y.
C
C   YNEW   -- real array containing the new approximate y.
C
C   ICNSTR -- INTEGER array of length NEQ containing flags indicating
C             which entries in YNEW are to be constrained.
C             if ICNSTR(I) =  2, then YNEW(I) must be .GT. 0,
C             if ICNSTR(I) =  1, then YNEW(I) must be .GE. 0,
C             if ICNSTR(I) = -1, then YNEW(I) must be .LE. 0, while
C             if ICNSTR(I) = -2, then YNEW(I) must be .LT. 0, while
C             if ICNSTR(I) =  0, then YNEW(I) is not constrained.
C
C   RLX    -- real scalar restricting update, if ICNSTR(I) = 2 or -2,
C             to ABS( (YNEW-Y)/Y ) < FAC2*RLX in component I.
C
C   TAU    -- the current size of the step length for the linesearch.
C
C On return
C
C   TAU    -- the adjusted size of the step length if a constraint
C             violation occurred (otherwise, it is unchanged).  it is
C             the step length to give to the linesearch routine.
C
C   IRET   -- output flag.
C             IRET=0 means that YNEW satisfied all constraints.
C             IRET=1 means that YNEW failed to satisfy all the
C                    constraints, and a new linesearch step
C                    must be computed.
C
C   IVAR   -- index of variable causing constraint to be violated.
C
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION(a-h,o-z)
      dimension y(neq), ynew(neq), icnstr(neq)
      SAVE fac, fac2, zero
      DATA fac /0.6d0/, fac2 /0.9d0/, zero/0.0d0/
C-----------------------------------------------------------------------
C Check constraints for proposed new step YNEW.  If a constraint has
C been violated, then calculate a new step length, TAU, to be
C used in the linesearch routine.
C-----------------------------------------------------------------------
      iret = 0
      rdymx = zero
      ivar = 0
      DO 100 i = 1,neq
C
         IF (icnstr(i) .EQ. 2) THEN
            rdy = abs( (ynew(i)-y(i))/y(i) )
            IF (rdy .GT. rdymx) THEN
               rdymx = rdy
               ivar = i
            ENDIF
            IF (ynew(i) .LE. zero) THEN
               tau = fac*tau
               ivar = i
               iret = 1
               RETURN
            ENDIF
C
         ELSEIF (icnstr(i) .EQ. 1) THEN
            IF (ynew(i) .LT. zero) THEN
               tau = fac*tau
               ivar = i
               iret = 1
               RETURN
            ENDIF
C
         ELSEIF (icnstr(i) .EQ. -1) THEN
            IF (ynew(i) .GT. zero) THEN
               tau = fac*tau
               ivar = i
               iret = 1
               RETURN
            ENDIF
C
         ELSEIF (icnstr(i) .EQ. -2) THEN
            rdy = abs( (ynew(i)-y(i))/y(i) )
            IF (rdy .GT. rdymx) THEN
               rdymx = rdy
               ivar = i
            ENDIF
            IF (ynew(i) .GE. zero) THEN
               tau = fac*tau
               ivar = i
               iret = 1
               RETURN
            ENDIF
C
         ENDIF
 100  CONTINUE

      IF(rdymx .GE. rlx) THEN
         tau = fac2*tau*rlx/rdymx
         iret = 1
      ENDIF
C
      RETURN
C----------------------- END OF SUBROUTINE DCNSTR ----------------------
      END