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