|
3911
|
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 |
