CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
ddebdf.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine ddebdf (DF, NEQ, T, Y, TOUT, INFO, RTOL, ATOL, IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, DJAC)
 
subroutine dlsod (DF, NEQ, T, Y, TOUT, RTOL, ATOL, IDID, YPOUT, YH, YH1, EWT, SAVF, ACOR, WM, IWM, DJAC, INTOUT, TSTOP, TOLFAC, DELSGN, RPAR, IPAR)
 
subroutine dstod (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, DF, DJAC, RPAR, IPAR)
 
subroutine dcfod (METH, ELCO, TESCO)
 
double precision function dvnrms (N, V, W)
 
subroutine dintyd (T, K, YH, NYH, DKY, IFLAG)
 
subroutine dpjac (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, DF, DJAC, RPAR, IPAR)
 
subroutine dslvs (WM, IWM, X, TEM)
 
subroutine dgbfa (ABD, LDA, N, ML, MU, IPVT, INFO)
 
subroutine dgbsl (ABD, LDA, N, ML, MU, IPVT, B, JOB)
 
subroutine dgefa (A, LDA, N, IPVT, INFO)
 
subroutine dgesl (A, LDA, N, IPVT, B, JOB)
 

Function/Subroutine Documentation

◆ dcfod()

subroutine dcfod ( integer  METH,
double precision, dimension(13,12)  ELCO,
double precision, dimension(3,12)  TESCO 
)
2107 C***BEGIN PROLOGUE DCFOD
2108 C***SUBSIDIARY
2109 C***PURPOSE Subsidiary to DDEBDF
2110 C***LIBRARY SLATEC
2111 C***TYPE DOUBLE PRECISION (CFOD-S, DCFOD-D)
2112 C***AUTHOR (UNKNOWN)
2113 C***DESCRIPTION
2114 C
2115 C DCFOD defines coefficients needed in the integrator package DDEBDF
2116 C
2117 C***SEE ALSO DDEBDF
2118 C***ROUTINES CALLED (NONE)
2119 C***REVISION HISTORY (YYMMDD)
2120 C 820301 DATE WRITTEN
2121 C 890911 Removed unnecessary intrinsics. (WRB)
2122 C 891214 Prologue converted to Version 4.0 format. (BAB)
2123 C 900328 Added TYPE section. (WRB)
2124 C***END PROLOGUE DCFOD
2125 C
2126 C
2127  INTEGER i, ib, meth, nq, nqm1, nqp1
2128  DOUBLE PRECISION agamq, elco, fnq, fnqm1, pc, pint, ragq,
2129  1 rq1fac, rqfac, tesco, tsign, xpin
2130  dimension elco(13,12),tesco(3,12)
2131 C ------------------------------------------------------------------
2132 C DCFOD IS CALLED BY THE INTEGRATOR ROUTINE TO SET COEFFICIENTS
2133 C NEEDED THERE. THE COEFFICIENTS FOR THE CURRENT METHOD, AS
2134 C GIVEN BY THE VALUE OF METH, ARE SET FOR ALL ORDERS AND SAVED.
2135 C THE MAXIMUM ORDER ASSUMED HERE IS 12 IF METH = 1 AND 5 IF METH =
2136 C 2. (A SMALLER VALUE OF THE MAXIMUM ORDER IS ALSO ALLOWED.)
2137 C DCFOD IS CALLED ONCE AT THE BEGINNING OF THE PROBLEM,
2138 C AND IS NOT CALLED AGAIN UNLESS AND UNTIL METH IS CHANGED.
2139 C
2140 C THE ELCO ARRAY CONTAINS THE BASIC METHOD COEFFICIENTS.
2141 C THE COEFFICIENTS EL(I), 1 .LE. I .LE. NQ+1, FOR THE METHOD OF
2142 C ORDER NQ ARE STORED IN ELCO(I,NQ). THEY ARE GIVEN BY A
2143 C GENERATING POLYNOMIAL, I.E.,
2144 C L(X) = EL(1) + EL(2)*X + ... + EL(NQ+1)*X**NQ.
2145 C FOR THE IMPLICIT ADAMS METHODS, L(X) IS GIVEN BY
2146 C DL/DX = (X+1)*(X+2)*...*(X+NQ-1)/FACTORIAL(NQ-1), L(-1) =
2147 C 0. FOR THE BDF METHODS, L(X) IS GIVEN BY
2148 C L(X) = (X+1)*(X+2)* ... *(X+NQ)/K,
2149 C WHERE K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ).
2150 C
2151 C THE TESCO ARRAY CONTAINS TEST CONSTANTS USED FOR THE
2152 C LOCAL ERROR TEST AND THE SELECTION OF STEP SIZE AND/OR ORDER.
2153 C AT ORDER NQ, TESCO(K,NQ) IS USED FOR THE SELECTION OF STEP
2154 C SIZE AT ORDER NQ - 1 IF K = 1, AT ORDER NQ IF K = 2, AND AT ORDER
2155 C NQ + 1 IF K = 3.
2156 C ------------------------------------------------------------------
2157  dimension pc(12)
2158 C
2159 C***FIRST EXECUTABLE STATEMENT DCFOD
2160  GO TO (10,60), meth
2161 C
2162  10 CONTINUE
2163  elco(1,1) = 1.0d0
2164  elco(2,1) = 1.0d0
2165  tesco(1,1) = 0.0d0
2166  tesco(2,1) = 2.0d0
2167  tesco(1,2) = 1.0d0
2168  tesco(3,12) = 0.0d0
2169  pc(1) = 1.0d0
2170  rqfac = 1.0d0
2171  DO 50 nq = 2, 12
2172 C ------------------------------------------------------------
2173 C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE
2174 C POLYNOMIAL P(X) = (X+1)*(X+2)*...*(X+NQ-1).
2175 C INITIALLY, P(X) = 1.
2176 C ------------------------------------------------------------
2177  rq1fac = rqfac
2178  rqfac = rqfac/nq
2179  nqm1 = nq - 1
2180  fnqm1 = nqm1
2181  nqp1 = nq + 1
2182 C FORM COEFFICIENTS OF P(X)*(X+NQ-1).
2183 C ----------------------------------
2184  pc(nq) = 0.0d0
2185  DO 20 ib = 1, nqm1
2186  i = nqp1 - ib
2187  pc(i) = pc(i-1) + fnqm1*pc(i)
2188  20 CONTINUE
2189  pc(1) = fnqm1*pc(1)
2190 C COMPUTE INTEGRAL, -1 TO 0, OF P(X) AND X*P(X).
2191 C -----------------------
2192  pint = pc(1)
2193  xpin = pc(1)/2.0d0
2194  tsign = 1.0d0
2195  DO 30 i = 2, nq
2196  tsign = -tsign
2197  pint = pint + tsign*pc(i)/i
2198  xpin = xpin + tsign*pc(i)/(i+1)
2199  30 CONTINUE
2200 C STORE COEFFICIENTS IN ELCO AND TESCO.
2201 C --------------------------------
2202  elco(1,nq) = pint*rq1fac
2203  elco(2,nq) = 1.0d0
2204  DO 40 i = 2, nq
2205  elco(i+1,nq) = rq1fac*pc(i)/i
2206  40 CONTINUE
2207  agamq = rqfac*xpin
2208  ragq = 1.0d0/agamq
2209  tesco(2,nq) = ragq
2210  IF (nq .LT. 12) tesco(1,nqp1) = ragq*rqfac/nqp1
2211  tesco(3,nqm1) = ragq
2212  50 CONTINUE
2213  GO TO 100
2214 C
2215  60 CONTINUE
2216  pc(1) = 1.0d0
2217  rq1fac = 1.0d0
2218  DO 90 nq = 1, 5
2219 C ------------------------------------------------------------
2220 C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE
2221 C POLYNOMIAL P(X) = (X+1)*(X+2)*...*(X+NQ).
2222 C INITIALLY, P(X) = 1.
2223 C ------------------------------------------------------------
2224  fnq = nq
2225  nqp1 = nq + 1
2226 C FORM COEFFICIENTS OF P(X)*(X+NQ).
2227 C ------------------------------------
2228  pc(nqp1) = 0.0d0
2229  DO 70 ib = 1, nq
2230  i = nq + 2 - ib
2231  pc(i) = pc(i-1) + fnq*pc(i)
2232  70 CONTINUE
2233  pc(1) = fnq*pc(1)
2234 C STORE COEFFICIENTS IN ELCO AND TESCO.
2235 C --------------------------------
2236  DO 80 i = 1, nqp1
2237  elco(i,nq) = pc(i)/pc(2)
2238  80 CONTINUE
2239  elco(2,nq) = 1.0d0
2240  tesco(1,nq) = rq1fac
2241  tesco(2,nq) = nqp1/elco(1,nq)
2242  tesco(3,nq) = (nq+2)/elco(1,nq)
2243  rq1fac = rq1fac/fnq
2244  90 CONTINUE
2245  100 CONTINUE
2246  RETURN
2247 C ----------------------- END OF SUBROUTINE DCFOD
2248 C -----------------------

◆ ddebdf()

