GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
ddasic.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 ddasic (X, Y, YPRIME, NEQ, ICOPT, ID, RES, JAC, PSOL,
6  * H, WT, NIC, IDID, RPAR, IPAR, PHI, SAVR, DELTA, E, YIC, YPIC,
7  * PWK, WM, IWM, HMIN, UROUND, EPLI, SQRTN, RSQRTN, EPCONI,
8  * STPTOL, JFLG, ICNFLG, ICNSTR, NLSIC)
9 C
10 C***BEGIN PROLOGUE DDASIC
11 C***REFER TO DDASPK
12 C***DATE WRITTEN 940628 (YYMMDD)
13 C***REVISION DATE 941206 (YYMMDD)
14 C***REVISION DATE 950714 (YYMMDD)
15 C
16 C-----------------------------------------------------------------------
17 C***DESCRIPTION
18 C
19 C DDASIC is a driver routine to compute consistent initial values
20 C for Y and YPRIME. There are two different options:
21 C Denoting the differential variables in Y by Y_d, and
22 C the algebraic variables by Y_a, the problem solved is either:
23 C 1. Given Y_d, calculate Y_a and Y_d', or
24 C 2. Given Y', calculate Y.
25 C In either case, initial values for the given components
26 C are input, and initial guesses for the unknown components
27 C must also be provided as input.
28 C
29 C The external routine NLSIC solves the resulting nonlinear system.
30 C
31 C The parameters represent
32 C
33 C X -- Independent variable.
34 C Y -- Solution vector at X.
35 C YPRIME -- Derivative of solution vector.
36 C NEQ -- Number of equations to be integrated.
37 C ICOPT -- Flag indicating initial condition option chosen.
38 C ICOPT = 1 for option 1 above.
39 C ICOPT = 2 for option 2.
40 C ID -- Array of dimension NEQ, which must be initialized
41 C if option 1 is chosen.
42 C ID(i) = +1 if Y_i is a differential variable,
43 C ID(i) = -1 if Y_i is an algebraic variable.
44 C RES -- External user-supplied subroutine to evaluate the
45 C residual. See RES description in DDASPK prologue.
46 C JAC -- External user-supplied routine to update Jacobian
47 C or preconditioner information in the nonlinear solver
48 C (optional). See JAC description in DDASPK prologue.
49 C PSOL -- External user-supplied routine to solve
50 C a linear system using preconditioning.
51 C See PSOL in DDASPK prologue.
52 C H -- Scaling factor in iteration matrix. DDASIC may
53 C reduce H to achieve convergence.
54 C WT -- Vector of weights for error criterion.
55 C NIC -- Input number of initial condition calculation call
56 C (= 1 or 2).
57 C IDID -- Completion code. See IDID in DDASPK prologue.
58 C RPAR,IPAR -- Real and integer parameter arrays that
59 C are used for communication between the
60 C calling program and external user routines.
61 C They are not altered by DNSK
62 C PHI -- Work space for DDASIC of length at least 2*NEQ.
63 C SAVR -- Work vector for DDASIC of length NEQ.
64 C DELTA -- Work vector for DDASIC of length NEQ.
65 C E -- Work vector for DDASIC of length NEQ.
66 C YIC,YPIC -- Work vectors for DDASIC, each of length NEQ.
67 C PWK -- Work vector for DDASIC of length NEQ.
68 C WM,IWM -- Real and integer arrays storing
69 C information required by the linear solver.
70 C EPCONI -- Test constant for Newton iteration convergence.
71 C ICNFLG -- Flag showing whether constraints on Y are to apply.
72 C ICNSTR -- Integer array of length NEQ with constraint types.
73 C
74 C The other parameters are for use internally by DDASIC.
75 C
76 C-----------------------------------------------------------------------
77 C***ROUTINES CALLED
78 C DCOPY, NLSIC
79 C
80 C***END PROLOGUE DDASIC
81 C
82 C
83  IMPLICIT DOUBLE PRECISION(a-h,o-z)
84  dimension y(*),yprime(*),id(*),wt(*),phi(neq,*)
85  dimension savr(*),delta(*),e(*),yic(*),ypic(*),pwk(*)
86  dimension wm(*),iwm(*), rpar(*),ipar(*), icnstr(*)
87  EXTERNAL res, jac, psol, nlsic
88 C
89  parameter(lcfn=15)
90  parameter(lmxnh=34)
91 C
92 C The following parameters are data-loaded here:
93 C RHCUT = factor by which H is reduced on retry of Newton solve.
94 C RATEMX = maximum convergence rate for which Newton iteration
95 C is considered converging.
96 C
97  SAVE rhcut, ratemx
98  DATA rhcut/0.1d0/, ratemx/0.8d0/
99 C
100 C
101 C-----------------------------------------------------------------------
102 C BLOCK 1.
103 C Initializations.
104 C JSKIP is a flag set to 1 when NIC = 2 and NH = 1, to signal that
105 C the initial call to the JAC routine is to be skipped then.
106 C Save Y and YPRIME in PHI. Initialize IDID, NH, and CJ.
107 C-----------------------------------------------------------------------
108 C
109  mxnh = iwm(lmxnh)
110  idid = 1
111  nh = 1
112  jskip = 0
113  IF (nic .EQ. 2) jskip = 1
114  CALL dcopy (neq, y, 1, phi(1,1), 1)
115  CALL dcopy (neq, yprime, 1, phi(1,2), 1)
116 C
117  IF (icopt .EQ. 2) THEN
118  cj = 0.0d0
119  ELSE
120  cj = 1.0d0/h
121  ENDIF
122 C
123 C-----------------------------------------------------------------------
124 C BLOCK 2
125 C Call the nonlinear system solver to obtain
126 C consistent initial values for Y and YPRIME.
127 C-----------------------------------------------------------------------
128 C
129  200 CONTINUE
130  CALL nlsic(x,y,yprime,neq,icopt,id,res,jac,psol,h,wt,jskip,
131  * rpar,ipar,savr,delta,e,yic,ypic,pwk,wm,iwm,cj,uround,
132  * epli,sqrtn,rsqrtn,epconi,ratemx,stptol,jflg,icnflg,icnstr,
133  * iernls)
134 C
135  IF (iernls .EQ. 0) RETURN
136 C
137 C-----------------------------------------------------------------------
138 C BLOCK 3
139 C The nonlinear solver was unsuccessful. Increment NCFN.
140 C Return with IDID = -12 if either
141 C IERNLS = -1: error is considered unrecoverable,
142 C ICOPT = 2: we are doing initialization problem type 2, or
143 C NH = MXNH: the maximum number of H values has been tried.
144 C Otherwise (problem 1 with IERNLS .GE. 1), reduce H and try again.
145 C If IERNLS > 1, restore Y and YPRIME to their original values.
146 C-----------------------------------------------------------------------
147 C
148  iwm(lcfn) = iwm(lcfn) + 1
149  jskip = 0
150 C
151  IF (iernls .EQ. -1) GO TO 350
152  IF (icopt .EQ. 2) GO TO 350
153  IF (nh .EQ. mxnh) GO TO 350
154 C
155  nh = nh + 1
156  h = h*rhcut
157  cj = 1.0d0/h
158 C
159  IF (iernls .EQ. 1) GO TO 200
160 C
161  CALL dcopy (neq, phi(1,1), 1, y, 1)
162  CALL dcopy (neq, phi(1,2), 1, yprime, 1)
163  GO TO 200
164 C
165  350 idid = -12
166  RETURN
167 C
168 C------END OF SUBROUTINE DDASIC-----------------------------------------
169  END
ddasic
subroutine ddasic(X, Y, YPRIME, NEQ, ICOPT, ID, RES, JAC, PSOL, H, WT, NIC, IDID, RPAR, IPAR, PHI, SAVR, DELTA, E, YIC, YPIC, PWK, WM, IWM, HMIN, UROUND, EPLI, SQRTN, RSQRTN, EPCONI, STPTOL, JFLG, ICNFLG, ICNSTR, NLSIC)
Definition: ddasic.f:9