GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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