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