GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dlinsk.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 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
subroutine dcnstr(NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR)
Definition: dcnstr.f:6
subroutine dfnrmk(NEQ, Y, T, YPRIME, SAVR, R, CJ, WT, SQRTN, RSQRTN, RES, IRES, PSOL, IRIN, IER, FNORM, EPLIN, WP, IWP, PWK, RPAR, IPAR)
Definition: dfnrmk.f:8
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)
Definition: dlinsk.f:9
subroutine dyypnw(NEQ, Y, YPRIME, CJ, RL, P, ICOPT, ID, YNEW, YPNEW)
Definition: dyypnw.f:7
subroutine xerrwd(MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
Definition: xerrwd.f:4