GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dnsik.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 dnsik(X,Y,YPRIME,NEQ,ICOPT,ID,RES,PSOL,WT,RPAR,IPAR,
6  * SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,SQRTN,RSQRTN,EPLIN,EPCON,
7  * RATEMX,MAXIT,STPTOL,ICNFLG,ICNSTR,IERNEW)
8 C
9 C***BEGIN PROLOGUE DNSIK
10 C***REFER TO DDASPK
11 C***DATE WRITTEN 940701 (YYMMDD)
12 C***REVISION DATE 950714 (YYMMDD)
13 C
14 C
15 C-----------------------------------------------------------------------
16 C***DESCRIPTION
17 C
18 C DNSIK solves a nonlinear system of algebraic equations of the
19 C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in
20 C the initial conditions.
21 C
22 C The method used is a Newton scheme combined with a linesearch
23 C algorithm, using Krylov iterative linear system methods.
24 C
25 C The parameters represent
26 C
27 C X -- Independent variable.
28 C Y -- Solution vector.
29 C YPRIME -- Derivative of solution vector.
30 C NEQ -- Number of unknowns.
31 C ICOPT -- Initial condition option chosen (1 or 2).
32 C ID -- Array of dimension NEQ, which must be initialized
33 C if ICOPT = 1. See DDASIC.
34 C RES -- External user-supplied subroutine
35 C to evaluate the residual. See RES description
36 C in DDASPK prologue.
37 C PSOL -- External user-supplied routine to solve
38 C a linear system using preconditioning.
39 C See explanation inside DDASPK.
40 C WT -- Vector of weights for error criterion.
41 C RPAR,IPAR -- Real and integer arrays used for communication
42 C between the calling program and external user
43 C routines. They are not altered within DASPK.
44 C SAVR -- Work vector for DNSIK of length NEQ.
45 C DELTA -- Residual vector on entry, and work vector of
46 C length NEQ for DNSIK.
47 C R -- Work vector for DNSIK of length NEQ.
48 C YIC,YPIC -- Work vectors for DNSIK, each of length NEQ.
49 C PWK -- Work vector for DNSIK of length NEQ.
50 C WM,IWM -- Real and integer arrays storing
51 C matrix information such as the matrix
52 C of partial derivatives, permutation
53 C vector, and various other information.
54 C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).
55 C SQRTN -- Square root of NEQ.
56 C RSQRTN -- reciprical of square root of NEQ.
57 C EPLIN -- Tolerance for linear system solver.
58 C EPCON -- Tolerance to test for convergence of the Newton
59 C iteration.
60 C RATEMX -- Maximum convergence rate for which Newton iteration
61 C is considered converging.
62 C MAXIT -- Maximum allowed number of Newton iterations.
63 C STPTOL -- Tolerance used in calculating the minimum lambda
64 C value allowed.
65 C ICNFLG -- Integer scalar. If nonzero, then constraint
66 C violations in the proposed new approximate solution
67 C will be checked for, and the maximum step length
68 C will be adjusted accordingly.
69 C ICNSTR -- Integer array of length NEQ containing flags for
70 C checking constraints.
71 C IERNEW -- Error flag for Newton iteration.
72 C 0 ==> Newton iteration converged.
73 C 1 ==> failed to converge, but RATE .lt. 1.
74 C 2 ==> failed to converge, RATE .gt. RATEMX.
75 C 3 ==> other recoverable error.
76 C -1 ==> unrecoverable error inside Newton iteration.
77 C-----------------------------------------------------------------------
78 C
79 C***ROUTINES CALLED
80 C DFNRMK, DSLVK, DDWNRM, DLINSK, DCOPY
81 C
82 C***END PROLOGUE DNSIK
83 C
84 C
85  IMPLICIT DOUBLE PRECISION(a-h,o-z)
86  dimension y(*),yprime(*),wt(*),id(*),delta(*),r(*),savr(*)
87  dimension yic(*),ypic(*),pwk(*),wm(*),iwm(*), rpar(*),ipar(*)
88  dimension icnstr(*)
89  EXTERNAL res, psol
90 C
91  parameter(lnni=19, lnps=21, llocwp=29, llciwp=30)
92  parameter(llsoff=35, lstol=14)
93 C
94 C
95 C Initializations. M is the Newton iteration counter.
96 C
97  lsoff = iwm(llsoff)
98  m = 0
99  rate = 1.0d0
100  lwp = iwm(llocwp)
101  liwp = iwm(llciwp)
102  rlx = 0.4d0
103 C
104 C Save residual in SAVR.
105 C
106  CALL dcopy (neq, delta, 1, savr, 1)
107 C
108 C Compute norm of (P-inverse)*(residual).
109 C
110  CALL dfnrmk (neq, y, x, yprime, savr, r, cj, wt, sqrtn, rsqrtn,
111  * res, ires, psol, 1, ier, fnrm, eplin, wm(lwp), iwm(liwp),
112  * pwk, rpar, ipar)
113  iwm(lnps) = iwm(lnps) + 1
114  IF (ier .NE. 0) THEN
115  iernew = 3
116  RETURN
117  ENDIF
118 C
119 C Return now if residual norm is .le. EPCON.
120 C
121  IF (fnrm .LE. epcon) RETURN
122 C
123 C Newton iteration loop.
124 C
125 300 CONTINUE
126  iwm(lnni) = iwm(lnni) + 1
127 C
128 C Compute a new step vector DELTA.
129 C
130  CALL dslvk (neq, y, x, yprime, savr, delta, wt, wm, iwm,
131  * res, ires, psol, iersl, cj, eplin, sqrtn, rsqrtn, rhok,
132  * rpar, ipar)
133  IF (ires .NE. 0 .OR. iersl .NE. 0) GO TO 390
134 C
135 C Get norm of DELTA. Return now if DELTA is zero.
136 C
137  delnrm = ddwnrm(neq,delta,wt,rpar,ipar)
138  IF (delnrm .EQ. 0.0d0) RETURN
139 C
140 C Call linesearch routine for global strategy and set RATE.
141 C
142  oldfnm = fnrm
143 C
144  CALL dlinsk (neq, y, x, yprime, savr, cj, delta, delnrm, wt,
145  * sqrtn, rsqrtn, lsoff, stptol, iret, res, ires, psol, wm, iwm,
146  * rhok, fnrm, icopt, id, wm(lwp), iwm(liwp), r, eplin, yic, ypic,
147  * pwk, icnflg, icnstr, rlx, rpar, ipar)
148 C
149  rate = fnrm/oldfnm
150 C
151 C Check for error condition from linesearch.
152  IF (iret .NE. 0) GO TO 390
153 C
154 C Test for convergence of the iteration, and return or loop.
155 C
156  IF (fnrm .LE. epcon) RETURN
157 C
158 C The iteration has not yet converged. Update M.
159 C Test whether the maximum number of iterations have been tried.
160 C
161  m=m+1
162  IF(m .GE. maxit) GO TO 380
163 C
164 C Copy the residual SAVR to DELTA and loop for another iteration.
165 C
166  CALL dcopy (neq, savr, 1, delta, 1)
167  GO TO 300
168 C
169 C The maximum number of iterations was done. Set IERNEW and return.
170 C
171 380 IF (rate .LE. ratemx) THEN
172  iernew = 1
173  ELSE
174  iernew = 2
175  ENDIF
176  RETURN
177 C
178 390 IF (ires .LE. -2 .OR. iersl .LT. 0) THEN
179  iernew = -1
180  ELSE
181  iernew = 3
182  IF (ires .EQ. 0 .AND. iersl .EQ. 1 .AND. m .GE. 2
183  1 .AND. rate .LT. 1.0d0) iernew = 1
184  ENDIF
185  RETURN
186 C
187 C
188 C----------------------- END OF SUBROUTINE DNSIK------------------------
189  END
ddwnrm
double precision function ddwnrm(NEQ, V, RWT, RPAR, IPAR)
Definition: ddwnrm.f:6
dfnrmk
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
dlinsk
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
dnsik
subroutine dnsik(X, Y, YPRIME, NEQ, ICOPT, ID, RES, PSOL, WT, RPAR, IPAR, SAVR, DELTA, R, YIC, YPIC, PWK, WM, IWM, CJ, SQRTN, RSQRTN, EPLIN, EPCON, RATEMX, MAXIT, STPTOL, ICNFLG, ICNSTR, IERNEW)
Definition: dnsik.f:8
dslvk
subroutine dslvk(NEQ, Y, TN, YPRIME, SAVR, X, EWT, WM, IWM, RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK, RPAR, IPAR)
Definition: dslvk.f:8