GNU Octave  4.2.1 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
i e
Definition: data.cc:2724
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
double h
Definition: graphics.cc:11205
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to