subroutine ddebdf ( external  DF,
integer  NEQ,
double precision  T,
double precision, dimension(*)  Y,
double precision  TOUT,
integer, dimension(15)  INFO,
double precision, dimension(*)  RTOL,
double precision, dimension(*)  ATOL,
integer  IDID,
double precision, dimension(*)  RWORK,
integer  LRW,
integer, dimension(*)  IWORK,
integer  LIW,
double precision, dimension(*)  RPAR,
integer, dimension(*)  IPAR,
external  DJAC 
)
7 C***BEGIN PROLOGUE DDEBDF
8 C***PURPOSE Solve an initial value problem in ordinary differential
9 C equations using backward differentiation formulas. It is
10 C intended primarily for stiff problems.
11 C***LIBRARY SLATEC (DEPAC)
12 C***CATEGORY I1A2
13 C***TYPE DOUBLE PRECISION (DEBDF-S, DDEBDF-D)
14 C***KEYWORDS BACKWARD DIFFERENTIATION FORMULAS, DEPAC,
15 C INITIAL VALUE PROBLEMS, ODE,
16 C ORDINARY DIFFERENTIAL EQUATIONS, STIFF
17 C***AUTHOR Shampine, L. F., (SNLA)
18 C Watts, H. A., (SNLA)
19 C***DESCRIPTION
20 C
21 C This is the backward differentiation code in the package of
22 C differential equation solvers DEPAC, consisting of the codes
23 C DDERKF, DDEABM, and DDEBDF. Design of the package was by
24 C L. F. Shampine and H. A. Watts. It is documented in
25 C SAND-79-2374 , DEPAC - Design of a User Oriented Package of ODE
26 C Solvers.
27 C DDEBDF is a driver for a modification of the code LSODE written by
28 C A. C. Hindmarsh
29 C Lawrence Livermore Laboratory
30 C Livermore, California 94550
31 C
32 C **********************************************************************
33 C ** DEPAC PACKAGE OVERVIEW **
34 C **********************************************************************
35 C
36 C You have a choice of three differential equation solvers from
37 C DEPAC. The following brief descriptions are meant to aid you
38 C in choosing the most appropriate code for your problem.
39 C
40 C DDERKF is a fifth order Runge-Kutta code. It is the simplest of
41 C the three choices, both algorithmically and in the use of the
42 C code. DDERKF is primarily designed to solve non-stiff and mild-
43 C ly stiff differential equations when derivative evaluations are
44 C not expensive. It should generally not be used to get high
45 C accuracy results nor answers at a great many specific points.
46 C Because DDERKF has very low overhead costs, it will usually
47 C result in the least expensive integration when solving
48 C problems requiring a modest amount of accuracy and having
49 C equations that are not costly to evaluate. DDERKF attempts to
50 C discover when it is not suitable for the task posed.
51 C
52 C DDEABM is a variable order (one through twelve) Adams code. Its
53 C complexity lies somewhere between that of DDERKF and DDEBDF.
54 C DDEABM is primarily designed to solve non-stiff and mildly
55 C stiff differential equations when derivative evaluations are
56 C expensive, high accuracy results are needed or answers at
57 C many specific points are required. DDEABM attempts to discover
58 C when it is not suitable for the task posed.
59 C
60 C DDEBDF is a variable order (one through five) backward
61 C differentiation formula code. It is the most complicated of
62 C the three choices. DDEBDF is primarily designed to solve stiff
63 C differential equations at crude to moderate tolerances.
64 C If the problem is very stiff at all, DDERKF and DDEABM will be
65 C quite inefficient compared to DDEBDF. However, DDEBDF will be
66 C inefficient compared to DDERKF and DDEABM on non-stiff problems
67 C because it uses much more storage, has a much larger overhead,
68 C and the low order formulas will not give high accuracies
69 C efficiently.
70 C
71 C The concept of stiffness cannot be described in a few words.
72 C If you do not know the problem to be stiff, try either DDERKF
73 C or DDEABM. Both of these codes will inform you of stiffness
74 C when the cost of solving such problems becomes important.
75 C
76 C **********************************************************************
77 C ** ABSTRACT **
78 C **********************************************************************
79 C
80 C Subroutine DDEBDF uses the backward differentiation formulas of
81 C orders one through five to integrate a system of NEQ first order
82 C ordinary differential equations of the form
83 C DU/DX = DF(X,U)
84 C when the vector Y(*) of initial values for U(*) at X=T is given.
85 C The subroutine integrates from T to TOUT. It is easy to continue the
86 C integration to get results at additional TOUT. This is the interval
87 C mode of operation. It is also easy for the routine to return with
88 C the solution at each intermediate step on the way to TOUT. This is
89 C the intermediate-output mode of operation.
90 C
91 C **********************************************************************
92 C * Description of The Arguments To DDEBDF (An Overview) *
93 C **********************************************************************
94 C
95 C The Parameters are:
96 C
97 C DF -- This is the name of a subroutine which you provide to
98 C define the differential equations.
99 C
100 C NEQ -- This is the number of (first order) differential
101 C equations to be integrated.
102 C
103 C T -- This is a DOUBLE PRECISION value of the independent
104 C variable.
105 C
106 C Y(*) -- This DOUBLE PRECISION array contains the solution
107 C components at T.
108 C
109 C TOUT -- This is a DOUBLE PRECISION point at which a solution is
110 C desired.
111 C
112 C INFO(*) -- The basic task of the code is to integrate the
113 C differential equations from T to TOUT and return an
114 C answer at TOUT. INFO(*) is an INTEGER array which is used
115 C to communicate exactly how you want this task to be
116 C carried out.
117 C
118 C RTOL, ATOL -- These DOUBLE PRECISION quantities
119 C represent relative and absolute error tolerances which you
120 C provide to indicate how accurately you wish the solution
121 C to be computed. You may choose them to be both scalars
122 C or else both vectors.
123 C
124 C IDID -- This scalar quantity is an indicator reporting what
125 C the code did. You must monitor this INTEGER variable to
126 C decide what action to take next.
127 C
128 C RWORK(*), LRW -- RWORK(*) is a DOUBLE PRECISION work array of
129 C length LRW which provides the code with needed storage
130 C space.
131 C
132 C IWORK(*), LIW -- IWORK(*) is an INTEGER work array of length LIW
133 C which provides the code with needed storage space and an
134 C across call flag.
135 C
136 C RPAR, IPAR -- These are DOUBLE PRECISION and INTEGER parameter
137 C arrays which you can use for communication between your
138 C calling program and the DF subroutine (and the DJAC
139 C subroutine).
140 C
141 C DJAC -- This is the name of a subroutine which you may choose to
142 C provide for defining the Jacobian matrix of partial
143 C derivatives DF/DU.
144 C
145 C Quantities which are used as input items are
146 C NEQ, T, Y(*), TOUT, INFO(*),
147 C RTOL, ATOL, RWORK(1), LRW,
148 C IWORK(1), IWORK(2), and LIW.
149 C
150 C Quantities which may be altered by the code are
151 C T, Y(*), INFO(1), RTOL, ATOL,
152 C IDID, RWORK(*) and IWORK(*).
153 C
154 C **********************************************************************
155 C * INPUT -- What To Do On The First Call To DDEBDF *
156 C **********************************************************************
157 C
158 C The first call of the code is defined to be the start of each new
159 C problem. Read through the descriptions of all the following items,
160 C provide sufficient storage space for designated arrays, set
161 C appropriate variables for the initialization of the problem, and
162 C give information about how you want the problem to be solved.
163 C
164 C
165 C DF -- Provide a subroutine of the form
166 C DF(X,U,UPRIME,RPAR,IPAR)
167 C to define the system of first order differential equations
168 C which is to be solved. For the given values of X and the
169 C vector U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must
170 C evaluate the NEQ components of the system of differential
171 C equations DU/DX=DF(X,U) and store the derivatives in the
172 C array UPRIME(*), that is, UPRIME(I) = * DU(I)/DX * for
173 C equations I=1,...,NEQ.
174 C
175 C Subroutine DF must not alter X or U(*). You must declare
176 C the name DF in an external statement in your program that
177 C calls DDEBDF. You must dimension U and UPRIME in DF.
178 C
179 C RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter
180 C arrays which you can use for communication between your
181 C calling program and subroutine DF. They are not used or
182 C altered by DDEBDF. If you do not need RPAR or IPAR,
183 C ignore these parameters by treating them as dummy
184 C arguments. If you do choose to use them, dimension them in
185 C your calling program and in DF as arrays of appropriate
186 C length.
187 C
188 C NEQ -- Set it to the number of differential equations.
189 C (NEQ .GE. 1)
190 C
191 C T -- Set it to the initial point of the integration.
192 C You must use a program variable for T because the code
193 C changes its value.
194 C
195 C Y(*) -- Set this vector to the initial values of the NEQ solution
196 C components at the initial point. You must dimension Y at
197 C least NEQ in your calling program.
198 C
199 C TOUT -- Set it to the first point at which a solution is desired.
200 C You can take TOUT = T, in which case the code
201 C will evaluate the derivative of the solution at T and
202 C return. Integration either forward in T (TOUT .GT. T)
203 C or backward in T (TOUT .LT. T) is permitted.
204 C
205 C The code advances the solution from T to TOUT using
206 C step sizes which are automatically selected so as to
207 C achieve the desired accuracy. If you wish, the code will
208 C return with the solution and its derivative following
209 C each intermediate step (intermediate-output mode) so that
210 C you can monitor them, but you still must provide TOUT in
211 C accord with the basic aim of the code.
212 C
213 C The first step taken by the code is a critical one
214 C because it must reflect how fast the solution changes near
215 C the initial point. The code automatically selects an
216 C initial step size which is practically always suitable for
217 C the problem. By using the fact that the code will not
218 C step past TOUT in the first step, you could, if necessary,
219 C restrict the length of the initial step size.
220 C
221 C For some problems it may not be permissible to integrate
222 C past a point TSTOP because a discontinuity occurs there
223 C or the solution or its derivative is not defined beyond
224 C TSTOP. When you have declared a TSTOP point (see INFO(4)
225 C and RWORK(1)), you have told the code not to integrate
226 C past TSTOP. In this case any TOUT beyond TSTOP is invalid
227 C input.
228 C
229 C INFO(*) -- Use the INFO array to give the code more details about
230 C how you want your problem solved. This array should be
231 C dimensioned of length 15 to accommodate other members of
232 C DEPAC or possible future extensions, though DDEBDF uses
233 C only the first six entries. You must respond to all of
234 C the following items which are arranged as questions. The
235 C simplest use of the code corresponds to answering all
236 C questions as YES ,i.e. setting all entries of INFO to 0.
237 C
238 C INFO(1) -- This parameter enables the code to initialize
239 C itself. You must set it to indicate the start of every
240 C new problem.
241 C
242 C **** Is this the first call for this problem ...
243 C YES -- Set INFO(1) = 0
244 C NO -- Not applicable here.
245 C See below for continuation calls. ****
246 C
247 C INFO(2) -- How much accuracy you want of your solution
248 C is specified by the error tolerances RTOL and ATOL.
249 C The simplest use is to take them both to be scalars.
250 C To obtain more flexibility, they can both be vectors.
251 C The code must be told your choice.
252 C
253 C **** Are both error tolerances RTOL, ATOL scalars ...
254 C YES -- Set INFO(2) = 0
255 C and input scalars for both RTOL and ATOL
256 C NO -- Set INFO(2) = 1
257 C and input arrays for both RTOL and ATOL ****
258 C
259 C INFO(3) -- The code integrates from T in the direction
260 C of TOUT by steps. If you wish, it will return the
261 C computed solution and derivative at the next
262 C intermediate step (the intermediate-output mode) or
263 C TOUT, whichever comes first. This is a good way to
264 C proceed if you want to see the behavior of the solution.
265 C If you must have solutions at a great many specific
266 C TOUT points, this code will compute them efficiently.
267 C
268 C **** Do you want the solution only at
269 C TOUT (and NOT at the next intermediate step) ...
270 C YES -- Set INFO(3) = 0
271 C NO -- Set INFO(3) = 1 ****
272 C
273 C INFO(4) -- To handle solutions at a great many specific
274 C values TOUT efficiently, this code may integrate past
275 C TOUT and interpolate to obtain the result at TOUT.
276 C Sometimes it is not possible to integrate beyond some
277 C point TSTOP because the equation changes there or it is
278 C not defined past TSTOP. Then you must tell the code
279 C not to go past.
280 C
281 C **** Can the integration be carried out without any
282 C restrictions on the independent variable T ...
283 C YES -- Set INFO(4)=0
284 C NO -- Set INFO(4)=1
285 C and define the stopping point TSTOP by
286 C setting RWORK(1)=TSTOP ****
287 C
288 C INFO(5) -- To solve stiff problems it is necessary to use the
289 C Jacobian matrix of partial derivatives of the system
290 C of differential equations. If you do not provide a
291 C subroutine to evaluate it analytically (see the
292 C description of the item DJAC in the call list), it will
293 C be approximated by numerical differencing in this code.
294 C Although it is less trouble for you to have the code
295 C compute partial derivatives by numerical differencing,
296 C the solution will be more reliable if you provide the
297 C derivatives via DJAC. Sometimes numerical differencing
298 C is cheaper than evaluating derivatives in DJAC and
299 C sometimes it is not - this depends on your problem.
300 C
301 C If your problem is linear, i.e. has the form
302 C DU/DX = DF(X,U) = J(X)*U + G(X) for some matrix J(X)
303 C and vector G(X), the Jacobian matrix DF/DU = J(X).
304 C Since you must provide a subroutine to evaluate DF(X,U)
305 C analytically, it is little extra trouble to provide
306 C subroutine DJAC for evaluating J(X) analytically.
307 C Furthermore, in such cases, numerical differencing is
308 C much more expensive than analytic evaluation.
309 C
310 C **** Do you want the code to evaluate the partial
311 C derivatives automatically by numerical differences ...
312 C YES -- Set INFO(5)=0
313 C NO -- Set INFO(5)=1
314 C and provide subroutine DJAC for evaluating the
315 C Jacobian matrix ****
316 C
317 C INFO(6) -- DDEBDF will perform much better if the Jacobian
318 C matrix is banded and the code is told this. In this
319 C case, the storage needed will be greatly reduced,
320 C numerical differencing will be performed more cheaply,
321 C and a number of important algorithms will execute much
322 C faster. The differential equation is said to have
323 C half-bandwidths ML (lower) and MU (upper) if equation I
324 C involves only unknowns Y(J) with
325 C I-ML .LE. J .LE. I+MU
326 C for all I=1,2,...,NEQ. Thus, ML and MU are the widths
327 C of the lower and upper parts of the band, respectively,
328 C with the main diagonal being excluded. If you do not
329 C indicate that the equation has a banded Jacobian,
330 C the code works with a full matrix of NEQ**2 elements
331 C (stored in the conventional way). Computations with
332 C banded matrices cost less time and storage than with
333 C full matrices if 2*ML+MU .LT. NEQ. If you tell the
334 C code that the Jacobian matrix has a banded structure and
335 C you want to provide subroutine DJAC to compute the
336 C partial derivatives, then you must be careful to store
337 C the elements of the Jacobian matrix in the special form
338 C indicated in the description of DJAC.
339 C
340 C **** Do you want to solve the problem using a full
341 C (dense) Jacobian matrix (and not a special banded
342 C structure) ...
343 C YES -- Set INFO(6)=0
344 C NO -- Set INFO(6)=1
345 C and provide the lower (ML) and upper (MU)
346 C bandwidths by setting
347 C IWORK(1)=ML
348 C IWORK(2)=MU ****
349 C
350 C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL)
351 C error tolerances to tell the code how accurately you want
352 C the solution to be computed. They must be defined as
353 C program variables because the code may change them. You
354 C have two choices --
355 C Both RTOL and ATOL are scalars. (INFO(2)=0)
356 C Both RTOL and ATOL are vectors. (INFO(2)=1)
357 C In either case all components must be non-negative.
358 C
359 C The tolerances are used by the code in a local error test
360 C at each step which requires roughly that
361 C ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL
362 C for each vector component.
363 C (More specifically, a root-mean-square norm is used to
364 C measure the size of vectors, and the error test uses the
365 C magnitude of the solution at the beginning of the step.)
366 C
367 C The true (global) error is the difference between the true
368 C solution of the initial value problem and the computed
369 C approximation. Practically all present day codes,
370 C including this one, control the local error at each step
371 C and do not even attempt to control the global error
372 C directly. Roughly speaking, they produce a solution Y(T)
373 C which satisfies the differential equations with a
374 C residual R(T), DY(T)/DT = DF(T,Y(T)) + R(T) ,
375 C and, almost always, R(T) is bounded by the error
376 C tolerances. Usually, but not always, the true accuracy of
377 C the computed Y is comparable to the error tolerances. This
378 C code will usually, but not always, deliver a more accurate
379 C solution if you reduce the tolerances and integrate again.
380 C By comparing two such solutions you can get a fairly
381 C reliable idea of the true error in the solution at the
382 C bigger tolerances.
383 C
384 C Setting ATOL=0. results in a pure relative error test on
385 C that component. Setting RTOL=0. results in a pure abso-
386 C lute error test on that component. A mixed test with non-
387 C zero RTOL and ATOL corresponds roughly to a relative error
388 C test when the solution component is much bigger than ATOL
389 C and to an absolute error test when the solution component
390 C is smaller than the threshold ATOL.
391 C
392 C Proper selection of the absolute error control parameters
393 C ATOL requires you to have some idea of the scale of the
394 C solution components. To acquire this information may mean
395 C that you will have to solve the problem more than once. In
396 C the absence of scale information, you should ask for some
397 C relative accuracy in all the components (by setting RTOL
398 C values non-zero) and perhaps impose extremely small
399 C absolute error tolerances to protect against the danger of
400 C a solution component becoming zero.
401 C
402 C The code will not attempt to compute a solution at an
403 C accuracy unreasonable for the machine being used. It will
404 C advise you if you ask for too much accuracy and inform
405 C you as to the maximum accuracy it believes possible.
406 C
407 C RWORK(*) -- Dimension this DOUBLE PRECISION work array of length
408 C LRW in your calling program.
409 C
410 C RWORK(1) -- If you have set INFO(4)=0, you can ignore this
411 C optional input parameter. Otherwise you must define a
412 C stopping point TSTOP by setting RWORK(1) = TSTOP.
413 C (For some problems it may not be permissible to integrate
414 C past a point TSTOP because a discontinuity occurs there
415 C or the solution or its derivative is not defined beyond
416 C TSTOP.)
417 C
418 C LRW -- Set it to the declared length of the RWORK array.
419 C You must have
420 C LRW .GE. 250+10*NEQ+NEQ**2
421 C for the full (dense) Jacobian case (when INFO(6)=0), or
422 C LRW .GE. 250+10*NEQ+(2*ML+MU+1)*NEQ
423 C for the banded Jacobian case (when INFO(6)=1).
424 C
425 C IWORK(*) -- Dimension this INTEGER work array of length LIW in
426 C your calling program.
427 C
428 C IWORK(1), IWORK(2) -- If you have set INFO(6)=0, you can ignore
429 C these optional input parameters. Otherwise you must define
430 C the half-bandwidths ML (lower) and MU (upper) of the
431 C Jacobian matrix by setting IWORK(1) = ML and
432 C IWORK(2) = MU. (The code will work with a full matrix
433 C of NEQ**2 elements unless it is told that the problem has
434 C a banded Jacobian, in which case the code will work with
435 C a matrix containing at most (2*ML+MU+1)*NEQ elements.)
436 C
437 C LIW -- Set it to the declared length of the IWORK array.
438 C You must have LIW .GE. 56+NEQ.
439 C
440 C RPAR, IPAR -- These are parameter arrays, of DOUBLE PRECISION and
441 C INTEGER type, respectively. You can use them for
442 C communication between your program that calls DDEBDF and
443 C the DF subroutine (and the DJAC subroutine). They are not
444 C used or altered by DDEBDF. If you do not need RPAR or
445 C IPAR, ignore these parameters by treating them as dummy
446 C arguments. If you do choose to use them, dimension them in
447 C your calling program and in DF (and in DJAC) as arrays of
448 C appropriate length.
449 C
450 C DJAC -- If you have set INFO(5)=0, you can ignore this parameter
451 C by treating it as a dummy argument. (For some compilers
452 C you may have to write a dummy subroutine named DJAC in
453 C order to avoid problems associated with missing external
454 C routine names.) Otherwise, you must provide a subroutine
455 C of the form
456 C DJAC(X,U,PD,NROWPD,RPAR,IPAR)
457 C to define the Jacobian matrix of partial derivatives DF/DU
458 C of the system of differential equations DU/DX = DF(X,U).
459 C For the given values of X and the vector
460 C U(*)=(U(1),U(2),...,U(NEQ)), the subroutine must evaluate
461 C the non-zero partial derivatives DF(I)/DU(J) for each
462 C differential equation I=1,...,NEQ and each solution
463 C component J=1,...,NEQ , and store these values in the
464 C matrix PD. The elements of PD are set to zero before each
465 C call to DJAC so only non-zero elements need to be defined.
466 C
467 C Subroutine DJAC must not alter X, U(*), or NROWPD. You
468 C must declare the name DJAC in an external statement in
469 C your program that calls DDEBDF. NROWPD is the row
470 C dimension of the PD matrix and is assigned by the code.
471 C Therefore you must dimension PD in DJAC according to
472 C DIMENSION PD(NROWPD,1)
473 C You must also dimension U in DJAC.
474 C
475 C The way you must store the elements into the PD matrix
476 C depends on the structure of the Jacobian which you
477 C indicated by INFO(6).
478 C *** INFO(6)=0 -- Full (Dense) Jacobian ***
479 C When you evaluate the (non-zero) partial derivative
480 C of equation I with respect to variable J, you must
481 C store it in PD according to
482 C PD(I,J) = * DF(I)/DU(J) *
483 C *** INFO(6)=1 -- Banded Jacobian with ML Lower and MU
484 C Upper Diagonal Bands (refer to INFO(6) description of
485 C ML and MU) ***
486 C When you evaluate the (non-zero) partial derivative
487 C of equation I with respect to variable J, you must
488 C store it in PD according to
489 C IROW = I - J + ML + MU + 1
490 C PD(IROW,J) = * DF(I)/DU(J) *
491 C
492 C RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter
493 C arrays which you can use for communication between your
494 C calling program and your Jacobian subroutine DJAC. They
495 C are not altered by DDEBDF. If you do not need RPAR or
496 C IPAR, ignore these parameters by treating them as dummy
497 C arguments. If you do choose to use them, dimension them
498 C in your calling program and in DJAC as arrays of
499 C appropriate length.
500 C
501 C **********************************************************************
502 C * OUTPUT -- After any return from DDEBDF *
503 C **********************************************************************
504 C
505 C The principal aim of the code is to return a computed solution at
506 C TOUT, although it is also possible to obtain intermediate results
507 C along the way. To find out whether the code achieved its goal
508 C or if the integration process was interrupted before the task was
509 C completed, you must check the IDID parameter.
510 C
511 C
512 C T -- The solution was successfully advanced to the
513 C output value of T.
514 C
515 C Y(*) -- Contains the computed solution approximation at T.
516 C You may also be interested in the approximate derivative
517 C of the solution at T. It is contained in
518 C RWORK(21),...,RWORK(20+NEQ).
519 C
520 C IDID -- Reports what the code did
521 C
522 C *** Task Completed ***
523 C Reported by positive values of IDID
524 C
525 C IDID = 1 -- A step was successfully taken in the
526 C intermediate-output mode. The code has not
527 C yet reached TOUT.
528 C
529 C IDID = 2 -- The integration to TOUT was successfully
530 C completed (T=TOUT) by stepping exactly to TOUT.
531 C
532 C IDID = 3 -- The integration to TOUT was successfully
533 C completed (T=TOUT) by stepping past TOUT.
534 C Y(*) is obtained by interpolation.
535 C
536 C *** Task Interrupted ***
537 C Reported by negative values of IDID
538 C
539 C IDID = -1 -- A large amount of work has been expended.
540 C (500 steps attempted)
541 C
542 C IDID = -2 -- The error tolerances are too stringent.
543 C
544 C IDID = -3 -- The local error test cannot be satisfied
545 C because you specified a zero component in ATOL
546 C and the corresponding computed solution
547 C component is zero. Thus, a pure relative error
548 C test is impossible for this component.
549 C
550 C IDID = -4,-5 -- Not applicable for this code but used
551 C by other members of DEPAC.
552 C
553 C IDID = -6 -- DDEBDF had repeated convergence test failures
554 C on the last attempted step.
555 C
556 C IDID = -7 -- DDEBDF had repeated error test failures on
557 C the last attempted step.
558 C
559 C IDID = -8,..,-32 -- Not applicable for this code but
560 C used by other members of DEPAC or possible
561 C future extensions.
562 C
563 C *** Task Terminated ***
564 C Reported by the value of IDID=-33
565 C
566 C IDID = -33 -- The code has encountered trouble from which
567 C it cannot recover. A message is printed
568 C explaining the trouble and control is returned
569 C to the calling program. For example, this
570 C occurs when invalid input is detected.
571 C
572 C RTOL, ATOL -- These quantities remain unchanged except when
573 C IDID = -2. In this case, the error tolerances have been
574 C increased by the code to values which are estimated to be
575 C appropriate for continuing the integration. However, the
576 C reported solution at T was obtained using the input values
577 C of RTOL and ATOL.
578 C
579 C RWORK, IWORK -- Contain information which is usually of no
580 C interest to the user but necessary for subsequent calls.
581 C However, you may find use for
582 C
583 C RWORK(11)--which contains the step size H to be
584 C attempted on the next step.
585 C
586 C RWORK(12)--If the tolerances have been increased by the
587 C code (IDID = -2) , they were multiplied by the
588 C value in RWORK(12).
589 C
590 C RWORK(13)--which contains the current value of the
591 C independent variable, i.e. the farthest point
592 C integration has reached. This will be
593 C different from T only when interpolation has
594 C been performed (IDID=3).
595 C
596 C RWORK(20+I)--which contains the approximate derivative
597 C of the solution component Y(I). In DDEBDF, it
598 C is never obtained by calling subroutine DF to
599 C evaluate the differential equation using T and
600 C Y(*), except at the initial point of
601 C integration.
602 C
603 C **********************************************************************
604 C ** INPUT -- What To Do To Continue The Integration **
605 C ** (calls after the first) **
606 C **********************************************************************
607 C
608 C This code is organized so that subsequent calls to continue the
609 C integration involve little (if any) additional effort on your
610 C part. You must monitor the IDID parameter in order to determine
611 C what to do next.
612 C
613 C Recalling that the principal task of the code is to integrate
614 C from T to TOUT (the interval mode), usually all you will need
615 C to do is specify a new TOUT upon reaching the current TOUT.
616 C
617 C Do not alter any quantity not specifically permitted below,
618 C in particular do not alter NEQ, T, Y(*), RWORK(*), IWORK(*) or
619 C the differential equation in subroutine DF. Any such alteration
620 C constitutes a new problem and must be treated as such, i.e.
621 C you must start afresh.
622 C
623 C You cannot change from vector to scalar error control or vice
624 C versa (INFO(2)) but you can change the size of the entries of
625 C RTOL, ATOL. Increasing a tolerance makes the equation easier
626 C to integrate. Decreasing a tolerance will make the equation
627 C harder to integrate and should generally be avoided.
628 C
629 C You can switch from the intermediate-output mode to the
630 C interval mode (INFO(3)) or vice versa at any time.
631 C
632 C If it has been necessary to prevent the integration from going
633 C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the
634 C code will not integrate to any TOUT beyond the currently
635 C specified TSTOP. Once TSTOP has been reached you must change
636 C the value of TSTOP or set INFO(4)=0. You may change INFO(4)
637 C or TSTOP at any time but you must supply the value of TSTOP in
638 C RWORK(1) whenever you set INFO(4)=1.
639 C
640 C Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2)
641 C unless you are going to restart the code.
642 C
643 C The parameter INFO(1) is used by the code to indicate the
644 C beginning of a new problem and to indicate whether integration
645 C is to be continued. You must input the value INFO(1) = 0
646 C when starting a new problem. You must input the value
647 C INFO(1) = 1 if you wish to continue after an interrupted task.
648 C Do not set INFO(1) = 0 on a continuation call unless you
649 C want the code to restart at the current T.
650 C
651 C *** Following a Completed Task ***
652 C If
653 C IDID = 1, call the code again to continue the integration
654 C another step in the direction of TOUT.
655 C
656 C IDID = 2 or 3, define a new TOUT and call the code again.
657 C TOUT must be different from T. You cannot change
658 C the direction of integration without restarting.
659 C
660 C *** Following an Interrupted Task ***
661 C To show the code that you realize the task was
662 C interrupted and that you want to continue, you
663 C must take appropriate action and reset INFO(1) = 1
664 C If
665 C IDID = -1, the code has attempted 500 steps.
666 C If you want to continue, set INFO(1) = 1 and
667 C call the code again. An additional 500 steps
668 C will be allowed.
669 C
670 C IDID = -2, the error tolerances RTOL, ATOL have been
671 C increased to values the code estimates appropriate
672 C for continuing. You may want to change them
673 C yourself. If you are sure you want to continue
674 C with relaxed error tolerances, set INFO(1)=1 and
675 C call the code again.
676 C
677 C IDID = -3, a solution component is zero and you set the
678 C corresponding component of ATOL to zero. If you
679 C are sure you want to continue, you must first
680 C alter the error criterion to use positive values
681 C for those components of ATOL corresponding to zero
682 C solution components, then set INFO(1)=1 and call
683 C the code again.
684 C
685 C IDID = -4,-5 --- cannot occur with this code but used
686 C by other members of DEPAC.
687 C
688 C IDID = -6, repeated convergence test failures occurred
689 C on the last attempted step in DDEBDF. An inaccu-
690 C rate Jacobian may be the problem. If you are
691 C absolutely certain you want to continue, restart
692 C the integration at the current T by setting
693 C INFO(1)=0 and call the code again.
694 C
695 C IDID = -7, repeated error test failures occurred on the
696 C last attempted step in DDEBDF. A singularity in
697 C the solution may be present. You should re-
698 C examine the problem being solved. If you are
699 C absolutely certain you want to continue, restart
700 C the integration at the current T by setting
701 C INFO(1)=0 and call the code again.
702 C
703 C IDID = -8,..,-32 --- cannot occur with this code but
704 C used by other members of DDEPAC or possible future
705 C extensions.
706 C
707 C *** Following a Terminated Task ***
708 C If
709 C IDID = -33, you cannot continue the solution of this
710 C problem. An attempt to do so will result in your
711 C run being terminated.
712 C
713 C **********************************************************************
714 C
715 C ***** Warning *****
716 C
717 C If DDEBDF is to be used in an overlay situation, you must save and
718 C restore certain items used internally by DDEBDF (values in the
719 C common block DDEBD1). This can be accomplished as follows.
720 C
721 C To save the necessary values upon return from DDEBDF, simply call
722 C DSVCO(RWORK(22+NEQ),IWORK(21+NEQ)).
723 C
724 C To restore the necessary values before the next call to DDEBDF,
725 C simply call DRSCO(RWORK(22+NEQ),IWORK(21+NEQ)).
726 C
727 C***REFERENCES L. F. Shampine and H. A. Watts, DEPAC - design of a user
728 C oriented package of ODE solvers, Report SAND79-2374,
729 C Sandia Laboratories, 1979.
730 C***ROUTINES CALLED DLSOD, XERMSG
731 C***COMMON BLOCKS DDEBD1
732 C***REVISION HISTORY (YYMMDD)
733 C 820301 DATE WRITTEN
734 C 890831 Modified array declarations. (WRB)
735 C 891024 Changed references from DVNORM to DHVNRM. (WRB)
736 C 891024 REVISION DATE from Version 3.2
737 C 891214 Prologue converted to Version 4.0 format. (BAB)
738 C 900326 Removed duplicate information from DESCRIPTION section.
739 C (WRB)
740 C 900510 Convert XERRWV calls to XERMSG calls, make Prologue comments
741 C consistent with DEBDF. (RWC)
742 C 920501 Reformatted the REFERENCES section. (WRB)
743 C***END PROLOGUE DDEBDF
744  INTEGER iacor, iband, ibegin, icomi, icomr, idelsn, idid, ier,
745  1 iewt, iinout, iinteg, ijac, ilrw, info, init,
746  2 iowns, ipar, iquit, isavf, itol, itstar, itstop, iwm,
747  3 iwork, iyh, iypout, jstart, kflag, ksteps, l, liw, lrw,
748  4 maxord, meth, miter, ml, mu, n, neq, nfe, nje, nq, nqu,
749  5 nst
750  DOUBLE PRECISION atol, el0, h, hmin, hmxi, hu, rowns, rpar,
751  1 rtol, rwork, t, tn, told, tout, uround, y
752  LOGICAL intout
753  CHARACTER*8 xern1, xern2
754  CHARACTER*16 xern3
755 C
756  dimension y(*),info(15),rtol(*),atol(*),rwork(*),iwork(*),
757  1 rpar(*),ipar(*)
758 C
759  COMMON /ddebd1/ told,rowns(210),el0,h,hmin,hmxi,hu,tn,uround,
760  1 iquit,init,iyh,iewt,iacor,isavf,iwm,ksteps,ibegin,
761  2 itol,iinteg,itstop,ijac,iband,iowns(6),ier,jstart,
762  3 kflag,l,meth,miter,maxord,n,nq,nst,nfe,nje,nqu
763 C
764  EXTERNAL df, djac
765 C
766 C CHECK FOR AN APPARENT INFINITE LOOP
767 C
768 C***FIRST EXECUTABLE STATEMENT DDEBDF
769  IF (info(1) .EQ. 0) iwork(liw) = 0
770 C
771  IF (iwork(liw).GE. 5) THEN
772  IF (t .EQ. rwork(21+neq)) THEN
773  WRITE (xern3, '(1PE15.6)') t
774  CALL xermsg ('SLATEC', 'DDEBDF',
775  * 'AN APPARENT INFINITE LOOP HAS BEEN DETECTED.$$' //
776  * 'YOU HAVE MADE REPEATED CALLS AT T = ' // xern3 //
777  * ' AND THE INTEGRATION HAS NOT ADVANCED. CHECK THE ' //
778  * 'WAY YOU HAVE SET PARAMETERS FOR THE CALL TO THE ' //
779  * 'CODE, PARTICULARLY INFO(1).', 13, 2)
780  RETURN
781  ENDIF
782  ENDIF
783 C
784  idid = 0
785 C
786 C CHECK VALIDITY OF INFO PARAMETERS
787 C
788  IF (info(1) .NE. 0 .AND. info(1) .NE. 1) THEN
789  WRITE (xern1, '(I8)') info(1)
790  CALL xermsg ('SLATEC', 'DDEBDF', 'INFO(1) MUST BE SET TO 0 ' //
791  * 'FOR THE START OF A NEW PROBLEM, AND MUST BE SET TO 1 ' //
792  * 'FOLLOWING AN INTERRUPTED TASK. YOU ARE ATTEMPTING TO ' //
793  * 'CONTINUE THE INTEGRATION ILLEGALLY BY CALLING THE ' //
794  * 'CODE WITH INFO(1) = ' // xern1, 3, 1)
795  idid = -33
796  ENDIF
797 C
798  IF (info(2) .NE. 0 .AND. info(2) .NE. 1) THEN
799  WRITE (xern1, '(I8)') info(2)
800  CALL xermsg ('SLATEC', 'DDEBDF', 'INFO(2) MUST BE 0 OR 1 ' //
801  * 'INDICATING SCALAR AND VECTOR ERROR TOLERANCES, ' //
802  * 'RESPECTIVELY. YOU HAVE CALLED THE CODE WITH INFO(2) = ' //
803  * xern1, 4, 1)
804  idid = -33
805  ENDIF
806 C
807  IF (info(3) .NE. 0 .AND. info(3) .NE. 1) THEN
808  WRITE (xern1, '(I8)') info(3)
809  CALL xermsg ('SLATEC', 'DDEBDF', 'INFO(3) MUST BE 0 OR 1 ' //
810  * 'INDICATING THE INTERVAL OR INTERMEDIATE-OUTPUT MODE OF ' //
811  * 'INTEGRATION, RESPECTIVELY. YOU HAVE CALLED THE CODE ' //
812  * 'WITH INFO(3) = ' // xern1, 5, 1)
813  idid = -33
814  ENDIF
815 C
816  IF (info(4) .NE. 0 .AND. info(4) .NE. 1) THEN
817  WRITE (xern1, '(I8)') info(4)
818  CALL xermsg ('SLATEC', 'DDEBDF', 'INFO(4) MUST BE 0 OR 1 ' //
819  * 'INDICATING WHETHER OR NOT THE INTEGRATION INTERVAL IS ' //
820  * 'TO BE RESTRICTED BY A POINT TSTOP. YOU HAVE CALLED ' //
821  * 'THE CODE WITH INFO(4) = ' // xern1, 14, 1)
822  idid = -33
823  ENDIF
824 C
825  IF (info(5) .NE. 0 .AND. info(5) .NE. 1) THEN
826  WRITE (xern1, '(I8)') info(5)
827  CALL xermsg ('SLATEC', 'DDEBDF', 'INFO(5) MUST BE 0 OR 1 ' //
828  * 'INDICATING WHETHER THE CODE IS TOLD TO FORM THE ' //
829  * 'JACOBIAN MATRIX BY NUMERICAL DIFFERENCING OR YOU ' //
830  * 'PROVIDE A SUBROUTINE TO EVALUATE IT ANALYTICALLY. ' //
831  * 'YOU HAVE CALLED THE CODE WITH INFO(5) = ' // xern1, 15, 1)
832  idid = -33
833  ENDIF
834 C
835  IF (info(6) .NE. 0 .AND. info(6) .NE. 1) THEN
836  WRITE (xern1, '(I8)') info(6)
837  CALL xermsg ('SLATEC', 'DDEBDF', 'INFO(6) MUST BE 0 OR 1 ' //
838  * 'INDICATING WHETHER THE CODE IS TOLD TO TREAT THE ' //
839  * 'JACOBIAN AS A FULL (DENSE) MATRIX OR AS HAVING A ' //
840  * 'SPECIAL BANDED STRUCTURE. YOU HAVE CALLED THE CODE ' //
841  * 'WITH INFO(6) = ' // xern1, 16, 1)
842  idid = -33
843  ENDIF
844 C
845  ilrw = neq
846  IF (info(6) .NE. 0) THEN
847 C
848 C CHECK BANDWIDTH PARAMETERS
849 C
850  ml = iwork(1)
851  mu = iwork(2)
852  ilrw = 2*ml + mu + 1
853 C
854  IF (ml.LT.0 .OR. ml.GE.neq .OR. mu.LT.0 .OR. mu.GE.neq) THEN
855  WRITE (xern1, '(I8)') ml
856  WRITE (xern2, '(I8)') mu
857  CALL xermsg ('SLATEC', 'DDEBDF', 'YOU HAVE SET INFO(6) ' //
858  * '= 1, TELLING THE CODE THAT THE JACOBIAN MATRIX HAS ' //
859  * 'A SPECIAL BANDED STRUCTURE. HOWEVER, THE LOWER ' //
860  * '(UPPER) BANDWIDTHS ML (MU) VIOLATE THE CONSTRAINTS ' //
861  * .GE..LT.'ML,MU 0 AND ML,MU NEQ. YOU HAVE CALLED ' //
862  * 'THE CODE WITH ML = ' // xern1 // ' AND MU = ' // xern2,
863  * 17, 1)
864  idid = -33
865  ENDIF
866  ENDIF
867 C
868 C CHECK LRW AND LIW FOR SUFFICIENT STORAGE ALLOCATION
869 C
870  IF (lrw .LT. 250 + (10 + ilrw)*neq) THEN
871  WRITE (xern1, '(I8)') lrw
872  IF (info(6) .EQ. 0) THEN
873  CALL xermsg ('SLATEC', 'DDEBDF', 'LENGTH OF ARRAY RWORK ' //
874  * 'MUST BE AT LEAST 250 + 10*NEQ + NEQ*NEQ.$$' //
875  * 'YOU HAVE CALLED THE CODE WITH LRW = ' // xern1, 1, 1)
876  ELSE
877  CALL xermsg ('SLATEC', 'DDEBDF', 'LENGTH OF ARRAY RWORK ' //
878  * 'MUST BE AT LEAST 250 + 10*NEQ + (2*ML+MU+1)*NEQ.$$' //
879  * 'YOU HAVE CALLED THE CODE WITH LRW = ' // xern1, 18, 1)
880  ENDIF
881  idid = -33
882  ENDIF
883 C
884  IF (liw .LT. 56 + neq) THEN
885  WRITE (xern1, '(I8)') liw
886  CALL xermsg ('SLATEC', 'DDEBDF', 'LENGTH OF ARRAY IWORK ' //
887  * 'BE AT LEAST 56 + NEQ. YOU HAVE CALLED THE CODE WITH ' //
888  * 'LIW = ' // xern1, 2, 1)
889  idid = -33
890  ENDIF
891 C
892 C COMPUTE THE INDICES FOR THE ARRAYS TO BE STORED IN THE WORK
893 C ARRAY AND RESTORE COMMON BLOCK DATA
894 C
895  icomi = 21 + neq
896  iinout = icomi + 33
897 C
898  iypout = 21
899  itstar = 21 + neq
900  icomr = 22 + neq
901 C
902  IF (info(1) .NE. 0) intout = iwork(iinout) .NE. (-1)
903 C CALL DRSCO(RWORK(ICOMR),IWORK(ICOMI))
904 C
905  iyh = icomr + 218
906  iewt = iyh + 6*neq
907  isavf = iewt + neq
908  iacor = isavf + neq
909  iwm = iacor + neq
910  idelsn = iwm + 2 + ilrw*neq
911 C
912  ibegin = info(1)
913  itol = info(2)
914  iinteg = info(3)
915  itstop = info(4)
916  ijac = info(5)
917  iband = info(6)
918  rwork(itstar) = t
919 C
920  CALL dlsod(df,neq,t,y,tout,rtol,atol,idid,rwork(iypout),
921  1 rwork(iyh),rwork(iyh),rwork(iewt),rwork(isavf),
922  2 rwork(iacor),rwork(iwm),iwork(1),djac,intout,
923  3 rwork(1),rwork(12),rwork(idelsn),rpar,ipar)
924 C
925  iwork(iinout) = -1
926  IF (intout) iwork(iinout) = 1
927 C
928  IF (idid .NE. (-2)) iwork(liw) = iwork(liw) + 1
929  IF (t .NE. rwork(itstar)) iwork(liw) = 0
930 C CALL DSVCO(RWORK(ICOMR),IWORK(ICOMI))
931  rwork(11) = h
932  rwork(13) = tn
933  info(1) = ibegin
934 C
935  RETURN
init
subroutine init(nktet, inodfa, ipofa, netet_)
Definition: init.f:20
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
xermsg
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: ddeabm.f:1265
djac
subroutine djac(x, u, pd, nrowpd, rpar, nev)
Definition: subspace.f:163
dlsod
subroutine dlsod(DF, NEQ, T, Y, TOUT, RTOL, ATOL, IDID, YPOUT, YH, YH1, EWT, SAVF, ACOR, WM, IWM, DJAC, INTOUT, TSTOP, TOLFAC, DELSGN, RPAR, IPAR)
Definition: ddebdf.f:941

◆ dgbfa()

subroutine dgbfa ( double precision, dimension(lda,*)  ABD,
integer  LDA,
integer  N,
integer  ML,
integer  MU,
integer, dimension(*)  IPVT,
integer  INFO 
)
2730 C***BEGIN PROLOGUE DGBFA
2731 C***PURPOSE Factor a band matrix using Gaussian elimination.
2732 C***LIBRARY SLATEC (LINPACK)
2733 C***CATEGORY D2A2
2734 C***TYPE DOUBLE PRECISION (SGBFA-S, DGBFA-D, CGBFA-C)
2735 C***KEYWORDS BANDED, LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION
2736 C***AUTHOR Moler, C. B., (U. of New Mexico)
2737 C***DESCRIPTION
2738 C
2739 C DGBFA factors a double precision band matrix by elimination.
2740 C
2741 C DGBFA is usually called by DGBCO, but it can be called
2742 C directly with a saving in time if RCOND is not needed.
2743 C
2744 C On Entry
2745 C
2746 C ABD DOUBLE PRECISION(LDA, N)
2747 C contains the matrix in band storage. The columns
2748 C of the matrix are stored in the columns of ABD and
2749 C the diagonals of the matrix are stored in rows
2750 C ML+1 through 2*ML+MU+1 of ABD .
2751 C See the comments below for details.
2752 C
2753 C LDA INTEGER
2754 C the leading dimension of the array ABD .
2755 C LDA must be .GE. 2*ML + MU + 1 .
2756 C
2757 C N INTEGER
2758 C the order of the original matrix.
2759 C
2760 C ML INTEGER
2761 C number of diagonals below the main diagonal.
2762 C 0 .LE. ML .LT. N .
2763 C
2764 C MU INTEGER
2765 C number of diagonals above the main diagonal.
2766 C 0 .LE. MU .LT. N .
2767 C More efficient if ML .LE. MU .
2768 C On Return
2769 C
2770 C ABD an upper triangular matrix in band storage and
2771 C the multipliers which were used to obtain it.
2772 C The factorization can be written A = L*U where
2773 C L is a product of permutation and unit lower
2774 C triangular matrices and U is upper triangular.
2775 C
2776 C IPVT INTEGER(N)
2777 C an integer vector of pivot indices.
2778 C
2779 C INFO INTEGER
2780 C = 0 normal value.
2781 C = K if U(K,K) .EQ. 0.0 . This is not an error
2782 C condition for this subroutine, but it does
2783 C indicate that DGBSL will divide by zero if
2784 C called. Use RCOND in DGBCO for a reliable
2785 C indication of singularity.
2786 C
2787 C Band Storage
2788 C
2789 C If A is a band matrix, the following program segment
2790 C will set up the input.
2791 C
2792 C ML = (band width below the diagonal)
2793 C MU = (band width above the diagonal)
2794 C M = ML + MU + 1
2795 C DO 20 J = 1, N
2796 C I1 = MAX(1, J-MU)
2797 C I2 = MIN(N, J+ML)
2798 C DO 10 I = I1, I2
2799 C K = I - J + M
2800 C ABD(K,J) = A(I,J)
2801 C 10 CONTINUE
2802 C 20 CONTINUE
2803 C
2804 C This uses rows ML+1 through 2*ML+MU+1 of ABD .
2805 C In addition, the first ML rows in ABD are used for
2806 C elements generated during the triangularization.
2807 C The total number of rows needed in ABD is 2*ML+MU+1 .
2808 C The ML+MU by ML+MU upper left triangle and the
2809 C ML by ML lower right triangle are not referenced.
2810 C
2811 C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
2812 C Stewart, LINPACK Users' Guide, SIAM, 1979.
2813 C***ROUTINES CALLED DAXPY, DSCAL, IDAMAX
2814 C***REVISION HISTORY (YYMMDD)
2815 C 780814 DATE WRITTEN
2816 C 890531 Changed all specific intrinsics to generic. (WRB)
2817 C 890831 Modified array declarations. (WRB)
2818 C 890831 REVISION DATE from Version 3.2
2819 C 891214 Prologue converted to Version 4.0 format. (BAB)
2820 C 900326 Removed duplicate information from DESCRIPTION section.
2821 C (WRB)
2822 C 920501 Reformatted the REFERENCES section. (WRB)
2823 C***END PROLOGUE DGBFA
2824  INTEGER lda,n,ml,mu,ipvt(*),info
2825  DOUBLE PRECISION abd(lda,*)
2826 C
2827  DOUBLE PRECISION t
2828  INTEGER i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
2829 C
2830 C***FIRST EXECUTABLE STATEMENT DGBFA
2831  m = ml + mu + 1
2832  info = 0
2833 C
2834 C ZERO INITIAL FILL-IN COLUMNS
2835 C
2836  j0 = mu + 2
2837  j1 = min(n,m) - 1
2838  IF (j1 .LT. j0) GO TO 30
2839  DO 20 jz = j0, j1
2840  i0 = m + 1 - jz
2841  DO 10 i = i0, ml
2842  abd(i,jz) = 0.0d0
2843  10 CONTINUE
2844  20 CONTINUE
2845  30 CONTINUE
2846  jz = j1
2847  ju = 0
2848 C
2849 C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
2850 C
2851  nm1 = n - 1
2852  IF (nm1 .LT. 1) GO TO 130
2853  DO 120 k = 1, nm1
2854  kp1 = k + 1
2855 C
2856 C ZERO NEXT FILL-IN COLUMN
2857 C
2858  jz = jz + 1
2859  IF (jz .GT. n) GO TO 50
2860  IF (ml .LT. 1) GO TO 50
2861  DO 40 i = 1, ml
2862  abd(i,jz) = 0.0d0
2863  40 CONTINUE
2864  50 CONTINUE
2865 C
2866 C FIND L = PIVOT INDEX
2867 C
2868  lm = min(ml,n-k)
2869  l = idamax(lm+1,abd(m,k),1) + m - 1
2870  ipvt(k) = l + k - m
2871 C
2872 C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
2873 C
2874  IF (abd(l,k) .EQ. 0.0d0) GO TO 100
2875 C
2876 C INTERCHANGE IF NECESSARY
2877 C
2878  IF (l .EQ. m) GO TO 60
2879  t = abd(l,k)
2880  abd(l,k) = abd(m,k)
2881  abd(m,k) = t
2882  60 CONTINUE
2883 C
2884 C COMPUTE MULTIPLIERS
2885 C
2886  t = -1.0d0/abd(m,k)
2887  CALL dscal(lm,t,abd(m+1,k),1)
2888 C
2889 C ROW ELIMINATION WITH COLUMN INDEXING
2890 C
2891  ju = min(max(ju,mu+ipvt(k)),n)
2892  mm = m
2893  IF (ju .LT. kp1) GO TO 90
2894  DO 80 j = kp1, ju
2895  l = l - 1
2896  mm = mm - 1
2897  t = abd(l,j)
2898  IF (l .EQ. mm) GO TO 70
2899  abd(l,j) = abd(mm,j)
2900  abd(mm,j) = t
2901  70 CONTINUE
2902  CALL daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
2903  80 CONTINUE
2904  90 CONTINUE
2905  GO TO 110
2906  100 CONTINUE
2907  info = k
2908  110 CONTINUE
2909  120 CONTINUE
2910  130 CONTINUE
2911  ipvt(n) = n
2912  IF (abd(m,n) .EQ. 0.0d0) info = n
2913  RETURN
max
#define max(a, b)
Definition: cascade.c:32
min
#define min(a, b)
Definition: cascade.c:31
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
Definition: dgmres.f:1276
dscal
subroutine dscal(n, da, dx, incx)
Definition: dgesv.f:1746
idamax
integer function idamax(n, dx, incx)
Definition: dgesv.f:1448

◆ dgbsl()

subroutine dgbsl ( double precision, dimension(lda,*)  ABD,
integer  LDA,
integer  N,
integer  ML,
integer  MU,
integer, dimension(*)  IPVT,
double precision, dimension(*)  B,
integer  JOB 
)
2917 C***BEGIN PROLOGUE DGBSL
2918 C***PURPOSE Solve the real band system A*X=B or TRANS(A)*X=B using
2919 C the factors computed by DGBCO or DGBFA.
2920 C***LIBRARY SLATEC (LINPACK)
2921 C***CATEGORY D2A2
2922 C***TYPE DOUBLE PRECISION (SGBSL-S, DGBSL-D, CGBSL-C)
2923 C***KEYWORDS BANDED, LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE
2924 C***AUTHOR Moler, C. B., (U. of New Mexico)
2925 C***DESCRIPTION
2926 C
2927 C DGBSL solves the double precision band system
2928 C A * X = B or TRANS(A) * X = B
2929 C using the factors computed by DGBCO or DGBFA.
2930 C
2931 C On Entry
2932 C
2933 C ABD DOUBLE PRECISION(LDA, N)
2934 C the output from DGBCO or DGBFA.
2935 C
2936 C LDA INTEGER
2937 C the leading dimension of the array ABD .
2938 C
2939 C N INTEGER
2940 C the order of the original matrix.
2941 C
2942 C ML INTEGER
2943 C number of diagonals below the main diagonal.
2944 C
2945 C MU INTEGER
2946 C number of diagonals above the main diagonal.
2947 C
2948 C IPVT INTEGER(N)
2949 C the pivot vector from DGBCO or DGBFA.
2950 C
2951 C B DOUBLE PRECISION(N)
2952 C the right hand side vector.
2953 C
2954 C JOB INTEGER
2955 C = 0 to solve A*X = B ,
2956 C = nonzero to solve TRANS(A)*X = B , where
2957 C TRANS(A) is the transpose.
2958 C
2959 C On Return
2960 C
2961 C B the solution vector X .
2962 C
2963 C Error Condition
2964 C
2965 C A division by zero will occur if the input factor contains a
2966 C zero on the diagonal. Technically this indicates singularity
2967 C but it is often caused by improper arguments or improper
2968 C setting of LDA . It will not occur if the subroutines are
2969 C called correctly and if DGBCO has set RCOND .GT. 0.0
2970 C or DGBFA has set INFO .EQ. 0 .
2971 C
2972 C To compute INVERSE(A) * C where C is a matrix
2973 C with P columns
2974 C CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z)
2975 C IF (RCOND is too small) GO TO ...
2976 C DO 10 J = 1, P
2977 C CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0)
2978 C 10 CONTINUE
2979 C
2980 C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
2981 C Stewart, LINPACK Users' Guide, SIAM, 1979.
2982 C***ROUTINES CALLED DAXPY, DDOT
2983 C***REVISION HISTORY (YYMMDD)
2984 C 780814 DATE WRITTEN
2985 C 890531 Changed all specific intrinsics to generic. (WRB)
2986 C 890831 Modified array declarations. (WRB)
2987 C 890831 REVISION DATE from Version 3.2
2988 C 891214 Prologue converted to Version 4.0 format. (BAB)
2989 C 900326 Removed duplicate information from DESCRIPTION section.
2990 C (WRB)
2991 C 920501 Reformatted the REFERENCES section. (WRB)
2992 C***END PROLOGUE DGBSL
2993  INTEGER lda,n,ml,mu,ipvt(*),job
2994  DOUBLE PRECISION abd(lda,*),b(*)
2995 C
2996  DOUBLE PRECISION ddot,t
2997  INTEGER k,kb,l,la,lb,lm,m,nm1
2998 C***FIRST EXECUTABLE STATEMENT DGBSL
2999  m = mu + ml + 1
3000  nm1 = n - 1
3001  IF (job .NE. 0) GO TO 50
3002 C
3003 C JOB = 0 , SOLVE A * X = B
3004 C FIRST SOLVE L*Y = B
3005 C
3006  IF (ml .EQ. 0) GO TO 30
3007  IF (nm1 .LT. 1) GO TO 30
3008  DO 20 k = 1, nm1
3009  lm = min(ml,n-k)
3010  l = ipvt(k)
3011  t = b(l)
3012  IF (l .EQ. k) GO TO 10
3013  b(l) = b(k)
3014  b(k) = t
3015  10 CONTINUE
3016  CALL daxpy(lm,t,abd(m+1,k),1,b(k+1),1)
3017  20 CONTINUE
3018  30 CONTINUE
3019 C
3020 C NOW SOLVE U*X = Y
3021 C
3022  DO 40 kb = 1, n
3023  k = n + 1 - kb
3024  b(k) = b(k)/abd(m,k)
3025  lm = min(k,m) - 1
3026  la = m - lm
3027  lb = k - lm
3028  t = -b(k)
3029  CALL daxpy(lm,t,abd(la,k),1,b(lb),1)
3030  40 CONTINUE
3031  GO TO 100
3032  50 CONTINUE
3033 C
3034 C JOB = NONZERO, SOLVE TRANS(A) * X = B
3035 C FIRST SOLVE TRANS(U)*Y = B
3036 C
3037  DO 60 k = 1, n
3038  lm = min(k,m) - 1
3039  la = m - lm
3040  lb = k - lm
3041  t = ddot(lm,abd(la,k),1,b(lb),1)
3042  b(k) = (b(k) - t)/abd(m,k)
3043  60 CONTINUE
3044 C
3045 C NOW SOLVE TRANS(L)*X = Y
3046 C
3047  IF (ml .EQ. 0) GO TO 90
3048  IF (nm1 .LT. 1) GO TO 90
3049  DO 80 kb = 1, nm1
3050  k = n - kb
3051  lm = min(ml,n-k)
3052  b(k) = b(k) + ddot(lm,abd(m+1,k),1,b(k+1),1)
3053  l = ipvt(k)
3054  IF (l .EQ. k) GO TO 70
3055  t = b(l)
3056  b(l) = b(k)
3057  b(k) = t
3058  70 CONTINUE
3059  80 CONTINUE
3060  90 CONTINUE
3061  100 CONTINUE
3062  RETURN
min
#define min(a, b)
Definition: cascade.c:31
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
Definition: dgmres.f:1276
ddot
double precision function ddot(N, DX, INCX, DY, INCY)
Definition: dgmres.f:2469

◆ dgefa()

subroutine dgefa ( double precision, dimension(lda,*)  A,
integer  LDA,
integer  N,
integer, dimension(*)  IPVT,
integer  INFO 
)
3066 C***BEGIN PROLOGUE DGEFA
3067 C***PURPOSE Factor a matrix using Gaussian elimination.
3068 C***LIBRARY SLATEC (LINPACK)
3069 C***CATEGORY D2A1
3070 C***TYPE DOUBLE PRECISION (SGEFA-S, DGEFA-D, CGEFA-C)
3071 C***KEYWORDS GENERAL MATRIX, LINEAR ALGEBRA, LINPACK,
3072 C MATRIX FACTORIZATION
3073 C***AUTHOR Moler, C. B., (U. of New Mexico)
3074 C***DESCRIPTION
3075 C
3076 C DGEFA factors a double precision matrix by Gaussian elimination.
3077 C
3078 C DGEFA is usually called by DGECO, but it can be called
3079 C directly with a saving in time if RCOND is not needed.
3080 C (Time for DGECO) = (1 + 9/N)*(Time for DGEFA) .
3081 C
3082 C On Entry
3083 C
3084 C A DOUBLE PRECISION(LDA, N)
3085 C the matrix to be factored.
3086 C
3087 C LDA INTEGER
3088 C the leading dimension of the array A .
3089 C
3090 C N INTEGER
3091 C the order of the matrix A .
3092 C
3093 C On Return
3094 C
3095 C A an upper triangular matrix and the multipliers
3096 C which were used to obtain it.
3097 C The factorization can be written A = L*U where
3098 C L is a product of permutation and unit lower
3099 C triangular matrices and U is upper triangular.
3100 C
3101 C IPVT INTEGER(N)
3102 C an integer vector of pivot indices.
3103 C
3104 C INFO INTEGER
3105 C = 0 normal value.
3106 C = K if U(K,K) .EQ. 0.0 . This is not an error
3107 C condition for this subroutine, but it does
3108 C indicate that DGESL or DGEDI will divide by zero
3109 C if called. Use RCOND in DGECO for a reliable
3110 C indication of singularity.
3111 C
3112 C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
3113 C Stewart, LINPACK Users' Guide, SIAM, 1979.
3114 C***ROUTINES CALLED DAXPY, DSCAL, IDAMAX
3115 C***REVISION HISTORY (YYMMDD)
3116 C 780814 DATE WRITTEN
3117 C 890831 Modified array declarations. (WRB)
3118 C 890831 REVISION DATE from Version 3.2
3119 C 891214 Prologue converted to Version 4.0 format. (BAB)
3120 C 900326 Removed duplicate information from DESCRIPTION section.
3121 C (WRB)
3122 C 920501 Reformatted the REFERENCES section. (WRB)
3123 C***END PROLOGUE DGEFA
3124  INTEGER lda,n,ipvt(*),info
3125  DOUBLE PRECISION a(lda,*)
3126 C
3127  DOUBLE PRECISION t
3128  INTEGER idamax,j,k,kp1,l,nm1
3129 C
3130 C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
3131 C
3132 C***FIRST EXECUTABLE STATEMENT DGEFA
3133  info = 0
3134  nm1 = n - 1
3135  IF (nm1 .LT. 1) GO TO 70
3136  DO 60 k = 1, nm1
3137  kp1 = k + 1
3138 C
3139 C FIND L = PIVOT INDEX
3140 C
3141  l = idamax(n-k+1,a(k,k),1) + k - 1
3142  ipvt(k) = l
3143 C
3144 C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
3145 C
3146  IF (a(l,k) .EQ. 0.0d0) GO TO 40
3147 C
3148 C INTERCHANGE IF NECESSARY
3149 C
3150  IF (l .EQ. k) GO TO 10
3151  t = a(l,k)
3152  a(l,k) = a(k,k)
3153  a(k,k) = t
3154  10 CONTINUE
3155 C
3156 C COMPUTE MULTIPLIERS
3157 C
3158  t = -1.0d0/a(k,k)
3159  CALL dscal(n-k,t,a(k+1,k),1)
3160 C
3161 C ROW ELIMINATION WITH COLUMN INDEXING
3162 C
3163  DO 30 j = kp1, n
3164  t = a(l,j)
3165  IF (l .EQ. k) GO TO 20
3166  a(l,j) = a(k,j)
3167  a(k,j) = t
3168  20 CONTINUE
3169  CALL daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
3170  30 CONTINUE
3171  GO TO 50
3172  40 CONTINUE
3173  info = k
3174  50 CONTINUE
3175  60 CONTINUE
3176  70 CONTINUE
3177  ipvt(n) = n
3178  IF (a(n,n) .EQ. 0.0d0) info = n
3179  RETURN
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
Definition: dgmres.f:1276
dscal
subroutine dscal(n, da, dx, incx)
Definition: dgesv.f:1746
idamax
integer function idamax(n, dx, incx)
Definition: dgesv.f:1448

◆ dgesl()

subroutine dgesl ( double precision, dimension(lda,*)  A,
integer  LDA,
integer  N,
integer, dimension(*)  IPVT,
double precision, dimension(*)  B,
integer  JOB 
)
3183 C***BEGIN PROLOGUE DGESL
3184 C***PURPOSE Solve the real system A*X=B or TRANS(A)*X=B using the
3185 C factors computed by DGECO or DGEFA.
3186 C***LIBRARY SLATEC (LINPACK)
3187 C***CATEGORY D2A1
3188 C***TYPE DOUBLE PRECISION (SGESL-S, DGESL-D, CGESL-C)
3189 C***KEYWORDS LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE
3190 C***AUTHOR Moler, C. B., (U. of New Mexico)
3191 C***DESCRIPTION
3192 C
3193 C DGESL solves the double precision system
3194 C A * X = B or TRANS(A) * X = B
3195 C using the factors computed by DGECO or DGEFA.
3196 C
3197 C On Entry
3198 C
3199 C A DOUBLE PRECISION(LDA, N)
3200 C the output from DGECO or DGEFA.
3201 C
3202 C LDA INTEGER
3203 C the leading dimension of the array A .
3204 C
3205 C N INTEGER
3206 C the order of the matrix A .
3207 C
3208 C IPVT INTEGER(N)
3209 C the pivot vector from DGECO or DGEFA.
3210 C
3211 C B DOUBLE PRECISION(N)
3212 C the right hand side vector.
3213 C
3214 C JOB INTEGER
3215 C = 0 to solve A*X = B ,
3216 C = nonzero to solve TRANS(A)*X = B where
3217 C TRANS(A) is the transpose.
3218 C
3219 C On Return
3220 C
3221 C B the solution vector X .
3222 C
3223 C Error Condition
3224 C
3225 C A division by zero will occur if the input factor contains a
3226 C zero on the diagonal. Technically this indicates singularity
3227 C but it is often caused by improper arguments or improper
3228 C setting of LDA . It will not occur if the subroutines are
3229 C called correctly and if DGECO has set RCOND .GT. 0.0
3230 C or DGEFA has set INFO .EQ. 0 .
3231 C
3232 C To compute INVERSE(A) * C where C is a matrix
3233 C with P columns
3234 C CALL DGECO(A,LDA,N,IPVT,RCOND,Z)
3235 C IF (RCOND is too small) GO TO ...
3236 C DO 10 J = 1, P
3237 C CALL DGESL(A,LDA,N,IPVT,C(1,J),0)
3238 C 10 CONTINUE
3239 C
3240 C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
3241 C Stewart, LINPACK Users' Guide, SIAM, 1979.
3242 C***ROUTINES CALLED DAXPY, DDOT
3243 C***REVISION HISTORY (YYMMDD)
3244 C 780814 DATE WRITTEN
3245 C 890831 Modified array declarations. (WRB)
3246 C 890831 REVISION DATE from Version 3.2
3247 C 891214 Prologue converted to Version 4.0 format. (BAB)
3248 C 900326 Removed duplicate information from DESCRIPTION section.
3249 C (WRB)
3250 C 920501 Reformatted the REFERENCES section. (WRB)
3251 C***END PROLOGUE DGESL
3252  INTEGER lda,n,ipvt(*),job
3253  DOUBLE PRECISION a(lda,*),b(*)
3254 C
3255  DOUBLE PRECISION ddot,t
3256  INTEGER k,kb,l,nm1
3257 C***FIRST EXECUTABLE STATEMENT DGESL
3258  nm1 = n - 1
3259  IF (job .NE. 0) GO TO 50
3260 C
3261 C JOB = 0 , SOLVE A * X = B
3262 C FIRST SOLVE L*Y = B
3263 C
3264  IF (nm1 .LT. 1) GO TO 30
3265  DO 20 k = 1, nm1
3266  l = ipvt(k)
3267  t = b(l)
3268  IF (l .EQ. k) GO TO 10
3269  b(l) = b(k)
3270  b(k) = t
3271  10 CONTINUE
3272  CALL daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
3273  20 CONTINUE
3274  30 CONTINUE
3275 C
3276 C NOW SOLVE U*X = Y
3277 C
3278  DO 40 kb = 1, n
3279  k = n + 1 - kb
3280  b(k) = b(k)/a(k,k)
3281  t = -b(k)
3282  CALL daxpy(k-1,t,a(1,k),1,b(1),1)
3283  40 CONTINUE
3284  GO TO 100
3285  50 CONTINUE
3286 C
3287 C JOB = NONZERO, SOLVE TRANS(A) * X = B
3288 C FIRST SOLVE TRANS(U)*Y = B
3289 C
3290  DO 60 k = 1, n
3291  t = ddot(k-1,a(1,k),1,b(1),1)
3292  b(k) = (b(k) - t)/a(k,k)
3293  60 CONTINUE
3294 C
3295 C NOW SOLVE TRANS(L)*X = Y
3296 C
3297  IF (nm1 .LT. 1) GO TO 90
3298  DO 80 kb = 1, nm1
3299  k = n - kb
3300  b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
3301  l = ipvt(k)
3302  IF (l .EQ. k) GO TO 70
3303  t = b(l)
3304  b(l) = b(k)
3305  b(k) = t
3306  70 CONTINUE
3307  80 CONTINUE
3308  90 CONTINUE
3309  100 CONTINUE
3310  RETURN
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
Definition: dgmres.f:1276
ddot
double precision function ddot(N, DX, INCX, DY, INCY)
Definition: dgmres.f:2469

◆ dintyd()

subroutine dintyd ( double precision  T,
integer  K,
double precision, dimension(nyh,*)  YH,
integer  NYH,
double precision, dimension(*)  DKY,
integer  IFLAG 
)
2288 C***BEGIN PROLOGUE DINTYD
2289 C***SUBSIDIARY
2290 C***PURPOSE Subsidiary to DDEBDF
2291 C***LIBRARY SLATEC
2292 C***TYPE DOUBLE PRECISION (INTYD-S, DINTYD-D)
2293 C***AUTHOR Watts, H. A., (SNLA)
2294 C***DESCRIPTION
2295 C
2296 C DINTYD approximates the solution and derivatives at T by polynomial
2297 C interpolation. Must be used in conjunction with the integrator
2298 C package DDEBDF.
2299 C ----------------------------------------------------------------------
2300 C DINTYD computes interpolated values of the K-th derivative of the
2301 C dependent variable vector Y, and stores it in DKY.
2302 C This routine is called by DDEBDF with K = 0,1 and T = TOUT, but may
2303 C also be called by the user for any K up to the current order.
2304 C (see detailed instructions in LSODE usage documentation.)
2305 C ----------------------------------------------------------------------
2306 C The computed values in DKY are gotten by interpolation using the
2307 C Nordsieck history array YH. This array corresponds uniquely to a
2308 C vector-valued polynomial of degree NQCUR or less, and DKY is set
2309 C to the K-th derivative of this polynomial at T.
2310 C The formula for DKY is..
2311 C Q
2312 C DKY(I) = Sum C(J,K) * (T - TN)**(J-K) * H**(-J) * YH(I,J+1)
2313 C J=K
2314 C where C(J,K) = J*(J-1)*...*(J-K+1), Q = NQCUR, TN = TCUR, H = HCUR.
2315 C The quantities NQ = NQCUR, L = NQ+1, N = NEQ, TN, and H are
2316 C communicated by common. The above sum is done in reverse order.
2317 C IFLAG is returned negative if either K or T is out of bounds.
2318 C ----------------------------------------------------------------------
2319 C
2320 C***SEE ALSO DDEBDF
2321 C***ROUTINES CALLED (NONE)
2322 C***COMMON BLOCKS DDEBD1
2323 C***REVISION HISTORY (YYMMDD)
2324 C 820301 DATE WRITTEN
2325 C 890911 Removed unnecessary intrinsics. (WRB)
2326 C 891214 Prologue converted to Version 4.0 format. (BAB)
2327 C 900328 Added TYPE section. (WRB)
2328 C 910722 Updated AUTHOR section. (ALS)
2329 C***END PROLOGUE DINTYD
2330 C
2331  INTEGER i, ic, ier, iflag, iownd, iowns, j, jb, jb2, jj, jj1,
2332  1 jp1, jstart, k, kflag, l, maxord, meth, miter, n, nfe,
2333  2 nje, nq, nqu, nst, nyh
2334  DOUBLE PRECISION c, dky, el0, h, hmin, hmxi, hu, r, rownd,
2335  1 rowns, s, t, tn, tp, uround, yh
2336  dimension yh(nyh,*),dky(*)
2337  COMMON /ddebd1/ rownd,rowns(210),el0,h,hmin,hmxi,hu,tn,uround,
2338  1 iownd(14),iowns(6),ier,jstart,kflag,l,meth,miter,
2339  2 maxord,n,nq,nst,nfe,nje,nqu
2340 C
2341 C BEGIN BLOCK PERMITTING ...EXITS TO 130
2342 C***FIRST EXECUTABLE STATEMENT DINTYD
2343  iflag = 0
2344  IF (k .LT. 0 .OR. k .GT. nq) GO TO 110
2345  tp = tn - hu*(1.0d0 + 100.0d0*uround)
2346  IF ((t - tp)*(t - tn) .LE. 0.0d0) GO TO 10
2347  iflag = -2
2348 C .........EXIT
2349  GO TO 130
2350  10 CONTINUE
2351 C
2352  s = (t - tn)/h
2353  ic = 1
2354  IF (k .EQ. 0) GO TO 30
2355  jj1 = l - k
2356  DO 20 jj = jj1, nq
2357  ic = ic*jj
2358  20 CONTINUE
2359  30 CONTINUE
2360  c = ic
2361  DO 40 i = 1, n
2362  dky(i) = c*yh(i,l)
2363  40 CONTINUE
2364  IF (k .EQ. nq) GO TO 90
2365  jb2 = nq - k
2366  DO 80 jb = 1, jb2
2367  j = nq - jb
2368  jp1 = j + 1
2369  ic = 1
2370  IF (k .EQ. 0) GO TO 60
2371  jj1 = jp1 - k
2372  DO 50 jj = jj1, j
2373  ic = ic*jj
2374  50 CONTINUE
2375  60 CONTINUE
2376  c = ic
2377  DO 70 i = 1, n
2378  dky(i) = c*yh(i,jp1) + s*dky(i)
2379  70 CONTINUE
2380  80 CONTINUE
2381 C .........EXIT
2382  IF (k .EQ. 0) GO TO 130
2383  90 CONTINUE
2384  r = h**(-k)
2385  DO 100 i = 1, n
2386  dky(i) = r*dky(i)
2387  100 CONTINUE
2388  GO TO 120
2389  110 CONTINUE
2390 C
2391  iflag = -1
2392  120 CONTINUE
2393  130 CONTINUE
2394  RETURN
2395 C ----------------------- END OF SUBROUTINE DINTYD
2396 C -----------------------

◆ dlsod()

subroutine dlsod ( external  DF,
integer  NEQ,
double precision  T,
double precision, dimension(*)  Y,
double precision  TOUT,
double precision, dimension(*)  RTOL,
double precision, dimension(*)  ATOL,
integer  IDID,
double precision, dimension(*)  YPOUT,
double precision, dimension(neq,6)  YH,
double precision, dimension(*)  YH1,
double precision, dimension(*)  EWT,
double precision, dimension(*)  SAVF,
double precision, dimension(*)  ACOR,
double precision, dimension(*)  WM,
integer, dimension(*)  IWM,
external  DJAC,
logical  INTOUT,
double precision  TSTOP,
double precision  TOLFAC,
double precision  DELSGN,
double precision, dimension(*)  RPAR,
integer, dimension(*)  IPAR 
)
941 C***BEGIN PROLOGUE DLSOD
942 C***SUBSIDIARY
943 C***PURPOSE Subsidiary to DDEBDF
944 C***LIBRARY SLATEC
945 C***TYPE DOUBLE PRECISION (LSOD-S, DLSOD-D)
946 C***AUTHOR (UNKNOWN)
947 C***DESCRIPTION
948 C
949 C DDEBDF merely allocates storage for DLSOD to relieve the user of
950 C the inconvenience of a long call list. Consequently DLSOD is used
951 C as described in the comments for DDEBDF .
952 C
953 C***SEE ALSO DDEBDF
954 C***ROUTINES CALLED D1MACH, DHSTRT, DINTYD, DSTOD, DVNRMS, XERMSG
955 C***COMMON BLOCKS DDEBD1
956 C***REVISION HISTORY (YYMMDD)
957 C 820301 DATE WRITTEN
958 C 890531 Changed all specific intrinsics to generic. (WRB)
959 C 890831 Modified array declarations. (WRB)
960 C 891214 Prologue converted to Version 4.0 format. (BAB)
961 C 900328 Added TYPE section. (WRB)
962 C 900510 Convert XERRWV calls to XERMSG calls. (RWC)
963 C***END PROLOGUE DLSOD
964 C
965  INTEGER iband, ibegin, idid, ier, iinteg, ijac, init, intflg,
966  1 iowns, ipar, iquit, itol, itstop, iwm, jstart, k, kflag,
967  2 ksteps, l, lacor, ldum, lewt, lsavf, ltol, lwm, lyh, maxnum,
968  3 maxord, meth, miter, n, natolp, neq, nfe, nje, nq, nqu,
969  4 nrtolp, nst
970  DOUBLE PRECISION absdel, acor, atol, big, d1mach, del,
971  1 delsgn, dt, dvnrms, el0, ewt,
972  2 h, ha, hmin, hmxi, hu, rowns, rpar, rtol, savf, t, tol,
973  3 told, tolfac, tout, tstop, u, wm, x, y, yh, yh1, ypout
974  LOGICAL intout
975  CHARACTER*8 xern1
976  CHARACTER*16 xern3, xern4
977 C
978  dimension y(*),ypout(*),yh(neq,6),yh1(*),ewt(*),savf(*),
979  1 acor(*),wm(*),iwm(*),rtol(*),atol(*),rpar(*),ipar(*)
980 C
981 C
982  COMMON /ddebd1/ told,rowns(210),el0,h,hmin,hmxi,hu,x,u,iquit,init,
983  1 lyh,lewt,lacor,lsavf,lwm,ksteps,ibegin,itol,
984  2 iinteg,itstop,ijac,iband,iowns(6),ier,jstart,
985  3 kflag,ldum,meth,miter,maxord,n,nq,nst,nfe,nje,nqu
986 C
987  EXTERNAL df, djac
988 C
989 C ..................................................................
990 C
991 C THE EXPENSE OF SOLVING THE PROBLEM IS MONITORED BY COUNTING THE
992 C NUMBER OF STEPS ATTEMPTED. WHEN THIS EXCEEDS MAXNUM, THE
993 C COUNTER IS RESET TO ZERO AND THE USER IS INFORMED ABOUT POSSIBLE
994 C EXCESSIVE WORK.
995  SAVE maxnum
996 C
997  DATA maxnum /500/
998 C
999 C ..................................................................
1000 C
1001 C***FIRST EXECUTABLE STATEMENT DLSOD
1002  IF (ibegin .EQ. 0) THEN
1003 C
1004 C ON THE FIRST CALL , PERFORM INITIALIZATION --
1005 C DEFINE THE MACHINE UNIT ROUNDOFF QUANTITY U BY CALLING THE
1006 C FUNCTION ROUTINE D1MACH. THE USER MUST MAKE SURE THAT THE
1007 C VALUES SET IN D1MACH ARE RELEVANT TO THE COMPUTER BEING USED.
1008 C
1009  u = d1mach(4)
1010 C -- SET ASSOCIATED MACHINE DEPENDENT PARAMETER
1011  wm(1) = sqrt(u)
1012 C -- SET TERMINATION FLAG
1013  iquit = 0
1014 C -- SET INITIALIZATION INDICATOR
1015  init = 0
1016 C -- SET COUNTER FOR ATTEMPTED STEPS
1017  ksteps = 0
1018 C -- SET INDICATOR FOR INTERMEDIATE-OUTPUT
1019  intout = .false.
1020 C -- SET START INDICATOR FOR DSTOD CODE
1021  jstart = 0
1022 C -- SET BDF METHOD INDICATOR
1023  meth = 2
1024 C -- SET MAXIMUM ORDER FOR BDF METHOD
1025  maxord = 5
1026 C -- SET ITERATION MATRIX INDICATOR
1027 C
1028  IF (ijac .EQ. 0 .AND. iband .EQ. 0) miter = 2
1029  IF (ijac .EQ. 1 .AND. iband .EQ. 0) miter = 1
1030  IF (ijac .EQ. 0 .AND. iband .EQ. 1) miter = 5
1031  IF (ijac .EQ. 1 .AND. iband .EQ. 1) miter = 4
1032 C
1033 C -- SET OTHER NECESSARY ITEMS IN COMMON BLOCK
1034  n = neq
1035  nst = 0
1036  nje = 0
1037  hmxi = 0.0d0
1038  nq = 1
1039  h = 1.0d0
1040 C -- RESET IBEGIN FOR SUBSEQUENT CALLS
1041  ibegin = 1
1042  ENDIF
1043 C
1044 C ..................................................................
1045 C
1046 C CHECK VALIDITY OF INPUT PARAMETERS ON EACH ENTRY
1047 C
1048  IF (neq .LT. 1) THEN
1049  WRITE (xern1, '(I8)') neq
1050  CALL xermsg ('SLATEC', 'DLSOD',
1051  * 'IN DDEBDF, THE NUMBER OF EQUATIONS MUST BE A ' //
1052  * 'POSITIVE INTEGER.$$YOU HAVE CALLED THE CODE WITH NEQ = ' //
1053  * xern1, 6, 1)
1054  idid=-33
1055  ENDIF
1056 C
1057  nrtolp = 0
1058  natolp = 0
1059  DO 60 k = 1, neq
1060  IF (nrtolp .LE. 0) THEN
1061  IF (rtol(k) .LT. 0.) THEN
1062  WRITE (xern1, '(I8)') k
1063  WRITE (xern3, '(1PE15.6)') rtol(k)
1064  CALL xermsg ('SLATEC', 'DLSOD',
1065  * 'IN DDEBDF, THE RELATIVE ERROR TOLERANCES MUST ' //
1066  * 'BE NON-NEGATIVE.$$YOU HAVE CALLED THE CODE WITH ' //
1067  * 'RTOL(' // xern1 // ') = ' // xern3 // '$$IN THE ' //
1068  * 'CASE OF VECTOR ERROR TOLERANCES, NO FURTHER ' //
1069  * 'CHECKING OF RTOL COMPONENTS IS DONE.', 7, 1)
1070  idid = -33
1071  IF (natolp .GT. 0) GO TO 70
1072  nrtolp = 1
1073  ELSEIF (natolp .GT. 0) THEN
1074  GO TO 50
1075  ENDIF
1076  ENDIF
1077 C
1078  IF (atol(k) .LT. 0.) THEN
1079  WRITE (xern1, '(I8)') k
1080  WRITE (xern3, '(1PE15.6)') atol(k)
1081  CALL xermsg ('SLATEC', 'DLSOD',
1082  * 'IN DDEBDF, THE ABSOLUTE ERROR ' //
1083  * 'TOLERANCES MUST BE NON-NEGATIVE.$$YOU HAVE CALLED ' //
1084  * 'THE CODE WITH ATOL(' // xern1 // ') = ' // xern3 //
1085  * '$$IN THE CASE OF VECTOR ERROR TOLERANCES, NO FURTHER '
1086  * // 'CHECKING OF ATOL COMPONENTS IS DONE.', 8, 1)
1087  idid=-33
1088  IF (nrtolp .GT. 0) GO TO 70
1089  natolp=1
1090  ENDIF
1091  50 IF (itol .EQ. 0) GO TO 70
1092  60 CONTINUE
1093 C
1094  70 IF (itstop .EQ. 1) THEN
1095  IF (sign(1.0d0,tout-t) .NE. sign(1.0d0,tstop-t) .OR.
1096  1 abs(tout-t) .GT. abs(tstop-t)) THEN
1097  WRITE (xern3, '(1PE15.6)') tout
1098  WRITE (xern4, '(1PE15.6)') tstop
1099  CALL xermsg ('SLATEC', 'DLSOD',
1100  * 'IN DDEBDF, YOU HAVE CALLED THE ' //
1101  * 'CODE WITH TOUT = ' // xern3 // '$$BUT YOU HAVE ' //
1102  * 'ALSO TOLD THE CODE NOT TO INTEGRATE PAST THE POINT ' //
1103  * 'TSTOP = ' // xern4 // ' BY SETTING INFO(4) = 1.$$' //
1104  * 'THESE INSTRUCTIONS CONFLICT.', 14, 1)
1105  idid=-33
1106  ENDIF
1107  ENDIF
1108 C
1109 C CHECK SOME CONTINUATION POSSIBILITIES
1110 C
1111  IF (init .NE. 0) THEN
1112  IF (t .EQ. tout) THEN
1113  WRITE (xern3, '(1PE15.6)') t
1114  CALL xermsg ('SLATEC', 'DLSOD',
1115  * 'IN DDEBDF, YOU HAVE CALLED THE CODE WITH T = TOUT = ' //
1116  * xern3 // '$$THIS IS NOT ALLOWED ON CONTINUATION CALLS.',
1117  * 9, 1)
1118  idid=-33
1119  ENDIF
1120 C
1121  IF (t .NE. told) THEN
1122  WRITE (xern3, '(1PE15.6)') told
1123  WRITE (xern4, '(1PE15.6)') t
1124  CALL xermsg ('SLATEC', 'DLSOD',
1125  * 'IN DDEBDF, YOU HAVE CHANGED THE VALUE OF T FROM ' //
1126  * xern3 // ' TO ' // xern4 //
1127  * ' THIS IS NOT ALLOWED ON CONTINUATION CALLS.', 10, 1)
1128  idid=-33
1129  ENDIF
1130 C
1131  IF (init .NE. 1) THEN
1132  IF (delsgn*(tout-t) .LT. 0.0d0) THEN
1133  WRITE (xern3, '(1PE15.6)') tout
1134  CALL xermsg ('SLATEC', 'DLSOD',
1135  * 'IN DDEBDF, BY CALLING THE CODE WITH TOUT = ' //
1136  * xern3 // ' YOU ARE ATTEMPTING TO CHANGE THE ' //
1137  * 'DIRECTION OF INTEGRATION.$$THIS IS NOT ALLOWED ' //
1138  * 'WITHOUT RESTARTING.', 11, 1)
1139  idid=-33
1140  ENDIF
1141  ENDIF
1142  ENDIF
1143 C
1144  IF (idid .EQ. (-33)) THEN
1145  IF (iquit .NE. (-33)) THEN
1146 C INVALID INPUT DETECTED
1147  iquit=-33
1148  ibegin=-1
1149  ELSE
1150  CALL xermsg ('SLATEC', 'DLSOD',
1151  * 'IN DDEBDF, INVALID INPUT WAS DETECTED ON ' //
1152  * 'SUCCESSIVE ENTRIES. IT IS IMPOSSIBLE TO PROCEED ' //
1153  * 'BECAUSE YOU HAVE NOT CORRECTED THE PROBLEM, ' //
1154  * 'SO EXECUTION IS BEING TERMINATED.', 12, 2)
1155  ENDIF
1156  RETURN
1157  ENDIF
1158 C
1159 C ...............................................................
1160 C
1161 C RTOL = ATOL = 0. IS ALLOWED AS VALID INPUT AND INTERPRETED
1162 C AS ASKING FOR THE MOST ACCURATE SOLUTION POSSIBLE. IN THIS
1163 C CASE, THE RELATIVE ERROR TOLERANCE RTOL IS RESET TO THE
1164 C SMALLEST VALUE 100*U WHICH IS LIKELY TO BE REASONABLE FOR
1165 C THIS METHOD AND MACHINE
1166 C
1167  DO 180 k = 1, neq
1168  IF (rtol(k) + atol(k) .GT. 0.0d0) GO TO 170
1169  rtol(k) = 100.0d0*u
1170  idid = -2
1171  170 CONTINUE
1172 C ...EXIT
1173  IF (itol .EQ. 0) GO TO 190
1174  180 CONTINUE
1175  190 CONTINUE
1176 C
1177  IF (idid .NE. (-2)) GO TO 200
1178 C RTOL=ATOL=0 ON INPUT, SO RTOL IS CHANGED TO A
1179 C SMALL POSITIVE VALUE
1180  ibegin = -1
1181  GO TO 460
1182  200 CONTINUE
1183 C BEGIN BLOCK PERMITTING ...EXITS TO 450
1184 C BEGIN BLOCK PERMITTING ...EXITS TO 430
1185 C BEGIN BLOCK PERMITTING ...EXITS TO 260
1186 C BEGIN BLOCK PERMITTING ...EXITS TO 230
1187 C
1188 C BRANCH ON STATUS OF INITIALIZATION INDICATOR
1189 C INIT=0 MEANS INITIAL DERIVATIVES AND
1190 C NOMINAL STEP SIZE
1191 C AND DIRECTION NOT YET SET
1192 C INIT=1 MEANS NOMINAL STEP SIZE AND
1193 C DIRECTION NOT YET SET INIT=2 MEANS NO
1194 C FURTHER INITIALIZATION REQUIRED
1195 C
1196  IF (init .EQ. 0) GO TO 210
1197 C ......EXIT
1198  IF (init .EQ. 1) GO TO 230
1199 C .........EXIT
1200  GO TO 260
1201  210 CONTINUE
1202 C
1203 C ................................................
1204 C
1205 C MORE INITIALIZATION --
1206 C -- EVALUATE INITIAL
1207 C DERIVATIVES
1208 C
1209  init = 1
1210  CALL df(t,y,yh(1,2),rpar,ipar)
1211  nfe = 1
1212 C ...EXIT
1213  IF (t .NE. tout) GO TO 230
1214  idid = 2
1215  DO 220 l = 1, neq
1216  ypout(l) = yh(l,2)
1217  220 CONTINUE
1218  told = t
1219 C ............EXIT
1220  GO TO 450
1221  230 CONTINUE
1222 C
1223 C -- COMPUTE INITIAL STEP SIZE
1224 C -- SAVE SIGN OF INTEGRATION DIRECTION
1225 C -- SET INDEPENDENT AND DEPENDENT VARIABLES
1226 C X AND YH(*) FOR DSTOD
1227 C
1228  ltol = 1
1229  DO 240 l = 1, neq
1230  IF (itol .EQ. 1) ltol = l
1231  tol = rtol(ltol)*abs(y(l)) + atol(ltol)
1232  IF (tol .EQ. 0.0d0) GO TO 390
1233  ewt(l) = tol
1234  240 CONTINUE
1235 C
1236  big = sqrt(d1mach(2))
1237  CALL dhstrt(df,neq,t,tout,y,yh(1,2),ewt,1,u,big,
1238  1 yh(1,3),yh(1,4),yh(1,5),yh(1,6),rpar,
1239  2 ipar,h)
1240 C
1241  delsgn = sign(1.0d0,tout-t)
1242  x = t
1243  DO 250 l = 1, neq
1244  yh(l,1) = y(l)
1245  yh(l,2) = h*yh(l,2)
1246  250 CONTINUE
1247  init = 2
1248  260 CONTINUE
1249 C
1250 C ......................................................
1251 C
1252 C ON EACH CALL SET INFORMATION WHICH DETERMINES THE
1253 C ALLOWED INTERVAL OF INTEGRATION BEFORE RETURNING
1254 C WITH AN ANSWER AT TOUT
1255 C
1256  del = tout - t
1257  absdel = abs(del)
1258 C
1259 C ......................................................
1260 C
1261 C IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND
1262 C RETURN
1263 C
1264  270 CONTINUE
1265 C BEGIN BLOCK PERMITTING ...EXITS TO 400
1266 C BEGIN BLOCK PERMITTING ...EXITS TO 380
1267  IF (abs(x-t) .LT. absdel) GO TO 290
1268  CALL dintyd(tout,0,yh,neq,y,intflg)
1269  CALL dintyd(tout,1,yh,neq,ypout,intflg)
1270  idid = 3
1271  IF (x .NE. tout) GO TO 280
1272  idid = 2
1273  intout = .false.
1274  280 CONTINUE
1275  t = tout
1276  told = t
1277 C ..................EXIT
1278  GO TO 450
1279  290 CONTINUE
1280 C
1281 C IF CANNOT GO PAST TSTOP AND SUFFICIENTLY
1282 C CLOSE, EXTRAPOLATE AND RETURN
1283 C
1284  IF (itstop .NE. 1) GO TO 310
1285  IF (abs(tstop-x) .GE. 100.0d0*u*abs(x))
1286  1 GO TO 310
1287  dt = tout - x
1288  DO 300 l = 1, neq
1289  y(l) = yh(l,1) + (dt/h)*yh(l,2)
1290  300 CONTINUE
1291  CALL df(tout,y,ypout,rpar,ipar)
1292  nfe = nfe + 1
1293  idid = 3
1294  t = tout
1295  told = t
1296 C ..................EXIT
1297  GO TO 450
1298  310 CONTINUE
1299 C
1300  IF (iinteg .EQ. 0 .OR. .NOT.intout) GO TO 320
1301 C
1302 C INTERMEDIATE-OUTPUT MODE
1303 C
1304  idid = 1
1305  GO TO 370
1306  320 CONTINUE
1307 C
1308 C .............................................
1309 C
1310 C MONITOR NUMBER OF STEPS ATTEMPTED
1311 C
1312  IF (ksteps .LE. maxnum) GO TO 330
1313 C
1314 C A SIGNIFICANT AMOUNT OF WORK HAS BEEN
1315 C EXPENDED
1316  idid = -1
1317  ksteps = 0
1318  ibegin = -1
1319  GO TO 370
1320  330 CONTINUE
1321 C
1322 C ..........................................
1323 C
1324 C LIMIT STEP SIZE AND SET WEIGHT VECTOR
1325 C
1326  hmin = 100.0d0*u*abs(x)
1327  ha = max(abs(h),hmin)
1328  IF (itstop .EQ. 1)
1329  1 ha = min(ha,abs(tstop-x))
1330  h = sign(ha,h)
1331  ltol = 1
1332  DO 340 l = 1, neq
1333  IF (itol .EQ. 1) ltol = l
1334  ewt(l) = rtol(ltol)*abs(yh(l,1))
1335  1 + atol(ltol)
1336 C .........EXIT
1337  IF (ewt(l) .LE. 0.0d0) GO TO 380
1338  340 CONTINUE
1339  tolfac = u*dvnrms(neq,yh,ewt)
1340 C .........EXIT
1341  IF (tolfac .LE. 1.0d0) GO TO 400
1342 C
1343 C TOLERANCES TOO SMALL
1344  idid = -2
1345  tolfac = 2.0d0*tolfac
1346  rtol(1) = tolfac*rtol(1)
1347  atol(1) = tolfac*atol(1)
1348  IF (itol .EQ. 0) GO TO 360
1349  DO 350 l = 2, neq
1350  rtol(l) = tolfac*rtol(l)
1351  atol(l) = tolfac*atol(l)
1352  350 CONTINUE
1353  360 CONTINUE
1354  ibegin = -1
1355  370 CONTINUE
1356 C ............EXIT
1357  GO TO 430
1358  380 CONTINUE
1359 C
1360 C RELATIVE ERROR CRITERION INAPPROPRIATE
1361  390 CONTINUE
1362  idid = -3
1363  ibegin = -1
1364 C .........EXIT
1365  GO TO 430
1366  400 CONTINUE
1367 C
1368 C ...................................................
1369 C
1370 C TAKE A STEP
1371 C
1372  CALL dstod(neq,y,yh,neq,yh1,ewt,savf,acor,wm,iwm,
1373  1 df,djac,rpar,ipar)
1374 C
1375  jstart = -2
1376  intout = .true.
1377  IF (kflag .EQ. 0) GO TO 270
1378 C
1379 C ......................................................
1380 C
1381  IF (kflag .EQ. -1) GO TO 410
1382 C
1383 C REPEATED CORRECTOR CONVERGENCE FAILURES
1384  idid = -6
1385  ibegin = -1
1386  GO TO 420
1387  410 CONTINUE
1388 C
1389 C REPEATED ERROR TEST FAILURES
1390  idid = -7
1391  ibegin = -1
1392  420 CONTINUE
1393  430 CONTINUE
1394 C
1395 C .........................................................
1396 C
1397 C STORE VALUES BEFORE RETURNING TO
1398 C DDEBDF
1399  DO 440 l = 1, neq
1400  y(l) = yh(l,1)
1401  ypout(l) = yh(l,2)/h
1402  440 CONTINUE
1403  t = x
1404  told = t
1405  intout = .false.
1406  450 CONTINUE
1407  460 CONTINUE
1408  RETURN
init
subroutine init(nktet, inodfa, ipofa, netet_)
Definition: init.f:20
dintyd
subroutine dintyd(T, K, YH, NYH, DKY, IFLAG)
Definition: ddebdf.f:2288
max
#define max(a, b)
Definition: cascade.c:32
d1mach
double precision function d1mach(I)
Definition: ddeabm.f:2012
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
xermsg
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: ddeabm.f:1265
min
#define min(a, b)
Definition: cascade.c:31
dvnrms
double precision function dvnrms(N, V, W)
Definition: ddebdf.f:2252
djac
subroutine djac(x, u, pd, nrowpd, rpar, nev)
Definition: subspace.f:163
dstod
subroutine dstod(NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, DF, DJAC, RPAR, IPAR)
Definition: ddebdf.f:1413
dhstrt
subroutine dhstrt(DF, NEQ, A, B, Y, YPRIME, ETOL, MORDER, SMALL, BIG, SPY, PV, YP, SF, RPAR, IPAR, H)
Definition: ddeabm.f:4092

◆ dpjac()

subroutine dpjac ( integer  NEQ,
double precision, dimension(*)  Y,
double precision, dimension(nyh,*)  YH,
integer  NYH,
double precision, dimension(*)  EWT,
double precision, dimension(*)  FTEM,
double precision, dimension(*)  SAVF,
double precision, dimension(*)  WM,
integer, dimension(*)  IWM,
external  DF,
external  DJAC,
double precision, dimension(*)  RPAR,
integer, dimension(*)  IPAR 
)
2401 C***BEGIN PROLOGUE DPJAC
2402 C***SUBSIDIARY
2403 C***PURPOSE Subsidiary to DDEBDF
2404 C***LIBRARY SLATEC
2405 C***TYPE DOUBLE PRECISION (PJAC-S, DPJAC-D)
2406 C***AUTHOR Watts, H. A., (SNLA)
2407 C***DESCRIPTION
2408 C
2409 C DPJAC sets up the iteration matrix (involving the Jacobian) for the
2410 C integration package DDEBDF.
2411 C
2412 C***SEE ALSO DDEBDF
2413 C***ROUTINES CALLED DGBFA, DGEFA, DVNRMS
2414 C***COMMON BLOCKS DDEBD1
2415 C***REVISION HISTORY (YYMMDD)
2416 C 820301 DATE WRITTEN
2417 C 890531 Changed all specific intrinsics to generic. (WRB)
2418 C 890911 Removed unnecessary intrinsics. (WRB)
2419 C 891214 Prologue converted to Version 4.0 format. (BAB)
2420 C 900328 Added TYPE section. (WRB)
2421 C 910722 Updated AUTHOR section. (ALS)
2422 C 920422 Changed DIMENSION statement. (WRB)
2423 C***END PROLOGUE DPJAC
2424 C
2425  INTEGER i, i1, i2, ier, ii, iownd, iowns, ipar, iwm, j, j1,
2426  1 jj, jstart, kflag, l, lenp, maxord, mba, mband,
2427  2 meb1, meband, meth, miter, ml, ml3, mu, n, neq,
2428  3 nfe, nje, nq, nqu, nst, nyh
2429  DOUBLE PRECISION con, di, dvnrms, el0, ewt,
2430  1 fac, ftem, h, hl0, hmin, hmxi, hu, r, r0, rownd, rowns,
2431  2 rpar, savf, srur, tn, uround, wm, y, yh, yi, yj, yjj
2432  EXTERNAL df, djac
2433  dimension y(*),yh(nyh,*),ewt(*),ftem(*),savf(*),wm(*),iwm(*),
2434  1 rpar(*),ipar(*)
2435  COMMON /ddebd1/ rownd,rowns(210),el0,h,hmin,hmxi,hu,tn,uround,
2436  1 iownd(14),iowns(6),ier,jstart,kflag,l,meth,miter,
2437  2 maxord,n,nq,nst,nfe,nje,nqu
2438 C ------------------------------------------------------------------
2439 C DPJAC IS CALLED BY DSTOD TO COMPUTE AND PROCESS THE MATRIX
2440 C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
2441 C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE DJAC IF
2442 C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
2443 C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
2444 C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN
2445 C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION
2446 C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE
2447 C BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5.
2448 C
2449 C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
2450 C WITH DPJAC USES THE FOLLOWING..
2451 C Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
2452 C FTEM = WORK ARRAY OF LENGTH N (ACOR IN DSTOD ).
2453 C SAVF = ARRAY CONTAINING DF EVALUATED AT PREDICTED Y.
2454 C WM = DOUBLE PRECISION WORK SPACE FOR MATRICES. ON OUTPUT IT
2455 C CONTAINS THE
2456 C INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU
2457 C DECOMPOSITION OF P IF MITER IS 1, 2 , 4, OR 5.
2458 C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
2459 C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
2460 C WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN
2461 C INCREMENTS. WM(2) = H*EL0, SAVED FOR LATER USE IF MITER =
2462 C 3.
2463 C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING
2464 C AT IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS
2465 C THE BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER
2466 C IS 4 OR 5.
2467 C EL0 = EL(1) (INPUT).
2468 C IER = OUTPUT ERROR FLAG, = 0 IF NO TROUBLE, .NE. 0 IF
2469 C P MATRIX FOUND TO BE SINGULAR.
2470 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
2471 C MITER, N, NFE, AND NJE.
2472 C-----------------------------------------------------------------------
2473 C BEGIN BLOCK PERMITTING ...EXITS TO 240
2474 C BEGIN BLOCK PERMITTING ...EXITS TO 220
2475 C BEGIN BLOCK PERMITTING ...EXITS TO 130
2476 C BEGIN BLOCK PERMITTING ...EXITS TO 70
2477 C***FIRST EXECUTABLE STATEMENT DPJAC
2478  nje = nje + 1
2479  hl0 = h*el0
2480  GO TO (10,40,90,140,170), miter
2481 C IF MITER = 1, CALL DJAC AND MULTIPLY BY SCALAR.
2482 C -----------------------
2483  10 CONTINUE
2484  lenp = n*n
2485  DO 20 i = 1, lenp
2486  wm(i+2) = 0.0d0
2487  20 CONTINUE
2488  CALL djac(tn,y,wm(3),n,rpar,ipar)
2489  con = -hl0
2490  DO 30 i = 1, lenp
2491  wm(i+2) = wm(i+2)*con
2492  30 CONTINUE
2493 C ...EXIT
2494  GO TO 70
2495 C IF MITER = 2, MAKE N CALLS TO DF TO APPROXIMATE J.
2496 C --------------------
2497  40 CONTINUE
2498  fac = dvnrms(n,savf,ewt)
2499  r0 = 1000.0d0*abs(h)*uround*n*fac
2500  IF (r0 .EQ. 0.0d0) r0 = 1.0d0
2501  srur = wm(1)
2502  j1 = 2
2503  DO 60 j = 1, n
2504  yj = y(j)
2505  r = max(srur*abs(yj),r0*ewt(j))
2506  y(j) = y(j) + r
2507  fac = -hl0/r
2508  CALL df(tn,y,ftem,rpar,ipar)
2509  DO 50 i = 1, n
2510  wm(i+j1) = (ftem(i) - savf(i))*fac
2511  50 CONTINUE
2512  y(j) = yj
2513  j1 = j1 + n
2514  60 CONTINUE
2515  nfe = nfe + n
2516  70 CONTINUE
2517 C ADD IDENTITY MATRIX.
2518 C -------------------------------------------------
2519  j = 3
2520  DO 80 i = 1, n
2521  wm(j) = wm(j) + 1.0d0
2522  j = j + (n + 1)
2523  80 CONTINUE
2524 C DO LU DECOMPOSITION ON P.
2525 C --------------------------------------------
2526  CALL dgefa(wm(3),n,n,iwm(21),ier)
2527 C .........EXIT
2528  GO TO 240
2529 C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND
2530 C P. ---------
2531  90 CONTINUE
2532  wm(2) = hl0
2533  ier = 0
2534  r = el0*0.1d0
2535  DO 100 i = 1, n
2536  y(i) = y(i) + r*(h*savf(i) - yh(i,2))
2537  100 CONTINUE
2538  CALL df(tn,y,wm(3),rpar,ipar)
2539  nfe = nfe + 1
2540  DO 120 i = 1, n
2541  r0 = h*savf(i) - yh(i,2)
2542  di = 0.1d0*r0 - h*(wm(i+2) - savf(i))
2543  wm(i+2) = 1.0d0
2544  IF (abs(r0) .LT. uround*ewt(i)) GO TO 110
2545 C .........EXIT
2546  IF (abs(di) .EQ. 0.0d0) GO TO 130
2547  wm(i+2) = 0.1d0*r0/di
2548  110 CONTINUE
2549  120 CONTINUE
2550 C .........EXIT
2551  GO TO 240
2552  130 CONTINUE
2553  ier = -1
2554 C ......EXIT
2555  GO TO 240
2556 C IF MITER = 4, CALL DJAC AND MULTIPLY BY SCALAR.
2557 C -----------------------
2558  140 CONTINUE
2559  ml = iwm(1)
2560  mu = iwm(2)
2561  ml3 = 3
2562  mband = ml + mu + 1
2563  meband = mband + ml
2564  lenp = meband*n
2565  DO 150 i = 1, lenp
2566  wm(i+2) = 0.0d0
2567  150 CONTINUE
2568  CALL djac(tn,y,wm(ml3),meband,rpar,ipar)
2569  con = -hl0
2570  DO 160 i = 1, lenp
2571  wm(i+2) = wm(i+2)*con
2572  160 CONTINUE
2573 C ...EXIT
2574  GO TO 220
2575 C IF MITER = 5, MAKE MBAND CALLS TO DF TO APPROXIMATE J.
2576 C ----------------
2577  170 CONTINUE
2578  ml = iwm(1)
2579  mu = iwm(2)
2580  mband = ml + mu + 1
2581  mba = min(mband,n)
2582  meband = mband + ml
2583  meb1 = meband - 1
2584  srur = wm(1)
2585  fac = dvnrms(n,savf,ewt)
2586  r0 = 1000.0d0*abs(h)*uround*n*fac
2587  IF (r0 .EQ. 0.0d0) r0 = 1.0d0
2588  DO 210 j = 1, mba
2589  DO 180 i = j, n, mband
2590  yi = y(i)
2591  r = max(srur*abs(yi),r0*ewt(i))
2592  y(i) = y(i) + r
2593  180 CONTINUE
2594  CALL df(tn,y,ftem,rpar,ipar)
2595  DO 200 jj = j, n, mband
2596  y(jj) = yh(jj,1)
2597  yjj = y(jj)
2598  r = max(srur*abs(yjj),r0*ewt(jj))
2599  fac = -hl0/r
2600  i1 = max(jj-mu,1)
2601  i2 = min(jj+ml,n)
2602  ii = jj*meb1 - ml + 2
2603  DO 190 i = i1, i2
2604  wm(ii+i) = (ftem(i) - savf(i))*fac
2605  190 CONTINUE
2606  200 CONTINUE
2607  210 CONTINUE
2608  nfe = nfe + mba
2609  220 CONTINUE
2610 C ADD IDENTITY MATRIX.
2611 C -------------------------------------------------
2612  ii = mband + 2
2613  DO 230 i = 1, n
2614  wm(ii) = wm(ii) + 1.0d0
2615  ii = ii + meband
2616  230 CONTINUE
2617 C DO LU DECOMPOSITION OF P.
2618 C --------------------------------------------
2619  CALL dgbfa(wm(3),meband,n,ml,mu,iwm(21),ier)
2620  240 CONTINUE
2621  RETURN
2622 C ----------------------- END OF SUBROUTINE DPJAC
2623 C -----------------------
max
#define max(a, b)
Definition: cascade.c:32
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
min
#define min(a, b)
Definition: cascade.c:31
dgbfa
subroutine dgbfa(ABD, LDA, N, ML, MU, IPVT, INFO)
Definition: ddebdf.f:2730
dvnrms
double precision function dvnrms(N, V, W)
Definition: ddebdf.f:2252
djac
subroutine djac(x, u, pd, nrowpd, rpar, nev)
Definition: subspace.f:163
dgefa
subroutine dgefa(A, LDA, N, IPVT, INFO)
Definition: ddebdf.f:3066

◆ dslvs()

subroutine dslvs ( double precision, dimension(*)  WM,
integer, dimension(*)  IWM,
double precision, dimension(*)  X,
double precision, dimension(*)  TEM 
)
2627 C***BEGIN PROLOGUE DSLVS
2628 C***SUBSIDIARY
2629 C***PURPOSE Subsidiary to DDEBDF
2630 C***LIBRARY SLATEC
2631 C***TYPE DOUBLE PRECISION (SLVS-S, DSLVS-D)
2632 C***AUTHOR Watts, H. A., (SNLA)
2633 C***DESCRIPTION
2634 C
2635 C DSLVS solves the linear system in the iteration scheme for the
2636 C integrator package DDEBDF.
2637 C
2638 C***SEE ALSO DDEBDF
2639 C***ROUTINES CALLED DGBSL, DGESL
2640 C***COMMON BLOCKS DDEBD1
2641 C***REVISION HISTORY (YYMMDD)
2642 C 820301 DATE WRITTEN
2643 C 890531 Changed all specific intrinsics to generic. (WRB)
2644 C 891214 Prologue converted to Version 4.0 format. (BAB)
2645 C 900328 Added TYPE section. (WRB)
2646 C 910722 Updated AUTHOR section. (ALS)
2647 C 920422 Changed DIMENSION statement. (WRB)
2648 C***END PROLOGUE DSLVS
2649 C
2650  INTEGER i, ier, iownd, iowns, iwm, jstart, kflag, l, maxord,
2651  1 meband, meth, miter, ml, mu, n, nfe, nje, nq, nqu, nst
2652  DOUBLE PRECISION di, el0, h, hl0, hmin, hmxi, hu, phl0,
2653  1 r, rownd, rowns, tem, tn, uround, wm, x
2654  dimension wm(*), iwm(*), x(*), tem(*)
2655  COMMON /ddebd1/ rownd,rowns(210),el0,h,hmin,hmxi,hu,tn,uround,
2656  1 iownd(14),iowns(6),ier,jstart,kflag,l,meth,miter,
2657  2 maxord,n,nq,nst,nfe,nje,nqu
2658 C ------------------------------------------------------------------
2659 C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING
2660 C FROM A CHORD ITERATION. IT IS CALLED BY DSTOD IF MITER .NE. 0.
2661 C IF MITER IS 1 OR 2, IT CALLS DGESL TO ACCOMPLISH THIS.
2662 C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL
2663 C MATRIX, AND THEN COMPUTES THE SOLUTION.
2664 C IF MITER IS 4 OR 5, IT CALLS DGBSL.
2665 C COMMUNICATION WITH DSLVS USES THE FOLLOWING VARIABLES..
2666 C WM = DOUBLE PRECISION WORK SPACE CONTAINING THE INVERSE DIAGONAL
2667 C MATRIX IF MITER
2668 C IS 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE.
2669 C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
2670 C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
2671 C WM(1) = SQRT(UROUND) (NOT USED HERE),
2672 C WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER =
2673 C 3.
2674 C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING
2675 C AT IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS
2676 C THE BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS
2677 C 4 OR 5.
2678 C X = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION
2679 C VECTOR ON OUTPUT, OF LENGTH N.
2680 C TEM = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION.
2681 C IER = OUTPUT FLAG (IN COMMON). IER = 0 IF NO TROUBLE OCCURRED.
2682 C IER = -1 IF A SINGULAR MATRIX AROSE WITH MITER = 3.
2683 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N.
2684 C-----------------------------------------------------------------------
2685 C BEGIN BLOCK PERMITTING ...EXITS TO 80
2686 C BEGIN BLOCK PERMITTING ...EXITS TO 60
2687 C***FIRST EXECUTABLE STATEMENT DSLVS
2688  ier = 0
2689  GO TO (10,10,20,70,70), miter
2690  10 CONTINUE
2691  CALL dgesl(wm(3),n,n,iwm(21),x,0)
2692 C ......EXIT
2693  GO TO 80
2694 C
2695  20 CONTINUE
2696  phl0 = wm(2)
2697  hl0 = h*el0
2698  wm(2) = hl0
2699  IF (hl0 .EQ. phl0) GO TO 40
2700  r = hl0/phl0
2701  DO 30 i = 1, n
2702  di = 1.0d0 - r*(1.0d0 - 1.0d0/wm(i+2))
2703 C .........EXIT
2704  IF (abs(di) .EQ. 0.0d0) GO TO 60
2705  wm(i+2) = 1.0d0/di
2706  30 CONTINUE
2707  40 CONTINUE
2708  DO 50 i = 1, n
2709  x(i) = wm(i+2)*x(i)
2710  50 CONTINUE
2711 C ......EXIT
2712  GO TO 80
2713  60 CONTINUE
2714  ier = -1
2715 C ...EXIT
2716  GO TO 80
2717 C
2718  70 CONTINUE
2719  ml = iwm(1)
2720  mu = iwm(2)
2721  meband = 2*ml + mu + 1
2722  CALL dgbsl(wm(3),meband,n,ml,mu,iwm(21),x,0)
2723  80 CONTINUE
2724  RETURN
2725 C ----------------------- END OF SUBROUTINE DSLVS
2726 C -----------------------
dgesl
subroutine dgesl(A, LDA, N, IPVT, B, JOB)
Definition: ddebdf.f:3183
dgbsl
subroutine dgbsl(ABD, LDA, N, ML, MU, IPVT, B, JOB)
Definition: ddebdf.f:2917

◆ dstod()

subroutine dstod ( integer  NEQ,
double precision, dimension(*)  Y,
double precision, dimension(nyh,*)  YH,
integer  NYH,
double precision, dimension(*)  YH1,
double precision, dimension(*)  EWT,
double precision, dimension(*)  SAVF,
double precision, dimension(*)  ACOR,
double precision, dimension(*)  WM,
integer, dimension(*)  IWM,
external  DF,
external  DJAC,
double precision, dimension(*)  RPAR,
integer, dimension(*)  IPAR 
)
1413 C***BEGIN PROLOGUE DSTOD
1414 C***SUBSIDIARY
1415 C***PURPOSE Subsidiary to DDEBDF
1416 C***LIBRARY SLATEC
1417 C***TYPE DOUBLE PRECISION (STOD-S, DSTOD-D)
1418 C***AUTHOR Watts, H. A., (SNLA)
1419 C***DESCRIPTION
1420 C
1421 C DSTOD integrates a system of first order odes over one step in the
1422 C integrator package DDEBDF.
1423 C ----------------------------------------------------------------------
1424 C DSTOD performs one step of the integration of an initial value
1425 C problem for a system of ordinary differential equations.
1426 C Note.. DSTOD is independent of the value of the iteration method
1427 C indicator MITER, when this is .NE. 0, and hence is independent
1428 C of the type of chord method used, or the Jacobian structure.
1429 C Communication with DSTOD is done with the following variables..
1430 C
1431 C Y = An array of length .GE. N used as the Y argument in
1432 C all calls to DF and DJAC.
1433 C NEQ = Integer array containing problem size in NEQ(1), and
1434 C passed as the NEQ argument in all calls to DF and DJAC.
1435 C YH = An NYH by LMAX array containing the dependent variables
1436 C and their approximate scaled derivatives, where
1437 C LMAX = MAXORD + 1. YH(I,J+1) contains the approximate
1438 C J-th derivative of Y(I), scaled by H**J/FACTORIAL(J)
1439 C (J = 0,1,...,NQ). On entry for the first step, the first
1440 C two columns of YH must be set from the initial values.
1441 C NYH = A constant integer .GE. N, the first dimension of YH.
1442 C YH1 = A one-dimensional array occupying the same space as YH.
1443 C EWT = An array of N elements with which the estimated local
1444 C errors in YH are compared.
1445 C SAVF = An array of working storage, of length N.
1446 C ACOR = A work array of length N, used for the accumulated
1447 C corrections. On a successful return, ACOR(I) contains
1448 C the estimated one-step local error in Y(I).
1449 C WM,IWM = DOUBLE PRECISION and INTEGER work arrays associated with
1450 C matrix operations in chord iteration (MITER .NE. 0).
1451 C DPJAC = Name of routine to evaluate and preprocess Jacobian matrix
1452 C if a chord method is being used.
1453 C DSLVS = Name of routine to solve linear system in chord iteration.
1454 C H = The step size to be attempted on the next step.
1455 C H is altered by the error control algorithm during the
1456 C problem. H can be either positive or negative, but its
1457 C sign must remain constant throughout the problem.
1458 C HMIN = The minimum absolute value of the step size H to be used.
1459 C HMXI = Inverse of the maximum absolute value of H to be used.
1460 C HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
1461 C HMIN and HMXI may be changed at any time, but will not
1462 C take effect until the next change of H is considered.
1463 C TN = The independent variable. TN is updated on each step taken.
1464 C JSTART = An integer used for input only, with the following
1465 C values and meanings..
1466 C 0 Perform the first step.
1467 C .GT.0 Take a new step continuing from the last.
1468 C -1 Take the next step with a new value of H, MAXORD,
1469 C N, METH, MITER, and/or matrix parameters.
1470 C -2 Take the next step with a new value of H,
1471 C but with other inputs unchanged.
1472 C On return, JSTART is set to 1 to facilitate continuation.
1473 C KFLAG = a completion code with the following meanings..
1474 C 0 The step was successful.
1475 C -1 The requested error could not be achieved.
1476 C -2 Corrector convergence could not be achieved.
1477 C A return with KFLAG = -1 or -2 means either
1478 C ABS(H) = HMIN or 10 consecutive failures occurred.
1479 C On a return with KFLAG negative, the values of TN and
1480 C the YH array are as of the beginning of the last
1481 C step, and H is the last step size attempted.
1482 C MAXORD = The maximum order of integration method to be allowed.
1483 C METH/MITER = The method flags. See description in driver.
1484 C N = The number of first-order differential equations.
1485 C ----------------------------------------------------------------------
1486 C
1487 C***SEE ALSO DDEBDF
1488 C***ROUTINES CALLED DCFOD, DPJAC, DSLVS, DVNRMS
1489 C***COMMON BLOCKS DDEBD1
1490 C***REVISION HISTORY (YYMMDD)
1491 C 820301 DATE WRITTEN
1492 C 890531 Changed all specific intrinsics to generic. (WRB)
1493 C 890911 Removed unnecessary intrinsics. (WRB)
1494 C 891214 Prologue converted to Version 4.0 format. (BAB)
1495 C 900328 Added TYPE section. (WRB)
1496 C 910722 Updated AUTHOR section. (ALS)
1497 C 920422 Changed DIMENSION statement. (WRB)
1498 C***END PROLOGUE DSTOD
1499 C
1500  INTEGER i, i1, ialth, ier, iod, iownd, ipar, ipup, iredo, iret,
1501  1 iwm, j, jb, jstart, kflag, ksteps, l, lmax, m, maxord,
1502  2 meo, meth, miter, n, ncf, neq, newq, nfe, nje, nq, nqnyh,
1503  3 nqu, nst, nstepj, nyh
1504  DOUBLE PRECISION acor, conit, crate, dcon, ddn,
1505  1 del, delp, dsm, dup, dvnrms, el, el0, elco,
1506  2 ewt, exdn, exsm, exup, h, hmin, hmxi, hold, hu, r, rc,
1507  3 rh, rhdn, rhsm, rhup, rmax, rownd, rpar, savf, tesco,
1508  4 tn, told, uround, wm, y, yh, yh1
1509  EXTERNAL df, djac
1510 C
1511  dimension y(*),yh(nyh,*),yh1(*),ewt(*),savf(*),acor(*),wm(*),
1512  1 iwm(*),rpar(*),ipar(*)
1513  COMMON /ddebd1/ rownd,conit,crate,el(13),elco(13,12),hold,rc,rmax,
1514  1 tesco(3,12),el0,h,hmin,hmxi,hu,tn,uround,iownd(7),
1515  2 ksteps,iod(6),ialth,ipup,lmax,meo,nqnyh,nstepj,
1516  3 ier,jstart,kflag,l,meth,miter,maxord,n,nq,nst,nfe,
1517  4 nje,nqu
1518 C
1519 C
1520 C BEGIN BLOCK PERMITTING ...EXITS TO 690
1521 C BEGIN BLOCK PERMITTING ...EXITS TO 60
1522 C***FIRST EXECUTABLE STATEMENT DSTOD
1523  kflag = 0
1524  told = tn
1525  ncf = 0
1526  IF (jstart .GT. 0) GO TO 160
1527  IF (jstart .EQ. -1) GO TO 10
1528  IF (jstart .EQ. -2) GO TO 90
1529 C ---------------------------------------------------------
1530 C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER
1531 C VARIABLES ARE INITIALIZED. RMAX IS THE MAXIMUM RATIO BY
1532 C WHICH H CAN BE INCREASED IN A SINGLE STEP. IT IS
1533 C INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL INITIAL H,
1534 C BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE OCCURS
1535 C (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT
1536 C 2 FOR THE NEXT INCREASE.
1537 C ---------------------------------------------------------
1538  lmax = maxord + 1
1539  nq = 1
1540  l = 2
1541  ialth = 2
1542  rmax = 10000.0d0
1543  rc = 0.0d0
1544  el0 = 1.0d0
1545  crate = 0.7d0
1546  delp = 0.0d0
1547  hold = h
1548  meo = meth
1549  nstepj = 0
1550  iret = 3
1551  GO TO 50
1552  10 CONTINUE
1553 C BEGIN BLOCK PERMITTING ...EXITS TO 30
1554 C ------------------------------------------------------
1555 C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN
1556 C JSTART = -1. IPUP IS SET TO MITER TO FORCE A MATRIX
1557 C UPDATE. IF AN ORDER INCREASE IS ABOUT TO BE
1558 C CONSIDERED (IALTH = 1), IALTH IS RESET TO 2 TO
1559 C POSTPONE CONSIDERATION ONE MORE STEP. IF THE CALLER
1560 C HAS CHANGED METH, DCFOD IS CALLED TO RESET THE
1561 C COEFFICIENTS OF THE METHOD. IF THE CALLER HAS
1562 C CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT
1563 C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN
1564 C ACCORDINGLY. IF H IS TO BE CHANGED, YH MUST BE
1565 C RESCALED. IF H OR METH IS BEING CHANGED, IALTH IS
1566 C RESET TO L = NQ + 1 TO PREVENT FURTHER CHANGES IN H
1567 C FOR THAT MANY STEPS.
1568 C ------------------------------------------------------
1569  ipup = miter
1570  lmax = maxord + 1
1571  IF (ialth .EQ. 1) ialth = 2
1572  IF (meth .EQ. meo) GO TO 20
1573  CALL dcfod(meth,elco,tesco)
1574  meo = meth
1575 C ......EXIT
1576  IF (nq .GT. maxord) GO TO 30
1577  ialth = l
1578  iret = 1
1579 C ............EXIT
1580  GO TO 60
1581  20 CONTINUE
1582  IF (nq .LE. maxord) GO TO 90
1583  30 CONTINUE
1584  nq = maxord
1585  l = lmax
1586  DO 40 i = 1, l
1587  el(i) = elco(i,nq)
1588  40 CONTINUE
1589  nqnyh = nq*nyh
1590  rc = rc*el(1)/el0
1591  el0 = el(1)
1592  conit = 0.5d0/(nq+2)
1593  ddn = dvnrms(n,savf,ewt)/tesco(1,l)
1594  exdn = 1.0d0/l
1595  rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
1596  rh = min(rhdn,1.0d0)
1597  iredo = 3
1598  IF (h .EQ. hold) GO TO 660
1599  rh = min(rh,abs(h/hold))
1600  h = hold
1601  GO TO 100
1602  50 CONTINUE
1603 C ------------------------------------------------------------
1604 C DCFOD IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS
1605 C FOR THE CURRENT METH. THEN THE EL VECTOR AND RELATED
1606 C CONSTANTS ARE RESET WHENEVER THE ORDER NQ IS CHANGED, OR AT
1607 C THE START OF THE PROBLEM.
1608 C ------------------------------------------------------------
1609  CALL dcfod(meth,elco,tesco)
1610  60 CONTINUE
1611  70 CONTINUE
1612 C BEGIN BLOCK PERMITTING ...EXITS TO 680
1613  DO 80 i = 1, l
1614  el(i) = elco(i,nq)
1615  80 CONTINUE
1616  nqnyh = nq*nyh
1617  rc = rc*el(1)/el0
1618  el0 = el(1)
1619  conit = 0.5d0/(nq+2)
1620  GO TO (90,660,160), iret
1621 C ---------------------------------------------------------
1622 C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST
1623 C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH
1624 C IS SET TO L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT
1625 C MANY STEPS, UNLESS FORCED BY A CONVERGENCE OR ERROR TEST
1626 C FAILURE.
1627 C ---------------------------------------------------------
1628  90 CONTINUE
1629  IF (h .EQ. hold) GO TO 160
1630  rh = h/hold
1631  h = hold
1632  iredo = 3
1633  100 CONTINUE
1634  110 CONTINUE
1635  rh = min(rh,rmax)
1636  rh = rh/max(1.0d0,abs(h)*hmxi*rh)
1637  r = 1.0d0
1638  DO 130 j = 2, l
1639  r = r*rh
1640  DO 120 i = 1, n
1641  yh(i,j) = yh(i,j)*r
1642  120 CONTINUE
1643  130 CONTINUE
1644  h = h*rh
1645  rc = rc*rh
1646  ialth = l
1647  IF (iredo .NE. 0) GO TO 150
1648  rmax = 10.0d0
1649  r = 1.0d0/tesco(2,nqu)
1650  DO 140 i = 1, n
1651  acor(i) = acor(i)*r
1652  140 CONTINUE
1653 C ...............EXIT
1654  GO TO 690
1655  150 CONTINUE
1656 C ------------------------------------------------------
1657 C THIS SECTION COMPUTES THE PREDICTED VALUES BY
1658 C EFFECTIVELY MULTIPLYING THE YH ARRAY BY THE PASCAL
1659 C TRIANGLE MATRIX. RC IS THE RATIO OF NEW TO OLD
1660 C VALUES OF THE COEFFICIENT H*EL(1). WHEN RC DIFFERS
1661 C FROM 1 BY MORE THAN 30 PERCENT, IPUP IS SET TO MITER
1662 C TO FORCE DPJAC TO BE CALLED, IF A JACOBIAN IS
1663 C INVOLVED. IN ANY CASE, DPJAC IS CALLED AT LEAST
1664 C EVERY 20-TH STEP.
1665 C ------------------------------------------------------
1666  160 CONTINUE
1667  170 CONTINUE
1668 C BEGIN BLOCK PERMITTING ...EXITS TO 610
1669 C BEGIN BLOCK PERMITTING ...EXITS TO 490
1670  IF (abs(rc-1.0d0) .GT. 0.3d0) ipup = miter
1671  IF (nst .GE. nstepj + 20) ipup = miter
1672  tn = tn + h
1673  i1 = nqnyh + 1
1674  DO 190 jb = 1, nq
1675  i1 = i1 - nyh
1676  DO 180 i = i1, nqnyh
1677  yh1(i) = yh1(i) + yh1(i+nyh)
1678  180 CONTINUE
1679  190 CONTINUE
1680  ksteps = ksteps + 1
1681 C ---------------------------------------------
1682 C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A
1683 C CONVERGENCE TEST IS MADE ON THE R.M.S. NORM
1684 C OF EACH CORRECTION, WEIGHTED BY THE ERROR
1685 C WEIGHT VECTOR EWT. THE SUM OF THE
1686 C CORRECTIONS IS ACCUMULATED IN THE VECTOR
1687 C ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE
1688 C CORRECTOR LOOP.
1689 C ---------------------------------------------
1690  200 CONTINUE
1691  m = 0
1692  DO 210 i = 1, n
1693  y(i) = yh(i,1)
1694  210 CONTINUE
1695  CALL df(tn,y,savf,rpar,ipar)
1696  nfe = nfe + 1
1697  IF (ipup .LE. 0) GO TO 220
1698 C ---------------------------------------
1699 C IF INDICATED, THE MATRIX P = I -
1700 C H*EL(1)*J IS REEVALUATED AND
1701 C PREPROCESSED BEFORE STARTING THE
1702 C CORRECTOR ITERATION. IPUP IS SET TO 0
1703 C AS AN INDICATOR THAT THIS HAS BEEN
1704 C DONE.
1705 C ---------------------------------------
1706  ipup = 0
1707  rc = 1.0d0
1708  nstepj = nst
1709  crate = 0.7d0
1710  CALL dpjac(neq,y,yh,nyh,ewt,acor,savf,
1711  1 wm,iwm,df,djac,rpar,ipar)
1712 C ......EXIT
1713  IF (ier .NE. 0) GO TO 440
1714  220 CONTINUE
1715  DO 230 i = 1, n
1716  acor(i) = 0.0d0
1717  230 CONTINUE
1718  240 CONTINUE
1719  IF (miter .NE. 0) GO TO 270
1720 C ------------------------------------
1721 C IN THE CASE OF FUNCTIONAL
1722 C ITERATION, UPDATE Y DIRECTLY FROM
1723 C THE RESULT OF THE LAST FUNCTION
1724 C EVALUATION.
1725 C ------------------------------------
1726  DO 250 i = 1, n
1727  savf(i) = h*savf(i) - yh(i,2)
1728  y(i) = savf(i) - acor(i)
1729  250 CONTINUE
1730  del = dvnrms(n,y,ewt)
1731  DO 260 i = 1, n
1732  y(i) = yh(i,1) + el(1)*savf(i)
1733  acor(i) = savf(i)
1734  260 CONTINUE
1735  GO TO 300
1736  270 CONTINUE
1737 C ------------------------------------
1738 C IN THE CASE OF THE CHORD METHOD,
1739 C COMPUTE THE CORRECTOR ERROR, AND
1740 C SOLVE THE LINEAR SYSTEM WITH THAT
1741 C AS RIGHT-HAND SIDE AND P AS
1742 C COEFFICIENT MATRIX.
1743 C ------------------------------------
1744  DO 280 i = 1, n
1745  y(i) = h*savf(i)
1746  1 - (yh(i,2) + acor(i))
1747  280 CONTINUE
1748  CALL dslvs(wm,iwm,y,savf)
1749 C ......EXIT
1750  IF (ier .NE. 0) GO TO 430
1751  del = dvnrms(n,y,ewt)
1752  DO 290 i = 1, n
1753  acor(i) = acor(i) + y(i)
1754  y(i) = yh(i,1) + el(1)*acor(i)
1755  290 CONTINUE
1756  300 CONTINUE
1757 C ---------------------------------------
1758 C TEST FOR CONVERGENCE. IF M.GT.0, AN
1759 C ESTIMATE OF THE CONVERGENCE RATE
1760 C CONSTANT IS STORED IN CRATE, AND THIS
1761 C IS USED IN THE TEST.
1762 C ---------------------------------------
1763  IF (m .NE. 0)
1764  1 crate = max(0.2d0*crate,del/delp)
1765  dcon = del*min(1.0d0,1.5d0*crate)
1766  1 /(tesco(2,nq)*conit)
1767  IF (dcon .GT. 1.0d0) GO TO 420
1768 C ------------------------------------
1769 C THE CORRECTOR HAS CONVERGED. IPUP
1770 C IS SET TO -1 IF MITER .NE. 0, TO
1771 C SIGNAL THAT THE JACOBIAN INVOLVED
1772 C MAY NEED UPDATING LATER. THE LOCAL
1773 C ERROR TEST IS MADE AND CONTROL
1774 C PASSES TO STATEMENT 500 IF IT
1775 C FAILS.
1776 C ------------------------------------
1777  IF (miter .NE. 0) ipup = -1
1778  IF (m .EQ. 0) dsm = del/tesco(2,nq)
1779  IF (m .GT. 0)
1780  1 dsm = dvnrms(n,acor,ewt)
1781  2 /tesco(2,nq)
1782  IF (dsm .GT. 1.0d0) GO TO 380
1783 C BEGIN BLOCK
1784 C PERMITTING ...EXITS TO 360
1785 C ------------------------------
1786 C AFTER A SUCCESSFUL STEP,
1787 C UPDATE THE YH ARRAY.
1788 C CONSIDER CHANGING H IF IALTH
1789 C = 1. OTHERWISE DECREASE
1790 C IALTH BY 1. IF IALTH IS THEN
1791 C 1 AND NQ .LT. MAXORD, THEN
1792 C ACOR IS SAVED FOR USE IN A
1793 C POSSIBLE ORDER INCREASE ON
1794 C THE NEXT STEP. IF A CHANGE
1795 C IN H IS CONSIDERED, AN
1796 C INCREASE OR DECREASE IN ORDER
1797 C BY ONE IS CONSIDERED ALSO. A
1798 C CHANGE IN H IS MADE ONLY IF
1799 C IT IS BY A FACTOR OF AT LEAST
1800 C 1.1. IF NOT, IALTH IS SET TO
1801 C 3 TO PREVENT TESTING FOR THAT
1802 C MANY STEPS.
1803 C ------------------------------
1804  kflag = 0
1805  iredo = 0
1806  nst = nst + 1
1807  hu = h
1808  nqu = nq
1809  DO 320 j = 1, l
1810  DO 310 i = 1, n
1811  yh(i,j) = yh(i,j)
1812  1 + el(j)
1813  2 *acor(i)
1814  310 CONTINUE
1815  320 CONTINUE
1816  ialth = ialth - 1
1817  IF (ialth .NE. 0) GO TO 340
1818 C ---------------------------
1819 C REGARDLESS OF THE SUCCESS
1820 C OR FAILURE OF THE STEP,
1821 C FACTORS RHDN, RHSM, AND
1822 C RHUP ARE COMPUTED, BY
1823 C WHICH H COULD BE
1824 C MULTIPLIED AT ORDER NQ -
1825 C 1, ORDER NQ, OR ORDER NQ +
1826 C 1, RESPECTIVELY. IN THE
1827 C CASE OF FAILURE, RHUP =
1828 C 0.0 TO AVOID AN ORDER
1829 C INCREASE. THE LARGEST OF
1830 C THESE IS DETERMINED AND
1831 C THE NEW ORDER CHOSEN
1832 C ACCORDINGLY. IF THE ORDER
1833 C IS TO BE INCREASED, WE
1834 C COMPUTE ONE ADDITIONAL
1835 C SCALED DERIVATIVE.
1836 C ---------------------------
1837  rhup = 0.0d0
1838 C .....................EXIT
1839  IF (l .EQ. lmax) GO TO 490
1840  DO 330 i = 1, n
1841  savf(i) = acor(i)
1842  1 - yh(i,lmax)
1843  330 CONTINUE
1844  dup = dvnrms(n,savf,ewt)
1845  1 /tesco(3,nq)
1846  exup = 1.0d0/(l+1)
1847  rhup = 1.0d0
1848  1 /(1.4d0*dup**exup
1849  2 + 0.0000014d0)
1850 C .....................EXIT
1851  GO TO 490
1852  340 CONTINUE
1853 C ...EXIT
1854  IF (ialth .GT. 1) GO TO 360
1855 C ...EXIT
1856  IF (l .EQ. lmax) GO TO 360
1857  DO 350 i = 1, n
1858  yh(i,lmax) = acor(i)
1859  350 CONTINUE
1860  360 CONTINUE
1861  r = 1.0d0/tesco(2,nqu)
1862  DO 370 i = 1, n
1863  acor(i) = acor(i)*r
1864  370 CONTINUE
1865 C .................................EXIT
1866  GO TO 690
1867  380 CONTINUE
1868 C ------------------------------------
1869 C THE ERROR TEST FAILED. KFLAG KEEPS
1870 C TRACK OF MULTIPLE FAILURES.
1871 C RESTORE TN AND THE YH ARRAY TO
1872 C THEIR PREVIOUS VALUES, AND PREPARE
1873 C TO TRY THE STEP AGAIN. COMPUTE THE
1874 C OPTIMUM STEP SIZE FOR THIS OR ONE
1875 C LOWER ORDER. AFTER 2 OR MORE
1876 C FAILURES, H IS FORCED TO DECREASE
1877 C BY A FACTOR OF 0.2 OR LESS.
1878 C ------------------------------------
1879  kflag = kflag - 1
1880  tn = told
1881  i1 = nqnyh + 1
1882  DO 400 jb = 1, nq
1883  i1 = i1 - nyh
1884  DO 390 i = i1, nqnyh
1885  yh1(i) = yh1(i) - yh1(i+nyh)
1886  390 CONTINUE
1887  400 CONTINUE
1888  rmax = 2.0d0
1889  IF (abs(h) .GT. hmin*1.00001d0)
1890  1 GO TO 410
1891 C ---------------------------------
1892 C ALL RETURNS ARE MADE THROUGH
1893 C THIS SECTION. H IS SAVED IN
1894 C HOLD TO ALLOW THE CALLER TO
1895 C CHANGE H ON THE NEXT STEP.
1896 C ---------------------------------
1897  kflag = -1
1898 C .................................EXIT
1899  GO TO 690
1900  410 CONTINUE
1901 C ...............EXIT
1902  IF (kflag .LE. -3) GO TO 610
1903  iredo = 2
1904  rhup = 0.0d0
1905 C ............EXIT
1906  GO TO 490
1907  420 CONTINUE
1908  m = m + 1
1909 C ...EXIT
1910  IF (m .EQ. 3) GO TO 430
1911 C ...EXIT
1912  IF (m .GE. 2 .AND. del .GT. 2.0d0*delp)
1913  1 GO TO 430
1914  delp = del
1915  CALL df(tn,y,savf,rpar,ipar)
1916  nfe = nfe + 1
1917  GO TO 240
1918  430 CONTINUE
1919 C ------------------------------------------
1920 C THE CORRECTOR ITERATION FAILED TO
1921 C CONVERGE IN 3 TRIES. IF MITER .NE. 0 AND
1922 C THE JACOBIAN IS OUT OF DATE, DPJAC IS
1923 C CALLED FOR THE NEXT TRY. OTHERWISE THE
1924 C YH ARRAY IS RETRACTED TO ITS VALUES
1925 C BEFORE PREDICTION, AND H IS REDUCED, IF
1926 C POSSIBLE. IF H CANNOT BE REDUCED OR 10
1927 C FAILURES HAVE OCCURRED, EXIT WITH KFLAG =
1928 C -2.
1929 C ------------------------------------------
1930 C ...EXIT
1931  IF (ipup .EQ. 0) GO TO 440
1932  ipup = miter
1933  GO TO 200
1934  440 CONTINUE
1935  tn = told
1936  ncf = ncf + 1
1937  rmax = 2.0d0
1938  i1 = nqnyh + 1
1939  DO 460 jb = 1, nq
1940  i1 = i1 - nyh
1941  DO 450 i = i1, nqnyh
1942  yh1(i) = yh1(i) - yh1(i+nyh)
1943  450 CONTINUE
1944  460 CONTINUE
1945  IF (abs(h) .GT. hmin*1.00001d0) GO TO 470
1946  kflag = -2
1947 C ........................EXIT
1948  GO TO 690
1949  470 CONTINUE
1950  IF (ncf .NE. 10) GO TO 480
1951  kflag = -2
1952 C ........................EXIT
1953  GO TO 690
1954  480 CONTINUE
1955  rh = 0.25d0
1956  ipup = miter
1957  iredo = 1
1958 C .........EXIT
1959  GO TO 650
1960  490 CONTINUE
1961  exsm = 1.0d0/l
1962  rhsm = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
1963  rhdn = 0.0d0
1964  IF (nq .EQ. 1) GO TO 500
1965  ddn = dvnrms(n,yh(1,l),ewt)/tesco(1,nq)
1966  exdn = 1.0d0/nq
1967  rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
1968  500 CONTINUE
1969  IF (rhsm .GE. rhup) GO TO 550
1970  IF (rhup .LE. rhdn) GO TO 540
1971  newq = l
1972  rh = rhup
1973  IF (rh .GE. 1.1d0) GO TO 520
1974  ialth = 3
1975  r = 1.0d0/tesco(2,nqu)
1976  DO 510 i = 1, n
1977  acor(i) = acor(i)*r
1978  510 CONTINUE
1979 C ...........................EXIT
1980  GO TO 690
1981  520 CONTINUE
1982  r = el(l)/l
1983  DO 530 i = 1, n
1984  yh(i,newq+1) = acor(i)*r
1985  530 CONTINUE
1986  nq = newq
1987  l = nq + 1
1988  iret = 2
1989 C ..................EXIT
1990  GO TO 680
1991  540 CONTINUE
1992  GO TO 580
1993  550 CONTINUE
1994  IF (rhsm .LT. rhdn) GO TO 580
1995  newq = nq
1996  rh = rhsm
1997  IF (kflag .EQ. 0 .AND. rh .LT. 1.1d0)
1998  1 GO TO 560
1999  IF (kflag .LE. -2) rh = min(rh,0.2d0)
2000 C ------------------------------------------
2001 C IF THERE IS A CHANGE OF ORDER, RESET NQ,
2002 C L, AND THE COEFFICIENTS. IN ANY CASE H
2003 C IS RESET ACCORDING TO RH AND THE YH ARRAY
2004 C IS RESCALED. THEN EXIT FROM 680 IF THE
2005 C STEP WAS OK, OR REDO THE STEP OTHERWISE.
2006 C ------------------------------------------
2007 C ............EXIT
2008  IF (newq .EQ. nq) GO TO 650
2009  nq = newq
2010  l = nq + 1
2011  iret = 2
2012 C ..................EXIT
2013  GO TO 680
2014  560 CONTINUE
2015  ialth = 3
2016  r = 1.0d0/tesco(2,nqu)
2017  DO 570 i = 1, n
2018  acor(i) = acor(i)*r
2019  570 CONTINUE
2020 C .....................EXIT
2021  GO TO 690
2022  580 CONTINUE
2023  newq = nq - 1
2024  rh = rhdn
2025  IF (kflag .LT. 0 .AND. rh .GT. 1.0d0) rh = 1.0d0
2026  IF (kflag .EQ. 0 .AND. rh .LT. 1.1d0) GO TO 590
2027  IF (kflag .LE. -2) rh = min(rh,0.2d0)
2028 C ---------------------------------------------
2029 C IF THERE IS A CHANGE OF ORDER, RESET NQ, L,
2030 C AND THE COEFFICIENTS. IN ANY CASE H IS
2031 C RESET ACCORDING TO RH AND THE YH ARRAY IS
2032 C RESCALED. THEN EXIT FROM 680 IF THE STEP
2033 C WAS OK, OR REDO THE STEP OTHERWISE.
2034 C ---------------------------------------------
2035 C .........EXIT
2036  IF (newq .EQ. nq) GO TO 650
2037  nq = newq
2038  l = nq + 1
2039  iret = 2
2040 C ...............EXIT
2041  GO TO 680
2042  590 CONTINUE
2043  ialth = 3
2044  r = 1.0d0/tesco(2,nqu)
2045  DO 600 i = 1, n
2046  acor(i) = acor(i)*r
2047  600 CONTINUE
2048 C ..................EXIT
2049  GO TO 690
2050  610 CONTINUE
2051 C ---------------------------------------------------
2052 C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES
2053 C HAVE OCCURRED. IF 10 FAILURES HAVE OCCURRED, EXIT
2054 C WITH KFLAG = -1. IT IS ASSUMED THAT THE
2055 C DERIVATIVES THAT HAVE ACCUMULATED IN THE YH ARRAY
2056 C HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST
2057 C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO
2058 C 1. THEN H IS REDUCED BY A FACTOR OF 10, AND THE
2059 C STEP IS RETRIED, UNTIL IT SUCCEEDS OR H REACHES
2060 C HMIN.
2061 C ---------------------------------------------------
2062  IF (kflag .NE. -10) GO TO 620
2063 C ------------------------------------------------
2064 C ALL RETURNS ARE MADE THROUGH THIS SECTION. H
2065 C IS SAVED IN HOLD TO ALLOW THE CALLER TO CHANGE
2066 C H ON THE NEXT STEP.
2067 C ------------------------------------------------
2068  kflag = -1
2069 C ..................EXIT
2070  GO TO 690
2071  620 CONTINUE
2072  rh = 0.1d0
2073  rh = max(hmin/abs(h),rh)
2074  h = h*rh
2075  DO 630 i = 1, n
2076  y(i) = yh(i,1)
2077  630 CONTINUE
2078  CALL df(tn,y,savf,rpar,ipar)
2079  nfe = nfe + 1
2080  DO 640 i = 1, n
2081  yh(i,2) = h*savf(i)
2082  640 CONTINUE
2083  ipup = miter
2084  ialth = 5
2085 C ......EXIT
2086  IF (nq .NE. 1) GO TO 670
2087  GO TO 170
2088  650 CONTINUE
2089  660 CONTINUE
2090  rh = max(rh,hmin/abs(h))
2091  GO TO 110
2092  670 CONTINUE
2093  nq = 1
2094  l = 2
2095  iret = 3
2096  680 CONTINUE
2097  GO TO 70
2098  690 CONTINUE
2099  hold = h
2100  jstart = 1
2101  RETURN
2102 C ----------------------- END OF SUBROUTINE DSTOD
2103 C -----------------------
max
#define max(a, b)
Definition: cascade.c:32
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
min
#define min(a, b)
Definition: cascade.c:31
dvnrms
double precision function dvnrms(N, V, W)
Definition: ddebdf.f:2252
djac
subroutine djac(x, u, pd, nrowpd, rpar, nev)
Definition: subspace.f:163
dcfod
subroutine dcfod(METH, ELCO, TESCO)
Definition: ddebdf.f:2107
dpjac
subroutine dpjac(NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, DF, DJAC, RPAR, IPAR)
Definition: ddebdf.f:2401
dslvs
subroutine dslvs(WM, IWM, X, TEM)
Definition: ddebdf.f:2627

◆ dvnrms()

double precision function dvnrms ( integer  N,
double precision, dimension(*)  V,
double precision, dimension(*)  W 
)
2252 C***BEGIN PROLOGUE DVNRMS
2253 C***SUBSIDIARY
2254 C***PURPOSE Subsidiary to DDEBDF
2255 C***LIBRARY SLATEC
2256 C***TYPE DOUBLE PRECISION (VNWRMS-S, DVNRMS-D)
2257 C***AUTHOR (UNKNOWN)
2258 C***DESCRIPTION
2259 C
2260 C DVNRMS computes a weighted root-mean-square vector norm for the
2261 C integrator package DDEBDF.
2262 C
2263 C***SEE ALSO DDEBDF
2264 C***ROUTINES CALLED (NONE)
2265 C***REVISION HISTORY (YYMMDD)
2266 C 820301 DATE WRITTEN
2267 C 890531 Changed all specific intrinsics to generic. (WRB)
2268 C 890831 Modified array declarations. (WRB)
2269 C 890911 Removed unnecessary intrinsics. (WRB)
2270 C 891214 Prologue converted to Version 4.0 format. (BAB)
2271 C 900328 Added TYPE section. (WRB)
2272 C***END PROLOGUE DVNRMS
2273  INTEGER i, n
2274  DOUBLE PRECISION sum, v, w
2275  dimension v(*),w(*)
2276 C***FIRST EXECUTABLE STATEMENT DVNRMS
2277  sum = 0.0d0
2278  DO 10 i = 1, n
2279  sum = sum + (v(i)/w(i))**2
2280  10 CONTINUE
2281  dvnrms = sqrt(sum/n)
2282  RETURN
2283 C ----------------------- END OF FUNCTION DVNRMS
2284 C ------------------------
dvnrms
double precision function dvnrms(N, V, W)
Definition: ddebdf.f:2252
Hosted by OpenAircraft.com, (Michigan UAV, LLC)