GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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