GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dcnstr.f
Go to the documentation of this file.
1 C Work performed under the auspices of the U.S. Department of Energy
2 C by Lawrence Livermore National Laboratory under contract number
3 C W-7405-Eng-48.
4 C
5  SUBROUTINE dcnstr (NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR)
6 C
7 C***BEGIN PROLOGUE DCNSTR
8 C***DATE WRITTEN 950808 (YYMMDD)
9 C***REVISION DATE 950814 (YYMMDD)
10 C
11 C
12 C-----------------------------------------------------------------------
13 C***DESCRIPTION
14 C
15 C This subroutine checks for constraint violations in the proposed
16 C new approximate solution YNEW.
17 C If a constraint violation occurs, then a new step length, TAU,
18 C is calculated, and this value is to be given to the linesearch routine
19 C to calculate a new approximate solution YNEW.
20 C
21 C On entry:
22 C
23 C NEQ -- size of the nonlinear system, and the length of arrays
24 C Y, YNEW and ICNSTR.
25 C
26 C Y -- real array containing the current approximate y.
27 C
28 C YNEW -- real array containing the new approximate y.
29 C
30 C ICNSTR -- INTEGER array of length NEQ containing flags indicating
31 C which entries in YNEW are to be constrained.
32 C if ICNSTR(I) = 2, then YNEW(I) must be .GT. 0,
33 C if ICNSTR(I) = 1, then YNEW(I) must be .GE. 0,
34 C if ICNSTR(I) = -1, then YNEW(I) must be .LE. 0, while
35 C if ICNSTR(I) = -2, then YNEW(I) must be .LT. 0, while
36 C if ICNSTR(I) = 0, then YNEW(I) is not constrained.
37 C
38 C RLX -- real scalar restricting update, if ICNSTR(I) = 2 or -2,
39 C to ABS( (YNEW-Y)/Y ) < FAC2*RLX in component I.
40 C
41 C TAU -- the current size of the step length for the linesearch.
42 C
43 C On return
44 C
45 C TAU -- the adjusted size of the step length if a constraint
46 C violation occurred (otherwise, it is unchanged). it is
47 C the step length to give to the linesearch routine.
48 C
49 C IRET -- output flag.
50 C IRET=0 means that YNEW satisfied all constraints.
51 C IRET=1 means that YNEW failed to satisfy all the
52 C constraints, and a new linesearch step
53 C must be computed.
54 C
55 C IVAR -- index of variable causing constraint to be violated.
56 C
57 C-----------------------------------------------------------------------
58  IMPLICIT DOUBLE PRECISION(a-h,o-z)
59  dimension y(neq), ynew(neq), icnstr(neq)
60  SAVE fac, fac2, zero
61  DATA fac /0.6d0/, fac2 /0.9d0/, zero/0.0d0/
62 C-----------------------------------------------------------------------
63 C Check constraints for proposed new step YNEW. If a constraint has
64 C been violated, then calculate a new step length, TAU, to be
65 C used in the linesearch routine.
66 C-----------------------------------------------------------------------
67  iret = 0
68  rdymx = zero
69  ivar = 0
70  DO 100 i = 1,neq
71 C
72  IF (icnstr(i) .EQ. 2) THEN
73  rdy = abs( (ynew(i)-y(i))/y(i) )
74  IF (rdy .GT. rdymx) THEN
75  rdymx = rdy
76  ivar = i
77  ENDIF
78  IF (ynew(i) .LE. zero) THEN
79  tau = fac*tau
80  ivar = i
81  iret = 1
82  RETURN
83  ENDIF
84 C
85  ELSEIF (icnstr(i) .EQ. 1) THEN
86  IF (ynew(i) .LT. zero) THEN
87  tau = fac*tau
88  ivar = i
89  iret = 1
90  RETURN
91  ENDIF
92 C
93  ELSEIF (icnstr(i) .EQ. -1) THEN
94  IF (ynew(i) .GT. zero) THEN
95  tau = fac*tau
96  ivar = i
97  iret = 1
98  RETURN
99  ENDIF
100 C
101  ELSEIF (icnstr(i) .EQ. -2) THEN
102  rdy = abs( (ynew(i)-y(i))/y(i) )
103  IF (rdy .GT. rdymx) THEN
104  rdymx = rdy
105  ivar = i
106  ENDIF
107  IF (ynew(i) .GE. zero) THEN
108  tau = fac*tau
109  ivar = i
110  iret = 1
111  RETURN
112  ENDIF
113 C
114  ENDIF
115  100 CONTINUE
116 
117  IF(rdymx .GE. rlx) THEN
118  tau = fac2*tau*rlx/rdymx
119  iret = 1
120  ENDIF
121 C
122  RETURN
123 C----------------------- END OF SUBROUTINE DCNSTR ----------------------
124  END
subroutine dcnstr(NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR)
Definition: dcnstr.f:6