GNU Octave  3.8.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
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