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