|
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 DLINSD (NEQ, Y, T, YPRIME, CJ, P, PNRM, WT, LSOFF, |
|
|
6 * STPTOL, IRET, RES, IRES, WM, IWM, |
|
|
7 * FNRM, ICOPT, ID, R, YNEW, YPNEW, ICNFLG, |
|
|
8 * ICNSTR, RLX, RPAR, IPAR) |
|
|
9 C |
|
|
10 C***BEGIN PROLOGUE DLINSD |
|
|
11 C***REFER TO DNSID |
|
|
12 C***DATE WRITTEN 941025 (YYMMDD) |
|
|
13 C***REVISION DATE 941215 (YYMMDD) |
|
|
14 C***REVISION DATE 960129 Moved line RL = ONE to top block. |
|
|
15 C |
|
|
16 C |
|
|
17 C----------------------------------------------------------------------- |
|
|
18 C***DESCRIPTION |
|
|
19 C |
|
|
20 C DLINSD uses a linesearch algorithm to calculate a new (Y,YPRIME) |
|
|
21 C pair (YNEW,YPNEW) such that |
|
|
22 C |
|
|
23 C f(YNEW,YPNEW) .le. (1 - 2*ALPHA*RL)*f(Y,YPRIME) , |
|
|
24 C |
|
|
25 C where 0 < RL <= 1. Here, f(y,y') is defined as |
|
|
26 C |
|
|
27 C f(y,y') = (1/2)*norm( (J-inverse)*G(t,y,y') )**2 , |
|
|
28 C |
|
|
29 C where norm() is the weighted RMS vector norm, G is the DAE |
|
|
30 C system residual function, and J is the system iteration matrix |
|
|
31 C (Jacobian). |
|
|
32 C |
|
|
33 C In addition to the parameters defined elsewhere, we have |
|
|
34 C |
|
|
35 C P -- Approximate Newton step used in backtracking. |
|
|
36 C PNRM -- Weighted RMS norm of P. |
|
|
37 C LSOFF -- Flag showing whether the linesearch algorithm is |
|
|
38 C to be invoked. 0 means do the linesearch, and |
|
|
39 C 1 means turn off linesearch. |
|
|
40 C STPTOL -- Tolerance used in calculating the minimum lambda |
|
|
41 C value allowed. |
|
|
42 C ICNFLG -- Integer scalar. If nonzero, then constraint violations |
|
|
43 C in the proposed new approximate solution will be |
|
|
44 C checked for, and the maximum step length will be |
|
|
45 C adjusted accordingly. |
|
|
46 C ICNSTR -- Integer array of length NEQ containing flags for |
|
|
47 C checking constraints. |
|
|
48 C RLX -- Real scalar restricting update size in DCNSTR. |
|
|
49 C YNEW -- Array of length NEQ used to hold the new Y in |
|
|
50 C performing the linesearch. |
|
|
51 C YPNEW -- Array of length NEQ used to hold the new YPRIME in |
|
|
52 C performing the linesearch. |
|
|
53 C Y -- Array of length NEQ containing the new Y (i.e.,=YNEW). |
|
|
54 C YPRIME -- Array of length NEQ containing the new YPRIME |
|
|
55 C (i.e.,=YPNEW). |
|
|
56 C FNRM -- Real scalar containing SQRT(2*f(Y,YPRIME)) for the |
|
|
57 C current (Y,YPRIME) on input and output. |
|
|
58 C R -- Work array of length NEQ, containing the scaled |
|
|
59 C residual (J-inverse)*G(t,y,y') on return. |
|
|
60 C IRET -- Return flag. |
|
|
61 C IRET=0 means that a satisfactory (Y,YPRIME) was found. |
|
|
62 C IRET=1 means that the routine failed to find a new |
|
|
63 C (Y,YPRIME) that was sufficiently distinct from |
|
|
64 C the current (Y,YPRIME) pair. |
|
|
65 C IRET=2 means IRES .ne. 0 from RES. |
|
|
66 C----------------------------------------------------------------------- |
|
|
67 C |
|
|
68 C***ROUTINES CALLED |
|
|
69 C DFNRMD, DYYPNW, DCOPY |
|
|
70 C |
|
|
71 C***END PROLOGUE DLINSD |
|
|
72 C |
|
|
73 IMPLICIT DOUBLE PRECISION(A-H,O-Z) |
|
|
74 EXTERNAL RES |
|
|
75 DIMENSION Y(*), YPRIME(*), WT(*), R(*), ID(*) |
|
|
76 DIMENSION WM(*), IWM(*) |
|
|
77 DIMENSION YNEW(*), YPNEW(*), P(*), ICNSTR(*) |
|
|
78 DIMENSION RPAR(*), IPAR(*) |
|
|
79 CHARACTER MSG*80 |
|
|
80 C |
|
|
81 PARAMETER (LNRE=12, LKPRIN=31) |
|
|
82 C |
|
|
83 SAVE ALPHA, ONE, TWO |
|
|
84 DATA ALPHA/1.0D-4/, ONE/1.0D0/, TWO/2.0D0/ |
|
|
85 C |
|
|
86 KPRIN=IWM(LKPRIN) |
|
|
87 C |
|
|
88 F1NRM = (FNRM*FNRM)/TWO |
|
|
89 RATIO = ONE |
|
|
90 IF (KPRIN .GE. 2) THEN |
|
|
91 MSG = '------ IN ROUTINE DLINSD-- PNRM = (R1) )' |
|
|
92 CALL XERRWD(MSG, 40, 901, 0, 0, 0, 0, 1, PNRM, 0.0D0) |
|
|
93 ENDIF |
|
|
94 TAU = PNRM |
|
|
95 IVIO = 0 |
|
|
96 RL = ONE |
|
|
97 C----------------------------------------------------------------------- |
|
|
98 C Check for violations of the constraints, if any are imposed. |
|
|
99 C If any violations are found, the step vector P is rescaled, and the |
|
|
100 C constraint check is repeated, until no violations are found. |
|
|
101 C----------------------------------------------------------------------- |
|
|
102 IF (ICNFLG .NE. 0) THEN |
|
|
103 10 CONTINUE |
|
|
104 CALL DYYPNW (NEQ,Y,YPRIME,CJ,RL,P,ICOPT,ID,YNEW,YPNEW) |
|
|
105 CALL DCNSTR (NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR) |
|
|
106 IF (IRET .EQ. 1) THEN |
|
|
107 IVIO = 1 |
|
|
108 RATIO1 = TAU/PNRM |
|
|
109 RATIO = RATIO*RATIO1 |
|
|
110 DO 20 I = 1,NEQ |
|
|
111 20 P(I) = P(I)*RATIO1 |
|
|
112 PNRM = TAU |
|
|
113 IF (KPRIN .GE. 2) THEN |
|
|
114 MSG = '------ CONSTRAINT VIOL., PNRM = (R1), INDEX = (I1)' |
|
|
115 CALL XERRWD(MSG, 50, 902, 0, 1, IVAR, 0, 1, PNRM, 0.0D0) |
|
|
116 ENDIF |
|
|
117 IF (PNRM .LE. STPTOL) THEN |
|
|
118 IRET = 1 |
|
|
119 RETURN |
|
|
120 ENDIF |
|
|
121 GO TO 10 |
|
|
122 ENDIF |
|
|
123 ENDIF |
|
|
124 C |
|
|
125 SLPI = (-TWO*F1NRM)*RATIO |
|
|
126 RLMIN = STPTOL/PNRM |
|
|
127 IF (LSOFF .EQ. 0 .AND. KPRIN .GE. 2) THEN |
|
|
128 MSG = '------ MIN. LAMBDA = (R1)' |
|
|
129 CALL XERRWD(MSG, 25, 903, 0, 0, 0, 0, 1, RLMIN, 0.0D0) |
|
|
130 ENDIF |
|
|
131 C----------------------------------------------------------------------- |
|
|
132 C Begin iteration to find RL value satisfying alpha-condition. |
|
|
133 C If RL becomes less than RLMIN, then terminate with IRET = 1. |
|
|
134 C----------------------------------------------------------------------- |
|
|
135 100 CONTINUE |
|
|
136 CALL DYYPNW (NEQ,Y,YPRIME,CJ,RL,P,ICOPT,ID,YNEW,YPNEW) |
|
|
137 CALL DFNRMD (NEQ, YNEW, T, YPNEW, R, CJ, WT, RES, IRES, |
|
|
138 * FNRMP, WM, IWM, RPAR, IPAR) |
|
|
139 IWM(LNRE) = IWM(LNRE) + 1 |
|
|
140 IF (IRES .NE. 0) THEN |
|
|
141 IRET = 2 |
|
|
142 RETURN |
|
|
143 ENDIF |
|
|
144 IF (LSOFF .EQ. 1) GO TO 150 |
|
|
145 C |
|
|
146 F1NRMP = FNRMP*FNRMP/TWO |
|
|
147 IF (KPRIN .GE. 2) THEN |
|
|
148 MSG = '------ LAMBDA = (R1)' |
|
|
149 CALL XERRWD(MSG, 20, 904, 0, 0, 0, 0, 1, RL, 0.0D0) |
|
|
150 MSG = '------ NORM(F1) = (R1), NORM(F1NEW) = (R2)' |
|
|
151 CALL XERRWD(MSG, 43, 905, 0, 0, 0, 0, 2, F1NRM, F1NRMP) |
|
|
152 ENDIF |
|
|
153 IF (F1NRMP .GT. F1NRM + ALPHA*SLPI*RL) GO TO 200 |
|
|
154 C----------------------------------------------------------------------- |
|
|
155 C Alpha-condition is satisfied, or linesearch is turned off. |
|
|
156 C Copy YNEW,YPNEW to Y,YPRIME and return. |
|
|
157 C----------------------------------------------------------------------- |
|
|
158 150 IRET = 0 |
|
|
159 CALL DCOPY (NEQ, YNEW, 1, Y, 1) |
|
|
160 CALL DCOPY (NEQ, YPNEW, 1, YPRIME, 1) |
|
|
161 FNRM = FNRMP |
|
|
162 IF (KPRIN .GE. 1) THEN |
|
|
163 MSG = '------ LEAVING ROUTINE DLINSD, FNRM = (R1)' |
|
|
164 CALL XERRWD(MSG, 42, 906, 0, 0, 0, 0, 1, FNRM, 0.0D0) |
|
|
165 ENDIF |
|
|
166 RETURN |
|
|
167 C----------------------------------------------------------------------- |
|
|
168 C Alpha-condition not satisfied. Perform backtrack to compute new RL |
|
|
169 C value. If no satisfactory YNEW,YPNEW can be found sufficiently |
|
|
170 C distinct from Y,YPRIME, then return IRET = 1. |
|
|
171 C----------------------------------------------------------------------- |
|
|
172 200 CONTINUE |
|
|
173 IF (RL .LT. RLMIN) THEN |
|
|
174 IRET = 1 |
|
|
175 RETURN |
|
|
176 ENDIF |
|
|
177 C |
|
|
178 RL = RL/TWO |
|
|
179 GO TO 100 |
|
|
180 C |
|
|
181 C----------------------- END OF SUBROUTINE DLINSD ---------------------- |
|
|
182 END |
