GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dlinsd.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 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
subroutine dcnstr(NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR)
Definition: dcnstr.f:6
subroutine dfnrmd(NEQ, Y, T, YPRIME, R, CJ, WT, RES, IRES, FNORM, WM, IWM, RPAR, IPAR)
Definition: dfnrmd.f:7
subroutine dlinsd(NEQ, Y, T, YPRIME, CJ, P, PNRM, WT, LSOFF, STPTOL, IRET, RES, IRES, WM, IWM, FNRM, ICOPT, ID, R, YNEW, YPNEW, ICNFLG, ICNSTR, RLX, RPAR, IPAR)
Definition: dlinsd.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