1
SUBROUTINE DDASKR (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
2
* IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL,
5
C***BEGIN PROLOGUE DDASKR
6
C***REVISION HISTORY (YYMMDD)
8
C 021105 Changed yprime argument in DRCHEK calls to YPRIME.
9
C 021217 Modified error return for zeros found too close together.
10
C 021217 Added root direction output in JROOT.
11
C 031201 stuck root masking
12
c 040615 Removing Hmin requirement
13
c 040615 Separating the error message of Singular Jacobian in DDASID
16
C***KEYWORDS DIFFERENTIAL/ALGEBRAIC, BACKWARD DIFFERENTIATION FORMULAS,
17
C IMPLICIT DIFFERENTIAL SYSTEMS, KRYLOV ITERATION
18
C***AUTHORS Linda R. Petzold, Peter N. Brown, Alan C. Hindmarsh, and
20
C Center for Computational Sciences & Engineering, L-316
21
C Lawrence Livermore National Laboratory
24
C***PURPOSE This code solves a system of differential/algebraic
25
C equations of the form
27
C using a combination of Backward Differentiation Formula
28
C (BDF) methods and a choice of two linear system solution
29
C methods: direct (dense or band) or Krylov (iterative).
30
C This version is in double precision.
31
C-----------------------------------------------------------------------
36
C IMPLICIT DOUBLE PRECISION(A-H,O-Z)
37
C INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR(*)
38
C DOUBLE PRECISION T, Y(*), YPRIME(*), TOUT, RTOL(*), ATOL(*),
40
C EXTERNAL RES, JAC, PSOL, RT
42
C CALL DDASKR (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
43
C * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL,
46
C Quantities which may be altered by the code are:
47
C T, Y(*), YPRIME(*), INFO(1), RTOL, ATOL, IDID, RWORK(*), IWORK(*)
52
C RES:EXT This is the name of a subroutine which you
53
C provide to define the residual function G(t,y,y')
54
C of the differential/algebraic system.
56
C NEQ:IN This is the number of equations in the system.
58
C T:INOUT This is the current value of the independent
61
C Y(*):INOUT This array contains the solution components at T.
63
C YPRIME(*):INOUT This array contains the derivatives of the solution
66
C TOUT:IN This is a point at which a solution is desired.
68
C INFO(N):IN This is an integer array used to communicate details
69
C of how the solution is to be carried out, such as
70
C tolerance type, matrix structure, step size and
71
C order limits, and choice of nonlinear system method.
72
C N must be at least 20.
74
C RTOL,ATOL:INOUT These quantities represent absolute and relative
75
C error tolerances (on local error) which you provide
76
C to indicate how accurately you wish the solution to
77
C be computed. You may choose them to be both scalars
78
C or else both arrays of length NEQ.
80
C IDID:OUT This integer scalar is an indicator reporting what
81
C the code did. You must monitor this variable to
82
C decide what action to take next.
84
C RWORK:WORK A real work array of length LRW which provides the
85
C code with needed storage space.
87
C LRW:IN The length of RWORK.
89
C IWORK:WORK An integer work array of length LIW which provides
90
C the code with needed storage space.
92
C LIW:IN The length of IWORK.
94
C RPAR,IPAR:IN These are real and integer parameter arrays which
95
C you can use for communication between your calling
96
C program and the RES, JAC, and PSOL subroutines.
98
C JAC:EXT This is the name of a subroutine which you may
99
C provide (optionally) for calculating Jacobian
100
C (partial derivative) data involved in solving linear
101
C systems within DDASKR.
103
C PSOL:EXT This is the name of a subroutine which you must
104
C provide for solving linear systems if you selected
105
C a Krylov method. The purpose of PSOL is to solve
106
C linear systems involving a left preconditioner P.
108
C RT:EXT This is the name of the subroutine for defining
109
C constraint functions Ri(T,Y,Y')) whose roots are
110
C desired during the integration. This name must be
111
C declared external in the calling program.
113
C NRT:IN This is the number of constraint functions
114
C Ri(T,Y,Y'). If there are no constraints, set
115
C NRT = 0, and pass a dummy name for RT.
117
C JROOT:OUT This is an integer array of length NRT for output
118
C of root information.
122
C The DDASKR solver uses the backward differentiation formulas of
123
C orders one through five to solve a system of the form G(t,y,y') = 0
124
C for y = Y and y' = YPRIME. Values for Y and YPRIME at the initial
125
C time must be given as input. These values should be consistent,
126
C that is, if T, Y, YPRIME are the given initial values, they should
127
C satisfy G(T,Y,YPRIME) = 0. However, if consistent values are not
128
C known, in many cases you can have DDASKR solve for them -- see
129
C INFO(11). (This and other options are described in detail below.)
131
C Normally, DDASKR solves the system from T to TOUT. It is easy to
132
C continue the solution to get results at additional TOUT. This is
133
C the interval mode of operation. Intermediate results can also be
134
C obtained easily by specifying INFO(3).
136
C On each step taken by DDASKR, a sequence of nonlinear algebraic
137
C systems arises. These are solved by one of two types of
139
C * a Newton iteration with a direct method for the linear
140
C systems involved (INFO(12) = 0), or
141
C * a Newton iteration with a preconditioned Krylov iterative
142
C method for the linear systems involved (INFO(12) = 1).
144
C The direct method choices are dense and band matrix solvers,
145
C with either a user-supplied or an internal difference quotient
146
C Jacobian matrix, as specified by INFO(5) and INFO(6).
147
C In the band case, INFO(6) = 1, you must supply half-bandwidths
148
C in IWORK(1) and IWORK(2).
150
C The Krylov method is the Generalized Minimum Residual (GMRES)
151
C method, in either complete or incomplete form, and with
152
C scaling and preconditioning. The method is implemented
153
C in an algorithm called SPIGMR. Certain options in the Krylov
154
C method case are specified by INFO(13) and INFO(15).
156
C If the Krylov method is chosen, you may supply a pair of routines,
157
C JAC and PSOL, to apply preconditioning to the linear system.
158
C If the system is A*x = b, the matrix is A = dG/dY + CJ*dG/dYPRIME
159
C (of order NEQ). This system can then be preconditioned in the form
160
C (P-inverse)*A*x = (P-inverse)*b, with left preconditioner P.
161
C (DDASKR does not allow right preconditioning.)
162
C Then the Krylov method is applied to this altered, but equivalent,
163
C linear system, hopefully with much better performance than without
164
C preconditioning. (In addition, a diagonal scaling matrix based on
165
C the tolerances is also introduced into the altered system.)
167
C The JAC routine evaluates any data needed for solving systems
168
C with coefficient matrix P, and PSOL carries out that solution.
169
C In any case, in order to improve convergence, you should try to
170
C make P approximate the matrix A as much as possible, while keeping
171
C the system P*x = b reasonably easy and inexpensive to solve for x,
174
C While integrating the given DAE system, DDASKR also searches for
175
C roots of the given constraint functions Ri(T,Y,Y') given by RT.
176
C If DDASKR detects a sign change in any Ri(T,Y,Y'), it will return
177
C the intermediate value of T and Y for which Ri(T,Y,Y') = 0.
181
C------INPUT - WHAT TO DO ON THE FIRST CALL TO DDASKR-------------------
184
C The first call of the code is defined to be the start of each new
185
C problem. Read through the descriptions of all the following items,
186
C provide sufficient storage space for designated arrays, set
187
C appropriate variables for the initialization of the problem, and
188
C give information about how you want the problem to be solved.
191
C RES -- Provide a subroutine of the form
193
C SUBROUTINE RES (T, Y, YPRIME, CJ, DELTA, IRES, RPAR, IPAR)
195
C to define the system of differential/algebraic
196
C equations which is to be solved. For the given values
197
C of T, Y and YPRIME, the subroutine should return
198
C the residual of the differential/algebraic system
199
C DELTA = G(T,Y,YPRIME)
200
C DELTA is a vector of length NEQ which is output from RES.
202
C Subroutine RES must not alter T, Y, YPRIME, or CJ.
203
C You must declare the name RES in an EXTERNAL
204
C statement in your program that calls DDASKR.
205
C You must dimension Y, YPRIME, and DELTA in RES.
207
C The input argument CJ can be ignored, or used to rescale
208
C constraint equations in the system (see Ref. 2, p. 145).
209
C Note: In this respect, DDASKR is not downward-compatible
210
C with DDASSL, which does not have the RES argument CJ.
212
C IRES is an integer flag which is always equal to zero
213
C on input. Subroutine RES should alter IRES only if it
214
C encounters an illegal value of Y or a stop condition.
215
C Set IRES = -1 if an input value is illegal, and DDASKR
216
C will try to solve the problem without getting IRES = -1.
217
C If IRES = -2, DDASKR will return control to the calling
218
C program with IDID = -11.
220
C RPAR and IPAR are real and integer parameter arrays which
221
C you can use for communication between your calling program
222
C and subroutine RES. They are not altered by DDASKR. If you
223
C do not need RPAR or IPAR, ignore these parameters by treat-
224
C ing them as dummy arguments. If you do choose to use them,
225
C dimension them in your calling program and in RES as arrays
226
C of appropriate length.
228
C NEQ -- Set it to the number of equations in the system (NEQ .GE. 1).
230
C T -- Set it to the initial point of the integration. (T must be
233
C Y(*) -- Set this array to the initial values of the NEQ solution
234
C components at the initial point. You must dimension Y of
235
C length at least NEQ in your calling program.
237
C YPRIME(*) -- Set this array to the initial values of the NEQ first
238
C derivatives of the solution components at the initial
239
C point. You must dimension YPRIME at least NEQ in your
242
C TOUT - Set it to the first point at which a solution is desired.
243
C You cannot take TOUT = T. Integration either forward in T
244
C (TOUT .GT. T) or backward in T (TOUT .LT. T) is permitted.
246
C The code advances the solution from T to TOUT using step
247
C sizes which are automatically selected so as to achieve the
248
C desired accuracy. If you wish, the code will return with the
249
C solution and its derivative at intermediate steps (the
250
C intermediate-output mode) so that you can monitor them,
251
C but you still must provide TOUT in accord with the basic
254
C The first step taken by the code is a critical one because
255
C it must reflect how fast the solution changes near the
256
C initial point. The code automatically selects an initial
257
C step size which is practically always suitable for the
258
C problem. By using the fact that the code will not step past
259
C TOUT in the first step, you could, if necessary, restrict the
260
C length of the initial step.
262
C For some problems it may not be permissible to integrate
263
C past a point TSTOP, because a discontinuity occurs there
264
C or the solution or its derivative is not defined beyond
265
C TSTOP. When you have declared a TSTOP point (see INFO(4)
266
C and RWORK(1)), you have told the code not to integrate past
267
C TSTOP. In this case any tout beyond TSTOP is invalid input.
269
C INFO(*) - Use the INFO array to give the code more details about
270
C how you want your problem solved. This array should be
271
C dimensioned of length 20, though DDASKR uses only the
272
C first 15 entries. You must respond to all of the following
273
C items, which are arranged as questions. The simplest use
274
C of DDASKR corresponds to setting all entries of INFO to 0.
276
C INFO(1) - This parameter enables the code to initialize itself.
277
C You must set it to indicate the start of every new
280
C **** Is this the first call for this problem ...
281
C yes - set INFO(1) = 0
282
C no - not applicable here.
283
C See below for continuation calls. ****
285
C INFO(2) - How much accuracy you want of your solution
286
C is specified by the error tolerances RTOL and ATOL.
287
C The simplest use is to take them both to be scalars.
288
C To obtain more flexibility, they can both be arrays.
289
C The code must be told your choice.
291
C **** Are both error tolerances RTOL, ATOL scalars ...
292
C yes - set INFO(2) = 0
293
C and input scalars for both RTOL and ATOL
294
C no - set INFO(2) = 1
295
C and input arrays for both RTOL and ATOL ****
297
C INFO(3) - The code integrates from T in the direction of TOUT
298
C by steps. If you wish, it will return the computed
299
C solution and derivative at the next intermediate step
300
C (the intermediate-output mode) or TOUT, whichever comes
301
C first. This is a good way to proceed if you want to
302
C see the behavior of the solution. If you must have
303
C solutions at a great many specific TOUT points, this
304
C code will compute them efficiently.
306
C **** Do you want the solution only at
307
C TOUT (and not at the next intermediate step) ...
308
C yes - set INFO(3) = 0 (interval-output mode)
309
C no - set INFO(3) = 1 (intermediate-output mode) ****
311
C INFO(4) - To handle solutions at a great many specific
312
C values TOUT efficiently, this code may integrate past
313
C TOUT and interpolate to obtain the result at TOUT.
314
C Sometimes it is not possible to integrate beyond some
315
C point TSTOP because the equation changes there or it is
316
C not defined past TSTOP. Then you must tell the code
317
C this stop condition.
319
C **** Can the integration be carried out without any
320
C restrictions on the independent variable T ...
321
C yes - set INFO(4) = 0
322
C no - set INFO(4) = 1
323
C and define the stopping point TSTOP by
324
C setting RWORK(1) = TSTOP ****
326
C INFO(5) - used only when INFO(12) = 0 (direct methods).
327
C To solve differential/algebraic systems you may wish
328
C to use a matrix of partial derivatives of the
329
C system of differential equations. If you do not
330
C provide a subroutine to evaluate it analytically (see
331
C description of the item JAC in the call list), it will
332
C be approximated by numerical differencing in this code.
333
C Although it is less trouble for you to have the code
334
C compute partial derivatives by numerical differencing,
335
C the solution will be more reliable if you provide the
336
C derivatives via JAC. Usually numerical differencing is
337
C more costly than evaluating derivatives in JAC, but
338
C sometimes it is not - this depends on your problem.
340
C **** Do you want the code to evaluate the partial deriv-
341
C atives automatically by numerical differences ...
342
C yes - set INFO(5) = 0
343
C no - set INFO(5) = 1
344
C and provide subroutine JAC for evaluating the
345
C matrix of partial derivatives ****
347
C INFO(6) - used only when INFO(12) = 0 (direct methods).
348
C DDASKR will perform much better if the matrix of
349
C partial derivatives, dG/dY + CJ*dG/dYPRIME (here CJ is
350
C a scalar determined by DDASKR), is banded and the code
351
C is told this. In this case, the storage needed will be
352
C greatly reduced, numerical differencing will be performed
353
C much cheaper, and a number of important algorithms will
354
C execute much faster. The differential equation is said
355
C to have half-bandwidths ML (lower) and MU (upper) if
356
C equation i involves only unknowns Y(j) with
357
C i-ML .le. j .le. i+MU .
358
C For all i=1,2,...,NEQ. Thus, ML and MU are the widths
359
C of the lower and upper parts of the band, respectively,
360
C with the main diagonal being excluded. If you do not
361
C indicate that the equation has a banded matrix of partial
362
C derivatives the code works with a full matrix of NEQ**2
363
C elements (stored in the conventional way). Computations
364
C with banded matrices cost less time and storage than with
365
C full matrices if 2*ML+MU .lt. NEQ. If you tell the
366
C code that the matrix of partial derivatives has a banded
367
C structure and you want to provide subroutine JAC to
368
C compute the partial derivatives, then you must be careful
369
C to store the elements of the matrix in the special form
370
C indicated in the description of JAC.
372
C **** Do you want to solve the problem using a full (dense)
373
C matrix (and not a special banded structure) ...
374
C yes - set INFO(6) = 0
375
C no - set INFO(6) = 1
376
C and provide the lower (ML) and upper (MU)
377
C bandwidths by setting
381
C INFO(7) - You can specify a maximum (absolute value of)
382
C stepsize, so that the code will avoid passing over very
385
C **** Do you want the code to decide on its own the maximum
387
C yes - set INFO(7) = 0
388
C no - set INFO(7) = 1
389
C and define HMAX by setting
390
C RWORK(2) = HMAX ****
392
C INFO(8) - Differential/algebraic problems may occasionally
393
C suffer from severe scaling difficulties on the first
394
C step. If you know a great deal about the scaling of
395
C your problem, you can help to alleviate this problem
396
C by specifying an initial stepsize H0.
398
C **** Do you want the code to define its own initial
400
C yes - set INFO(8) = 0
401
C no - set INFO(8) = 1
402
C and define H0 by setting
405
C INFO(9) - If storage is a severe problem, you can save some
406
C storage by restricting the maximum method order MAXORD.
407
C The default value is 5. For each order decrease below 5,
408
C the code requires NEQ fewer locations, but it is likely
409
C to be slower. In any case, you must have
410
C 1 .le. MAXORD .le. 5.
411
C **** Do you want the maximum order to default to 5 ...
412
C yes - set INFO(9) = 0
413
C no - set INFO(9) = 1
414
C and define MAXORD by setting
415
C IWORK(3) = MAXORD ****
417
C INFO(10) - If you know that certain components of the
418
C solutions to your equations are always nonnegative
419
C (or nonpositive), it may help to set this
420
C parameter. There are three options that are
422
C 1. To have constraint checking only in the initial
423
C condition calculation.
424
C 2. To enforce nonnegativity in Y during the integration.
425
C 3. To enforce both options 1 and 2.
427
C When selecting option 2 or 3, it is probably best to try
428
C the code without using this option first, and only use
429
C this option if that does not work very well.
431
C **** Do you want the code to solve the problem without
432
C invoking any special inequality constraints ...
433
C yes - set INFO(10) = 0
434
C no - set INFO(10) = 1 to have option 1 enforced
435
C no - set INFO(10) = 2 to have option 2 enforced
436
C no - set INFO(10) = 3 to have option 3 enforced ****
438
C If you have specified INFO(10) = 1 or 3, then you
439
C will also need to identify how each component of Y
440
C in the initial condition calculation is constrained.
442
C IWORK(40+I) = +1 if Y(I) must be .GE. 0,
443
C IWORK(40+I) = +2 if Y(I) must be .GT. 0,
444
C IWORK(40+I) = -1 if Y(I) must be .LE. 0, while
445
C IWORK(40+I) = -2 if Y(I) must be .LT. 0, while
446
C IWORK(40+I) = 0 if Y(I) is not constrained.
448
C INFO(11) - DDASKR normally requires the initial T, Y, and
449
C YPRIME to be consistent. That is, you must have
450
C G(T,Y,YPRIME) = 0 at the initial T. If you do not know
451
C the initial conditions precisely, in some cases
452
C DDASKR may be able to compute it.
454
C Denoting the differential variables in Y by Y_d
455
C and the algebraic variables by Y_a, DDASKR can solve
456
C one of two initialization problems:
457
C 1. Given Y_d, calculate Y_a and Y'_d, or
458
C 2. Given Y', calculate Y.
459
C In either case, initial values for the given
460
C components are input, and initial guesses for
461
C the unknown components must also be provided as input.
463
C **** Are the initial T, Y, YPRIME consistent ...
465
C yes - set INFO(11) = 0
466
C no - set INFO(11) = 1 to calculate option 1 above,
467
C or set INFO(11) = 2 to calculate option 2 ****
469
C If you have specified INFO(11) = 1, then you
470
C will also need to identify which are the
471
C differential and which are the algebraic
472
C components (algebraic components are components
473
C whose derivatives do not appear explicitly
474
C in the function G(T,Y,YPRIME)). You must set:
475
C IWORK(LID+I) = +1 if Y(I) is a differential variable
476
C IWORK(LID+I) = -1 if Y(I) is an algebraic variable,
477
C where LID = 40 if INFO(10) = 0 or 2 and LID = 40+NEQ
478
C if INFO(10) = 1 or 3.
480
C INFO(12) - Except for the addition of the RES argument CJ,
481
C DDASKR by default is downward-compatible with DDASSL,
482
C which uses only direct (dense or band) methods to solve
483
C the linear systems involved. You must set INFO(12) to
484
C indicate whether you want the direct methods or the
485
C Krylov iterative method.
486
C **** Do you want DDASKR to use standard direct methods
487
C (dense or band) or the Krylov (iterative) method ...
488
C direct methods - set INFO(12) = 0.
489
C Krylov method - set INFO(12) = 1,
490
C and check the settings of INFO(13) and INFO(15).
492
C INFO(13) - used when INFO(12) = 1 (Krylov methods).
493
C DDASKR uses scalars MAXL, KMP, NRMAX, and EPLI for the
494
C iterative solution of linear systems. INFO(13) allows
495
C you to override the default values of these parameters.
496
C These parameters and their defaults are as follows:
497
C MAXL = maximum number of iterations in the SPIGMR
498
C algorithm (MAXL .le. NEQ). The default is
500
C KMP = number of vectors on which orthogonalization is
501
C done in the SPIGMR algorithm. The default is
502
C KMP = MAXL, which corresponds to complete GMRES
503
C iteration, as opposed to the incomplete form.
504
C NRMAX = maximum number of restarts of the SPIGMR
505
C algorithm per nonlinear iteration. The default is
507
C EPLI = convergence test constant in SPIGMR algorithm.
508
C The default is EPLI = 0.05.
509
C Note that the length of RWORK depends on both MAXL
510
C and KMP. See the definition of LRW below.
511
C **** Are MAXL, KMP, and EPLI to be given their
513
C yes - set INFO(13) = 0
514
C no - set INFO(13) = 1,
515
C and set all of the following:
516
C IWORK(24) = MAXL (1 .le. MAXL .le. NEQ)
517
C IWORK(25) = KMP (1 .le. KMP .le. MAXL)
518
C IWORK(26) = NRMAX (NRMAX .ge. 0)
519
C RWORK(10) = EPLI (0 .lt. EPLI .lt. 1.0) ****
521
C INFO(14) - used with INFO(11) > 0 (initial condition
522
C calculation is requested). In this case, you may
523
C request control to be returned to the calling program
524
C immediately after the initial condition calculation,
525
C before proceeding to the integration of the system
526
C (e.g. to examine the computed Y and YPRIME).
527
C If this is done, and if the initialization succeeded
528
C (IDID = 4), you should reset INFO(11) to 0 for the
529
C next call, to prevent the solver from repeating the
530
C initialization (and to avoid an infinite loop).
531
C **** Do you want to proceed to the integration after
532
C the initial condition calculation is done ...
533
C yes - set INFO(14) = 0
534
C no - set INFO(14) = 1 ****
536
C INFO(15) - used when INFO(12) = 1 (Krylov methods).
537
C When using preconditioning in the Krylov method,
538
C you must supply a subroutine, PSOL, which solves the
539
C associated linear systems using P.
540
C The usage of DDASKR is simpler if PSOL can carry out
541
C the solution without any prior calculation of data.
542
C However, if some partial derivative data is to be
543
C calculated in advance and used repeatedly in PSOL,
544
C then you must supply a JAC routine to do this,
545
C and set INFO(15) to indicate that JAC is to be called
546
C for this purpose. For example, P might be an
547
C approximation to a part of the matrix A which can be
548
C calculated and LU-factored for repeated solutions of
549
C the preconditioner system. The arrays WP and IWP
550
C (described under JAC and PSOL) can be used to
551
C communicate data between JAC and PSOL.
552
C **** Does PSOL operate with no prior preparation ...
553
C yes - set INFO(15) = 0 (no JAC routine)
554
C no - set INFO(15) = 1
555
C and supply a JAC routine to evaluate and
556
C preprocess any required Jacobian data. ****
558
C INFO(16) - option to exclude algebraic variables from
560
C **** Do you wish to control errors locally on
561
C all the variables...
562
C yes - set INFO(16) = 0
563
C no - set INFO(16) = 1
564
C If you have specified INFO(16) = 1, then you
565
C will also need to identify which are the
566
C differential and which are the algebraic
567
C components (algebraic components are components
568
C whose derivatives do not appear explicitly
569
C in the function G(T,Y,YPRIME)). You must set:
570
C IWORK(LID+I) = +1 if Y(I) is a differential
572
C IWORK(LID+I) = -1 if Y(I) is an algebraic
574
C where LID = 40 if INFO(10) = 0 or 2 and
575
C LID = 40 + NEQ if INFO(10) = 1 or 3.
577
C INFO(17) - used when INFO(11) > 0 (DDASKR is to do an
578
C initial condition calculation).
579
C DDASKR uses several heuristic control quantities in the
580
C initial condition calculation. They have default values,
581
C but can also be set by the user using INFO(17).
582
C These parameters and their defaults are as follows:
583
C MXNIT = maximum number of Newton iterations
584
C per Jacobian or preconditioner evaluation.
586
C MXNIT = 5 in the direct case (INFO(12) = 0), and
587
C MXNIT = 15 in the Krylov case (INFO(12) = 1).
588
C MXNJ = maximum number of Jacobian or preconditioner
589
C evaluations. The default is:
590
C MXNJ = 6 in the direct case (INFO(12) = 0), and
591
C MXNJ = 2 in the Krylov case (INFO(12) = 1).
592
C MXNH = maximum number of values of the artificial
593
C stepsize parameter H to be tried if INFO(11) = 1.
594
C The default is MXNH = 5.
595
C NOTE: the maximum number of Newton iterations
596
C allowed in all is MXNIT*MXNJ*MXNH if INFO(11) = 1,
597
C and MXNIT*MXNJ if INFO(11) = 2.
598
C LSOFF = flag to turn off the linesearch algorithm
599
C (LSOFF = 0 means linesearch is on, LSOFF = 1 means
600
C it is turned off). The default is LSOFF = 0.
601
C STPTOL = minimum scaled step in linesearch algorithm.
602
C The default is STPTOL = (unit roundoff)**(2/3).
603
C EPINIT = swing factor in the Newton iteration convergence
604
C test. The test is applied to the residual vector,
605
C premultiplied by the approximate Jacobian (in the
606
C direct case) or the preconditioner (in the Krylov
607
C case). For convergence, the weighted RMS norm of
608
C this vector (scaled by the error weights) must be
609
C less than EPINIT*EPCON, where EPCON = .33 is the
610
C analogous test constant used in the time steps.
611
C The default is EPINIT = .01.
612
C **** Are the initial condition heuristic controls to be
613
C given their default values...
614
C yes - set INFO(17) = 0
615
C no - set INFO(17) = 1,
616
C and set all of the following:
617
C IWORK(32) = MXNIT (.GT. 0)
618
C IWORK(33) = MXNJ (.GT. 0)
619
C IWORK(34) = MXNH (.GT. 0)
620
C IWORK(35) = LSOFF ( = 0 or 1)
621
C RWORK(14) = STPTOL (.GT. 0.0)
622
C RWORK(15) = EPINIT (.GT. 0.0) ****
624
C INFO(18) - option to get extra printing in initial condition
626
C **** Do you wish to have extra printing...
627
C no - set INFO(18) = 0
628
C yes - set INFO(18) = 1 for minimal printing, or
629
C set INFO(18) = 2 for full printing.
630
C If you have specified INFO(18) .ge. 1, data
631
C will be printed with the error handler routines.
632
C To print to a non-default unit number L, include
633
C the line CALL XSETUN(L) in your program. ****
635
C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL)
636
C error tolerances to tell the code how accurately you
637
C want the solution to be computed. They must be defined
638
C as variables because the code may change them.
639
C you have two choices --
640
C Both RTOL and ATOL are scalars (INFO(2) = 0), or
641
C both RTOL and ATOL are vectors (INFO(2) = 1).
642
C In either case all components must be non-negative.
644
C The tolerances are used by the code in a local error
645
C test at each step which requires roughly that
646
C abs(local error in Y(i)) .le. EWT(i) ,
647
C where EWT(i) = RTOL*abs(Y(i)) + ATOL is an error weight
648
C quantity, for each vector component.
649
C (More specifically, a root-mean-square norm is used to
650
C measure the size of vectors, and the error test uses the
651
C magnitude of the solution at the beginning of the step.)
653
C The true (global) error is the difference between the
654
C true solution of the initial value problem and the
655
C computed approximation. Practically all present day
656
C codes, including this one, control the local error at
657
C each step and do not even attempt to control the global
660
C Usually, but not always, the true accuracy of
661
C the computed Y is comparable to the error tolerances.
662
C This code will usually, but not always, deliver a more
663
C accurate solution if you reduce the tolerances and
664
C integrate again. By comparing two such solutions you
665
C can get a fairly reliable idea of the true error in the
666
C solution at the larger tolerances.
668
C Setting ATOL = 0. results in a pure relative error test
669
C on that component. Setting RTOL = 0. results in a pure
670
C absolute error test on that component. A mixed test
671
C with non-zero RTOL and ATOL corresponds roughly to a
672
C relative error test when the solution component is
673
C much bigger than ATOL and to an absolute error test
674
C when the solution component is smaller than the
677
C The code will not attempt to compute a solution at an
678
C accuracy unreasonable for the machine being used. It
679
C will advise you if you ask for too much accuracy and
680
C inform you as to the maximum accuracy it believes
683
C RWORK(*) -- a real work array, which should be dimensioned in your
684
C calling program with a length equal to the value of
687
C LRW -- Set it to the declared length of the RWORK array. The
688
C minimum length depends on the options you have selected,
689
C given by a base value plus additional storage as
692
C If INFO(12) = 0 (standard direct method), the base value
693
C is BASE = 60 + max(MAXORD+4,7)*NEQ + 3*NRT.
694
C The default value is MAXORD = 5 (see INFO(9)). With the
695
C default MAXORD, BASE = 60 + 9*NEQ + 3*NRT.
696
C Additional storage must be added to the base value for
697
C any or all of the following options:
698
C If INFO(6) = 0 (dense matrix), add NEQ**2.
699
C If INFO(6) = 1 (banded matrix), then:
700
C if INFO(5) = 0, add (2*ML+MU+1)*NEQ
701
C + 2*[NEQ/(ML+MU+1) + 1], and
702
C if INFO(5) = 1, add (2*ML+MU+1)*NEQ.
703
C If INFO(16) = 1, add NEQ.
705
C If INFO(12) = 1 (Krylov method), the base value is
706
C BASE = 60 + (MAXORD+5)*NEQ + 3*NRT
707
C + [MAXL + 3 + min(1,MAXL-KMP)]*NEQ
708
C + (MAXL+3)*MAXL + 1 + LENWP.
709
C See PSOL for description of LENWP. The default values
710
C are: MAXORD = 5 (see INFO(9)), MAXL = min(5,NEQ) and
711
C KMP = MAXL (see INFO(13)). With these default values,
712
C BASE = 101 + 18*NEQ + 3*NRT + LENWP.
713
C Additional storage must be added to the base value for
714
C the following option:
715
C If INFO(16) = 1, add NEQ.
718
C IWORK(*) -- an integer work array, which should be dimensioned in
719
C your calling program with a length equal to the value
720
C of LIW (or greater).
722
C LIW -- Set it to the declared length of the IWORK array. The
723
C minimum length depends on the options you have selected,
724
C given by a base value plus additions as described below.
726
C If INFO(12) = 0 (standard direct method), the base value
727
C is BASE = 40 + NEQ.
728
C IF INFO(10) = 1 or 3, add NEQ to the base value.
729
C If INFO(11) = 1 or INFO(16) =1, add NEQ to the base value.
731
C If INFO(12) = 1 (Krylov method), the base value is
732
C BASE = 40 + LENIWP. See PSOL for description of LENIWP.
733
C If INFO(10) = 1 or 3, add NEQ to the base value.
734
C If INFO(11) = 1 or INFO(16) =1, add NEQ to the base value.
735
C >> Due to introduction of Mask in DASKR, NRT has been added
739
C RPAR, IPAR -- These are arrays of double precision and integer type,
740
C respectively, which are available for you to use
741
C for communication between your program that calls
742
C DDASKR and the RES subroutine (and the JAC and PSOL
743
C subroutines). They are not altered by DDASKR.
744
C If you do not need RPAR or IPAR, ignore these
745
C parameters by treating them as dummy arguments.
746
C If you do choose to use them, dimension them in
747
C your calling program and in RES (and in JAC and PSOL)
748
C as arrays of appropriate length.
750
C JAC -- This is the name of a routine that you may supply
751
C (optionally) that relates to the Jacobian matrix of the
752
C nonlinear system that the code must solve at each T step.
753
C The role of JAC (and its call sequence) depends on whether
754
C a direct (INFO(12) = 0) or Krylov (INFO(12) = 1) method
757
C **** INFO(12) = 0 (direct methods):
758
C If you are letting the code generate partial derivatives
759
C numerically (INFO(5) = 0), then JAC can be absent
760
C (or perhaps a dummy routine to satisfy the loader).
761
C Otherwise you must supply a JAC routine to compute
762
C the matrix A = dG/dY + CJ*dG/dYPRIME. It must have
765
C SUBROUTINE JAC (T, Y, YPRIME, PD, CJ, RPAR, IPAR)
767
C The JAC routine must dimension Y, YPRIME, and PD (and RPAR
768
C and IPAR if used). CJ is a scalar which is input to JAC.
769
C For the given values of T, Y, and YPRIME, the JAC routine
770
C must evaluate the nonzero elements of the matrix A, and
771
C store these values in the array PD. The elements of PD are
772
C set to zero before each call to JAC, so that only nonzero
773
C elements need to be defined.
774
C The way you store the elements into the PD array depends
775
C on the structure of the matrix indicated by INFO(6).
776
C *** INFO(6) = 0 (full or dense matrix) ***
777
C Give PD a first dimension of NEQ. When you evaluate the
778
C nonzero partial derivatives of equation i (i.e. of G(i))
779
C with respect to component j (of Y and YPRIME), you must
780
C store the element in PD according to
781
C PD(i,j) = dG(i)/dY(j) + CJ*dG(i)/dYPRIME(j).
782
C *** INFO(6) = 1 (banded matrix with half-bandwidths ML, MU
783
C as described under INFO(6)) ***
784
C Give PD a first dimension of 2*ML+MU+1. When you
785
C evaluate the nonzero partial derivatives of equation i
786
C (i.e. of G(i)) with respect to component j (of Y and
787
C YPRIME), you must store the element in PD according to
788
C IROW = i - j + ML + MU + 1
789
C PD(IROW,j) = dG(i)/dY(j) + CJ*dG(i)/dYPRIME(j).
791
C **** INFO(12) = 1 (Krylov method):
792
C If you are not calculating Jacobian data in advance for use
793
C in PSOL (INFO(15) = 0), JAC can be absent (or perhaps a
794
C dummy routine to satisfy the loader). Otherwise, you may
795
C supply a JAC routine to compute and preprocess any parts of
796
C of the Jacobian matrix A = dG/dY + CJ*dG/dYPRIME that are
797
C involved in the preconditioner matrix P.
798
C It is to have the form
800
C SUBROUTINE JAC (RES, IRES, NEQ, T, Y, YPRIME, REWT, SAVR,
801
C WK, H, CJ, WP, IWP, IER, RPAR, IPAR)
803
C The JAC routine must dimension Y, YPRIME, REWT, SAVR, WK,
804
C and (if used) WP, IWP, RPAR, and IPAR.
805
C The Y, YPRIME, and SAVR arrays contain the current values
806
C of Y, YPRIME, and the residual G, respectively.
807
C The array WK is work space of length NEQ.
808
C H is the step size. CJ is a scalar, input to JAC, that is
809
C normally proportional to 1/H. REWT is an array of
810
C reciprocal error weights, 1/EWT(i), where EWT(i) is
811
C RTOL*abs(Y(i)) + ATOL (unless you supplied routine DDAWTS
812
C instead), for use in JAC if needed. For example, if JAC
813
C computes difference quotient approximations to partial
814
C derivatives, the REWT array may be useful in setting the
815
C increments used. The JAC routine should do any
816
C factorization operations called for, in preparation for
817
C solving linear systems in PSOL. The matrix P should
818
C be an approximation to the Jacobian,
819
C A = dG/dY + CJ*dG/dYPRIME.
821
C WP and IWP are real and integer work arrays which you may
822
C use for communication between your JAC routine and your
823
C PSOL routine. These may be used to store elements of the
824
C preconditioner P, or related matrix data (such as factored
825
C forms). They are not altered by DDASKR.
826
C If you do not need WP or IWP, ignore these parameters by
827
C treating them as dummy arguments. If you do use them,
828
C dimension them appropriately in your JAC and PSOL routines.
829
C See the PSOL description for instructions on setting
830
C the lengths of WP and IWP.
832
C On return, JAC should set the error flag IER as follows..
833
C IER = 0 if JAC was successful,
834
C IER .ne. 0 if JAC was unsuccessful (e.g. if Y or YPRIME
835
C was illegal, or a singular matrix is found).
836
C (If IER .ne. 0, a smaller stepsize will be tried.)
837
C IER = 0 on entry to JAC, so need be reset only on a failure.
838
C If RES is used within JAC, then a nonzero value of IRES will
839
C override any nonzero value of IER (see the RES description).
841
C Regardless of the method type, subroutine JAC must not
842
C alter T, Y(*), YPRIME(*), H, CJ, or REWT(*).
843
C You must declare the name JAC in an EXTERNAL statement in
844
C your program that calls DDASKR.
846
C PSOL -- This is the name of a routine you must supply if you have
847
C selected a Krylov method (INFO(12) = 1) with preconditioning.
848
C In the direct case (INFO(12) = 0), PSOL can be absent
849
C (a dummy routine may have to be supplied to satisfy the
850
C loader). Otherwise, you must provide a PSOL routine to
851
C solve linear systems arising from preconditioning.
852
C When supplied with INFO(12) = 1, the PSOL routine is to
855
C SUBROUTINE PSOL (NEQ, T, Y, YPRIME, SAVR, WK, CJ, WGHT,
856
C WP, IWP, B, EPLIN, IER, RPAR, IPAR)
858
C The PSOL routine must solve linear systems of the form
859
C P*x = b where P is the left preconditioner matrix.
861
C The right-hand side vector b is in the B array on input, and
862
C PSOL must return the solution vector x in B.
863
C The Y, YPRIME, and SAVR arrays contain the current values
864
C of Y, YPRIME, and the residual G, respectively.
866
C Work space required by JAC and/or PSOL, and space for data to
867
C be communicated from JAC to PSOL is made available in the form
868
C of arrays WP and IWP, which are parts of the RWORK and IWORK
869
C arrays, respectively. The lengths of these real and integer
870
C work spaces WP and IWP must be supplied in LENWP and LENIWP,
871
C respectively, as follows..
872
C IWORK(27) = LENWP = length of real work space WP
873
C IWORK(28) = LENIWP = length of integer work space IWP.
875
C WK is a work array of length NEQ for use by PSOL.
876
C CJ is a scalar, input to PSOL, that is normally proportional
877
C to 1/H (H = stepsize). If the old value of CJ
878
C (at the time of the last JAC call) is needed, it must have
879
C been saved by JAC in WP.
881
C WGHT is an array of weights, to be used if PSOL uses an
882
C iterative method and performs a convergence test. (In terms
883
C of the argument REWT to JAC, WGHT is REWT/sqrt(NEQ).)
884
C If PSOL uses an iterative method, it should use EPLIN
885
C (a heuristic parameter) as the bound on the weighted norm of
886
C the residual for the computed solution. Specifically, the
887
C residual vector R should satisfy
888
C SQRT (SUM ( (R(i)*WGHT(i))**2 ) ) .le. EPLIN
890
C PSOL must not alter NEQ, T, Y, YPRIME, SAVR, CJ, WGHT, EPLIN.
892
C On return, PSOL should set the error flag IER as follows..
893
C IER = 0 if PSOL was successful,
894
C IER .lt. 0 if an unrecoverable error occurred, meaning
895
C control will be passed to the calling routine,
896
C IER .gt. 0 if a recoverable error occurred, meaning that
897
C the step will be retried with the same step size
898
C but with a call to JAC to update necessary data,
899
C unless the Jacobian data is current, in which case
900
C the step will be retried with a smaller step size.
901
C IER = 0 on entry to PSOL so need be reset only on a failure.
903
C You must declare the name PSOL in an EXTERNAL statement in
904
C your program that calls DDASKR.
906
C RT -- This is the name of the subroutine for defining the vector
907
C R(T,Y,Y') of constraint functions Ri(T,Y,Y'), whose roots
908
C are desired during the integration. It is to have the form
909
C SUBROUTINE RT(NEQ, T, Y, YP, NRT, RVAL, RPAR, IPAR)
910
C DIMENSION Y(NEQ), YP(NEQ), RVAL(NRT),
911
C where NEQ, T, Y and NRT are INPUT, and the array RVAL is
912
C output. NEQ, T, Y, and YP have the same meaning as in the
913
C RES routine, and RVAL is an array of length NRT.
914
C For i = 1,...,NRT, this routine is to load into RVAL(i) the
915
C value at (T,Y,Y') of the i-th constraint function Ri(T,Y,Y').
916
C DDASKR will find roots of the Ri of odd multiplicity
917
C (that is, sign changes) as they occur during the integration.
918
C RT must be declared EXTERNAL in the calling program.
920
C CAUTION.. Because of numerical errors in the functions Ri
921
C due to roundoff and integration error, DDASKR may return
922
C false roots, or return the same root at two or more nearly
923
C equal values of T. If such false roots are suspected,
924
C the user should consider smaller error tolerances and/or
925
C higher precision in the evaluation of the Ri.
927
C If a root of some Ri defines the end of the problem,
928
C the input to DDASKR should nevertheless allow
929
C integration to a point slightly past that root, so
930
C that DDASKR can locate the root by interpolation.
932
C NRT -- The number of constraint functions Ri(T,Y,Y'). If there are
933
C no constraints, set NRT = 0 and pass a dummy name for RT.
935
C JROOT -- This is an integer array of length NRT, used only for output.
936
C On a return where one or more roots were found (IDID = 5),
937
C JROOT(i) = 1 or -1 if Ri(T,Y,Y') has a root at T, and
938
C JROOT(i) = 0 if not. If nonzero, JROOT(i) shows the direction
939
C of the sign change in Ri in the direction of integration:
940
C JROOT(i) = 1 means Ri changed from negative to positive.
941
C JROOT(i) = -1 means Ri changed from positive to negative.
944
C OPTIONALLY REPLACEABLE SUBROUTINE:
946
C DDASKR uses a weighted root-mean-square norm to measure the
947
C size of various error vectors. The weights used in this norm
948
C are set in the following subroutine:
950
C SUBROUTINE DDAWTS1 (NEQ, IWT, RTOL, ATOL, Y, EWT, RPAR, IPAR)
951
C DIMENSION RTOL(*), ATOL(*), Y(*), EWT(*), RPAR(*), IPAR(*)
953
C A DDAWTS routine has been included with DDASKR which sets the
954
C weights according to
955
C EWT(I) = RTOL*ABS(Y(I)) + ATOL
956
C in the case of scalar tolerances (IWT = 0) or
957
C EWT(I) = RTOL(I)*ABS(Y(I)) + ATOL(I)
958
C in the case of array tolerances (IWT = 1). (IWT is INFO(2).)
959
C In some special cases, it may be appropriate for you to define
960
C your own error weights by writing a subroutine DDAWTS to be
961
C called instead of the version supplied. However, this should
962
C be attempted only after careful thought and consideration.
963
C If you supply this routine, you may use the tolerances and Y
964
C as appropriate, but do not overwrite these variables. You
965
C may also use RPAR and IPAR to communicate data as appropriate.
966
C ***Note: Aside from the values of the weights, the choice of
967
C norm used in DDASKR (weighted root-mean-square) is not subject
968
C to replacement by the user. In this respect, DDASKR is not
969
C downward-compatible with the original DDASSL solver (in which
970
C the norm routine was optionally user-replaceable).
973
C------OUTPUT - AFTER ANY RETURN FROM DDASKR----------------------------
975
C The principal aim of the code is to return a computed solution at
976
C T = TOUT, although it is also possible to obtain intermediate
977
C results along the way. To find out whether the code achieved its
978
C goal or if the integration process was interrupted before the task
979
C was completed, you must check the IDID parameter.
982
C T -- The output value of T is the point to which the solution
983
C was successfully advanced.
985
C Y(*) -- contains the computed solution approximation at T.
987
C YPRIME(*) -- contains the computed derivative approximation at T.
989
C IDID -- reports what the code did, described as follows:
991
C *** TASK COMPLETED ***
992
C Reported by positive values of IDID
994
C IDID = 1 -- A step was successfully taken in the
995
C interval-output mode. The code has not
998
C IDID = 2 -- The integration to TSTOP was successfully
999
C completed (T = TSTOP) by stepping exactly to TSTOP.
1001
C IDID = 3 -- The integration to TOUT was successfully
1002
C completed (T = TOUT) by stepping past TOUT.
1003
C Y(*) and YPRIME(*) are obtained by interpolation.
1005
C IDID = 4 -- The initial condition calculation, with
1006
C INFO(11) > 0, was successful, and INFO(14) = 1.
1007
C No integration steps were taken, and the solution
1008
C is not considered to have been started.
1010
C IDID = 5 -- The integration was successfully completed
1011
C by finding one or more roots of R(T,Y,Y') at T.
1012
C IDID = 6 -- The integration was successfully completed
1013
C by finding A ROOT, LIFTED FROM ZERO.
1016
C *** TASK INTERRUPTED ***
1017
C Reported by negative values of IDID
1019
C IDID = -1 -- A large amount of work has been expended
1020
C (about 500 steps).
1022
C IDID = -2 -- The error tolerances are too stringent.
1024
C IDID = -3 -- The local error test cannot be satisfied
1025
C because you specified a zero component in ATOL
1026
C and the corresponding computed solution component
1027
C is zero. Thus, a pure relative error test is
1028
C impossible for this component.
1030
C IDID = -5 -- There were repeated failures in the evaluation
1031
C or processing of the preconditioner (in JAC).
1033
C IDID = -6 -- DDASKR had repeated error test failures on the
1034
C last attempted step.
1036
C IDID = -7 -- The nonlinear system solver in the time
1037
C integration could not converge.
1039
C IDID = -8 -- The matrix of partial derivatives appears
1040
C to be singular (direct method).
1042
C IDID = -9 -- The nonlinear system solver in the integration
1043
C failed to achieve convergence, and there were
1044
C repeated error test failures in this step.
1046
C IDID =-10 -- The nonlinear system solver in the integration
1047
C failed to achieve convergence because IRES was
1050
C IDID =-11 -- IRES = -2 was encountered and control is
1051
C being returned to the calling program.
1053
C IDID =-12 -- DDASKR failed to compute the initial Y, YPRIME.
1055
C IDID =-13 -- An unrecoverable error was encountered inside
1056
C the user's PSOL routine, and control is being
1057
C returned to the calling program.
1059
C IDID =-14 -- The Krylov linear system solver could not
1060
C achieve convergence.
1063
C IDID =-15,..,-32 -- Not applicable for this code.
1066
C *** TASK TERMINATED ***
1067
C reported by the value of IDID=-33
1069
C IDID = -33 -- The code has encountered trouble from which
1070
C it cannot recover. A message is printed
1071
C explaining the trouble and control is returned
1072
C to the calling program. For example, this occurs
1073
C when invalid input is detected.
1075
C RTOL, ATOL -- these quantities remain unchanged except when
1076
C IDID = -2. In this case, the error tolerances have been
1077
C increased by the code to values which are estimated to
1078
C be appropriate for continuing the integration. However,
1079
C the reported solution at T was obtained using the input
1080
C values of RTOL and ATOL.
1082
C RWORK, IWORK -- contain information which is usually of no interest
1083
C to the user but necessary for subsequent calls.
1084
C However, you may be interested in the performance data
1085
C listed below. These quantities are accessed in RWORK
1086
C and IWORK but have internal mnemonic names, as follows..
1088
C RWORK(3)--contains H, the step size h to be attempted
1091
C RWORK(4)--contains TN, the current value of the
1092
C independent variable, i.e. the farthest point
1093
C integration has reached. This will differ
1094
C from T if interpolation has been performed
1097
C RWORK(7)--contains HOLD, the stepsize used on the last
1098
C successful step. If INFO(11) = INFO(14) = 1,
1099
C this contains the value of H used in the
1100
C initial condition calculation.
1102
C IWORK(7)--contains K, the order of the method to be
1103
C attempted on the next step.
1105
C IWORK(8)--contains KOLD, the order of the method used
1108
C IWORK(11)--contains NST, the number of steps (in T)
1111
C IWORK(12)--contains NRE, the number of calls to RES
1114
C IWORK(13)--contains NJE, the number of calls to JAC so
1115
C far (Jacobian or preconditioner evaluations).
1117
C IWORK(14)--contains NETF, the total number of error test
1120
C IWORK(15)--contains NCFN, the total number of nonlinear
1121
C convergence failures so far (includes counts
1122
C of singular iteration matrix or singular
1125
C IWORK(16)--contains NCFL, the number of convergence
1126
C failures of the linear iteration so far.
1128
C IWORK(17)--contains LENIW, the length of IWORK actually
1129
C required. This is defined on normal returns
1130
C and on an illegal input return for
1131
C insufficient storage.
1133
C IWORK(18)--contains LENRW, the length of RWORK actually
1134
C required. This is defined on normal returns
1135
C and on an illegal input return for
1136
C insufficient storage.
1138
C IWORK(19)--contains NNI, the total number of nonlinear
1139
C iterations so far (each of which calls a
1142
C IWORK(20)--contains NLI, the total number of linear
1143
C (Krylov) iterations so far.
1145
C IWORK(21)--contains NPS, the number of PSOL calls so
1146
C far, for preconditioning solve operations or
1147
C for solutions with the user-supplied method.
1149
C IWORK(36)--contains the total number of calls to the
1150
C constraint function routine RT so far.
1152
C Note: The various counters in IWORK do not include
1153
C counts during a prior call made with INFO(11) > 0 and
1157
C------INPUT - WHAT TO DO TO CONTINUE THE INTEGRATION -----------------
1158
C (CALLS AFTER THE FIRST)
1160
C This code is organized so that subsequent calls to continue the
1161
C integration involve little (if any) additional effort on your
1162
C part. You must monitor the IDID parameter in order to determine
1165
C Recalling that the principal task of the code is to integrate
1166
C from T to TOUT (the interval mode), usually all you will need
1167
C to do is specify a new TOUT upon reaching the current TOUT.
1169
C Do not alter any quantity not specifically permitted below. In
1170
C particular do not alter NEQ, T, Y(*), YPRIME(*), RWORK(*),
1171
C IWORK(*), or the differential equation in subroutine RES. Any
1172
C such alteration constitutes a new problem and must be treated
1173
C as such, i.e. you must start afresh.
1175
C You cannot change from array to scalar error control or vice
1176
C versa (INFO(2)), but you can change the size of the entries of
1177
C RTOL or ATOL. Increasing a tolerance makes the equation easier
1178
C to integrate. Decreasing a tolerance will make the equation
1179
C harder to integrate and should generally be avoided.
1181
C You can switch from the intermediate-output mode to the
1182
C interval mode (INFO(3)) or vice versa at any time.
1184
C If it has been necessary to prevent the integration from going
1185
C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the
1186
C code will not integrate to any TOUT beyond the currently
1187
C specified TSTOP. Once TSTOP has been reached, you must change
1188
C the value of TSTOP or set INFO(4) = 0. You may change INFO(4)
1189
C or TSTOP at any time but you must supply the value of TSTOP in
1190
C RWORK(1) whenever you set INFO(4) = 1.
1192
C Do not change INFO(5), INFO(6), INFO(12-17) or their associated
1193
C IWORK/RWORK locations unless you are going to restart the code.
1195
C *** FOLLOWING A COMPLETED TASK ***
1198
C IDID = 1, call the code again to continue the integration
1199
C another step in the direction of TOUT.
1201
C IDID = 2 or 3, define a new TOUT and call the code again.
1202
C TOUT must be different from T. You cannot change
1203
C the direction of integration without restarting.
1205
C IDID = 4, reset INFO(11) = 0 and call the code again to begin
1206
C the integration. (If you leave INFO(11) > 0 and
1207
C INFO(14) = 1, you may generate an infinite loop.)
1208
C In this situation, the next call to DDASKR is
1209
C considered to be the first call for the problem,
1210
C in that all initializations are done.
1212
C IDID = 5, call the code again to continue the integration in the
1213
C direction of TOUT. You may change the functions
1214
C Ri defined by RT after a return with IDID = 5, but
1215
C the number of constraint functions NRT must remain
1216
C the same. If you wish to change the functions in
1217
C RES or in RT, then you must restart the code.
1219
C *** FOLLOWING AN INTERRUPTED TASK ***
1221
C To show the code that you realize the task was interrupted and
1222
C that you want to continue, you must take appropriate action and
1226
C IDID = -1, the code has taken about 500 steps. If you want to
1227
C continue, set INFO(1) = 1 and call the code again.
1228
C An additional 500 steps will be allowed.
1231
C IDID = -2, the error tolerances RTOL, ATOL have been increased
1232
C to values the code estimates appropriate for
1233
C continuing. You may want to change them yourself.
1234
C If you are sure you want to continue with relaxed
1235
C error tolerances, set INFO(1) = 1 and call the code
1238
C IDID = -3, a solution component is zero and you set the
1239
C corresponding component of ATOL to zero. If you
1240
C are sure you want to continue, you must first alter
1241
C the error criterion to use positive values of ATOL
1242
C for those components corresponding to zero solution
1243
C components, then set INFO(1) = 1 and call the code
1246
C IDID = -4 --- cannot occur with this code.
1248
C IDID = -5, your JAC routine failed with the Krylov method. Check
1249
C for errors in JAC and restart the integration.
1251
C IDID = -6, repeated error test failures occurred on the last
1252
C attempted step in DDASKR. A singularity in the
1253
C solution may be present. If you are absolutely
1254
C certain you want to continue, you should restart
1255
C the integration. (Provide initial values of Y and
1256
C YPRIME which are consistent.)
1258
C IDID = -7, repeated convergence test failures occurred on the last
1259
C attempted step in DDASKR. An inaccurate or ill-
1260
C conditioned Jacobian or preconditioner may be the
1261
C problem. If you are absolutely certain you want
1262
C to continue, you should restart the integration.
1265
C IDID = -8, the matrix of partial derivatives is singular, with
1266
C the use of direct methods. Some of your equations
1267
C may be redundant. DDASKR cannot solve the problem
1268
C as stated. It is possible that the redundant
1269
C equations could be removed, and then DDASKR could
1270
C solve the problem. It is also possible that a
1271
C solution to your problem either does not exist
1274
C IDID = -9, DDASKR had multiple convergence test failures, preceded
1275
C by multiple error test failures, on the last
1276
C attempted step. It is possible that your problem is
1277
C ill-posed and cannot be solved using this code. Or,
1278
C there may be a discontinuity or a singularity in the
1279
C solution. If you are absolutely certain you want to
1280
C continue, you should restart the integration.
1282
C IDID = -10, DDASKR had multiple convergence test failures
1283
C because IRES was equal to -1. If you are
1284
C absolutely certain you want to continue, you
1285
C should restart the integration.
1287
C IDID = -11, there was an unrecoverable error (IRES = -2) from RES
1288
C inside the nonlinear system solver. Determine the
1289
C cause before trying again.
1291
C IDID = -12, DDASKR failed to compute the initial Y and YPRIME
1292
C vectors. This could happen because the initial
1293
C approximation to Y or YPRIME was not very good, or
1294
C because no consistent values of these vectors exist.
1295
C The problem could also be caused by an inaccurate or
1296
C singular iteration matrix, or a poor preconditioner.
1298
C IDID = -13, there was an unrecoverable error encountered inside
1299
C your PSOL routine. Determine the cause before
1302
C IDID = -14, the Krylov linear system solver failed to achieve
1303
C convergence. This may be due to ill-conditioning
1304
C in the iteration matrix, or a singularity in the
1305
C preconditioner (if one is being used).
1306
C Another possibility is that there is a better
1307
C choice of Krylov parameters (see INFO(13)).
1308
C Possibly the failure is caused by redundant equations
1309
C in the system, or by inconsistent equations.
1310
C In that case, reformulate the system to make it
1311
C consistent and non-redundant.
1313
C IDID = -15,..,-32 --- Cannot occur with this code.
1315
C *** FOLLOWING A TERMINATED TASK ***
1317
C If IDID = -33, you cannot continue the solution of this problem.
1318
C An attempt to do so will result in your run being
1321
C ---------------------------------------------------------------------
1324
C 1. L. R. Petzold, A Description of DASSL: A Differential/Algebraic
1325
C System Solver, in Scientific Computing, R. S. Stepleman et al.
1326
C (Eds.), North-Holland, Amsterdam, 1983, pp. 65-68.
1327
C 2. K. E. Brenan, S. L. Campbell, and L. R. Petzold, Numerical
1328
C Solution of Initial-Value Problems in Differential-Algebraic
1329
C Equations, Elsevier, New York, 1989.
1330
C 3. P. N. Brown and A. C. Hindmarsh, Reduced Storage Matrix Methods
1331
C in Stiff ODE Systems, J. Applied Mathematics and Computation,
1332
C 31 (1989), pp. 40-91.
1333
C 4. P. N. Brown, A. C. Hindmarsh, and L. R. Petzold, Using Krylov
1334
C Methods in the Solution of Large-Scale Differential-Algebraic
1335
C Systems, SIAM J. Sci. Comp., 15 (1994), pp. 1467-1488.
1336
C 5. P. N. Brown, A. C. Hindmarsh, and L. R. Petzold, Consistent
1337
C Initial Condition Calculation for Differential-Algebraic
1338
C Systems, SIAM J. Sci. Comp. 19 (1998), pp. 1495-1512.
1342
C The following are all the subordinate routines used by DDASKR.
1344
C DRCHEK does preliminary checking for roots, and serves as an
1345
C interface between Subroutine DDASKR and Subroutine DROOTS.
1346
C DROOTS finds the leftmost root of a set of functions.
1347
C DDASIC computes consistent initial conditions.
1348
C DYYPNW updates Y and YPRIME in linesearch for initial condition
1350
C DDSTP carries out one step of the integration.
1351
C DCNSTR/DCNST0 check the current solution for constraint violations.
1352
C DDAWTS sets error weight quantities.
1353
C DINVWT tests and inverts the error weights.
1354
C DDATRP performs interpolation to get an output solution.
1355
C DDWNRM computes the weighted root-mean-square norm of a vector.
1356
C D1MACH provides the unit roundoff of the computer.
1357
C XERRWD/XSETF/XSETUN/IXSAV is a package to handle error messages.
1358
C DDASID nonlinear equation driver to initialize Y and YPRIME using
1359
C direct linear system solver methods. Interfaces to Newton
1360
C solver (direct case).
1361
C DNSID solves the nonlinear system for unknown initial values by
1362
C modified Newton iteration and direct linear system methods.
1363
C DLINSD carries out linesearch algorithm for initial condition
1364
C calculation (direct case).
1365
C DFNRMD calculates weighted norm of preconditioned residual in
1366
C initial condition calculation (direct case).
1367
C DNEDD nonlinear equation driver for direct linear system solver
1368
C methods. Interfaces to Newton solver (direct case).
1369
C DMATD assembles the iteration matrix (direct case).
1370
C DNSD solves the associated nonlinear system by modified
1371
C Newton iteration and direct linear system methods.
1372
C DSLVD interfaces to linear system solver (direct case).
1373
C DDASIK nonlinear equation driver to initialize Y and YPRIME using
1374
C Krylov iterative linear system methods. Interfaces to
1375
C Newton solver (Krylov case).
1376
C DNSIK solves the nonlinear system for unknown initial values by
1377
C Newton iteration and Krylov iterative linear system methods.
1378
C DLINSK carries out linesearch algorithm for initial condition
1379
C calculation (Krylov case).
1380
C DFNRMK calculates weighted norm of preconditioned residual in
1381
C initial condition calculation (Krylov case).
1382
C DNEDK nonlinear equation driver for iterative linear system solver
1383
C methods. Interfaces to Newton solver (Krylov case).
1384
C DNSK solves the associated nonlinear system by Inexact Newton
1385
C iteration and (linear) Krylov iteration.
1386
C DSLVK interfaces to linear system solver (Krylov case).
1387
C DSPIGM solves a linear system by SPIGMR algorithm.
1388
C DATV computes matrix-vector product in Krylov algorithm.
1389
C DORTH performs orthogonalization of Krylov basis vectors.
1390
C DHEQR performs QR factorization of Hessenberg matrix.
1391
C DHELS finds least-squares solution of Hessenberg linear system.
1392
C DGEFA, DGESL, DGBFA, DGBSL are LINPACK routines for solving
1393
C linear systems (dense or band direct methods).
1394
C DAXPY, DCOPY, DDOT, DNRM2, DSCAL are Basic Linear Algebra (BLAS)
1397
C The routines called directly by DDASKR are:
1398
C DCNST0, DDAWTS, DINVWT, D1MACH, DDWNRM, DDASIC, DDATRP, DDSTP,
1401
C***END PROLOGUE DDASKR
1404
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1405
LOGICAL DONE, LAVL, LCFN, LCFL, LWARN
1406
DIMENSION Y(*),YPRIME(*)
1408
DIMENSION RWORK(LRW),IWORK(LIW)
1409
DIMENSION RTOL(*),ATOL(*)
1410
DIMENSION RPAR(*),IPAR(*)
1412
EXTERNAL RES, JAC, PSOL, RT, DDASID, DDASIK, DNEDD, DNEDK
1414
C Set pointers into IWORK.
1416
PARAMETER (LML=1, LMU=2, LMTYPE=4,
1417
* LIWM=1, LMXORD=3, LJCALC=5, LPHASE=6, LK=7, LKOLD=8,
1418
* LNS=9, LNSTL=10, LNST=11, LNRE=12, LNJE=13, LETF=14, LNCFN=15,
1419
* LNCFL=16, LNIW=17, LNRW=18, LNNI=19, LNLI=20, LNPS=21,
1420
* LNPD=22, LMITER=23, LMAXL=24, LKMP=25, LNRMAX=26, LLNWP=27,
1421
* LLNIWP=28, LLOCWP=29, LLCIWP=30, LKPRIN=31, LMXNIT=32,
1422
* LMXNJ=33, LMXNH=34, LLSOFF=35, LNRTE=36, LIRFND=37, LICNS=41)
1424
C Set pointers into RWORK.
1426
PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4, LCJ=5, LCJOLD=6,
1427
* LHOLD=7, LS=8, LROUND=9, LEPLI=10, LSQRN=11, LRSQRN=12,
1428
* LEPCON=13, LSTOL=14, LEPIN=15, LALPHA=21, LBETA=27,
1429
* LGAMMA=33, LPSI=39, LSIGMA=45, LT0=51, LTLAST=52, LDELTA=61)
1431
SAVE LID, LENID, NONNEG, NCPHI
1434
C***FIRST EXECUTABLE STATEMENT DDASKR
1438
IF(INFO(1).NE.0) GO TO 100
1441
C-----------------------------------------------------------------------
1442
C This block is executed for the initial call only.
1443
C It contains checking of inputs and initializations.
1444
C-----------------------------------------------------------------------
1446
C First check INFO array to make sure all elements of INFO
1447
C Are within the proper range. (INFO(1) is checked later, because
1448
C it must be tested on every call.) ITEMP holds the location
1449
C within INFO which may be out of range.
1453
IF (INFO(I) .NE. 0 .AND. INFO(I) .NE. 1) GO TO 701
1456
IF(INFO(10).LT.0 .OR. INFO(10).GT.3) GO TO 701
1458
IF(INFO(11).LT.0 .OR. INFO(11).GT.2) GO TO 701
1461
IF (INFO(I) .NE. 0 .AND. INFO(I) .NE. 1) GO TO 701
1464
IF(INFO(18).LT.0 .OR. INFO(18).GT.2) GO TO 701
1467
C Check NEQ to see if it is positive.
1469
IF (NEQ .LE. 0) GO TO 702
1471
C Check and compute maximum order.
1474
IF (INFO(9) .NE. 0) THEN
1476
IF (MXORD .LT. 1 .OR. MXORD .GT. 5) GO TO 703
1480
C Set and/or check inputs for constraint checking (INFO(10) .NE. 0).
1481
C Set values for ICNFLG, NONNEG, and pointer LID.
1486
IF (INFO(10) .EQ. 0) GO TO 20
1487
IF (INFO(10) .EQ. 1) THEN
1491
ELSEIF (INFO(10) .EQ. 2) THEN
1502
C Set and/or check inputs for Krylov solver (INFO(12) .NE. 0).
1503
C If indicated, set default values for MAXL, KMP, NRMAX, and EPLI.
1504
C Otherwise, verify inputs required for iterative solver.
1506
IF (INFO(12) .EQ. 0) GO TO 25
1508
IWORK(LMITER) = INFO(12)
1509
IF (INFO(13) .EQ. 0) THEN
1510
IWORK(LMAXL) = MIN(5,NEQ)
1511
IWORK(LKMP) = IWORK(LMAXL)
1513
RWORK(LEPLI) = 0.05D0
1515
IF(IWORK(LMAXL) .LT. 1 .OR. IWORK(LMAXL) .GT. NEQ) GO TO 720
1516
IF(IWORK(LKMP) .LT. 1 .OR. IWORK(LKMP) .GT. IWORK(LMAXL))
1518
IF(IWORK(LNRMAX) .LT. 0) GO TO 722
1519
IF(RWORK(LEPLI).LE.0.0D0 .OR. RWORK(LEPLI).GE.1.0D0)GO TO 723
1524
C Set and/or check controls for the initial condition calculation
1525
C (INFO(11) .GT. 0). If indicated, set default values.
1526
C Otherwise, verify inputs required for iterative solver.
1528
IF (INFO(11) .EQ. 0) GO TO 30
1529
IF (INFO(17) .EQ. 0) THEN
1531
IF (INFO(12) .GT. 0) IWORK(LMXNIT) = 15
1533
IF (INFO(12) .GT. 0) IWORK(LMXNJ) = 2
1536
RWORK(LEPIN) = 0.01D0
1538
IF (IWORK(LMXNIT) .LE. 0) GO TO 725
1539
IF (IWORK(LMXNJ) .LE. 0) GO TO 725
1540
IF (IWORK(LMXNH) .LE. 0) GO TO 725
1541
LSOFF = IWORK(LLSOFF)
1542
IF (LSOFF .LT. 0 .OR. LSOFF .GT. 1) GO TO 725
1543
IF (RWORK(LEPIN) .LE. 0.0D0) GO TO 725
1548
C Below is the computation and checking of the work array lengths
1549
C LENIW and LENRW, using direct methods (INFO(12) = 0) or
1550
C the Krylov methods (INFO(12) = 1).
1553
IF (INFO(10) .EQ. 1 .OR. INFO(10) .EQ. 3) LENIC = NEQ
1555
IF (INFO(11) .EQ. 1 .OR. INFO(16) .EQ. 1) LENID = NEQ
1556
IF (INFO(12) .EQ. 0) THEN
1558
C Compute MTYPE, etc. Check ML and MU.
1560
NCPHI = MAX(MXORD + 1, 4)
1561
IF(INFO(6).EQ.0) THEN
1563
LENRW = 60 + 3*NRT + (NCPHI+3)*NEQ + LENPD
1564
IF(INFO(5).EQ.0) THEN
1570
IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717
1571
IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718
1572
LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ
1573
IF(INFO(5).EQ.0) THEN
1575
MBAND=IWORK(LML)+IWORK(LMU)+1
1577
LENRW = 60 + 3*NRT + (NCPHI+3)*NEQ + LENPD + 2*MSAVE
1580
LENRW = 60 + 3*NRT + (NCPHI+3)*NEQ + LENPD
1584
C Compute LENIW, LENWP, LENIWP.
1586
LENIW = 40 + LENIC + LENID + NEQ
1590
ELSE IF (INFO(12) .EQ. 1) THEN
1593
LENWP = IWORK(LLNWP)
1594
LENIWP = IWORK(LLNIWP)
1595
LENPD = (MAXL+3+MIN0(1,MAXL-IWORK(LKMP)))*NEQ
1596
1 + (MAXL+3)*MAXL + 1 + LENWP
1597
LENRW = 60 + 3*NRT + (MXORD+5)*NEQ + LENPD
1598
LENIW = 40 + LENIC + LENID + LENIWP
1601
IF(INFO(16) .NE. 0) LENRW = LENRW + NEQ
1603
C Check lengths of RWORK and IWORK.
1605
c -------------- memory allocation for masking ----------
1607
c -------------- masking ------------------------------
1611
IWORK(LLOCWP) = LENPD-LENWP+1
1612
IF(LRW.LT.LENRW)GO TO 704
1613
IF(LIW.LT.LENIW)GO TO 705
1615
C Check ICNSTR for legality.
1617
IF (LENIC .GT. 0) THEN
1619
ICI = IWORK(LICNS-1+I)
1620
IF (ICI .LT. -2 .OR. ICI .GT. 2) GO TO 726
1624
C Check Y for consistency with constraints.
1626
IF (LENIC .GT. 0) THEN
1627
CALL DCNST0(NEQ,Y,IWORK(LICNS),IRET)
1628
IF (IRET .NE. 0) GO TO 727
1631
C Check ID for legality and set INDEX = 0 or 1.
1634
IF (LENID .GT. 0) THEN
1637
IDI = IWORK(LID-1+I)
1638
IF (IDI .NE. 1 .AND. IDI .NE. -1) GO TO 724
1639
IF (IDI .EQ. -1) INDEX = 1
1643
C Check to see that TOUT is different from T, and NRT .ge. 0.
1645
IF(TOUT .EQ. T)GO TO 719
1646
IF(NRT .LT. 0) GO TO 730
1650
IF(INFO(7) .NE. 0) THEN
1652
IF (HMAX .LE. 0.0D0) GO TO 710
1655
C Initialize counters and other flags.
1667
IWORK(LKPRIN)=INFO(18)
1671
C-----------------------------------------------------------------------
1672
C This block is for continuation calls only.
1673
C Here we check INFO(1), and if the last step was interrupted,
1674
C we check whether appropriate action was taken.
1675
C-----------------------------------------------------------------------
1678
IF(INFO(1).EQ.1)GO TO 110
1680
IF(INFO(1).NE.-1)GO TO 701
1682
C If we are here, the last step was interrupted by an error
1683
C condition from DDSTP, and appropriate action was not taken.
1684
C This is a fatal error.
1686
MSG = 'DASKR-- THE LAST STEP TERMINATED WITH A NEGATIVE'
1687
CALL XERRWD(MSG,49,201,0,0,0,0,0,0.0D0,0.0D0)
1688
MSG = 'DASKR-- VALUE (=I1) OF IDID AND NO APPROPRIATE'
1689
CALL XERRWD(MSG,47,202,0,1,IDID,0,0,0.0D0,0.0D0)
1690
MSG = 'DASKR-- ACTION WAS TAKEN. RUN TERMINATED'
1691
CALL XERRWD(MSG,41,203,1,0,0,0,0,0.0D0,0.0D0)
1695
C-----------------------------------------------------------------------
1696
C This block is executed on all calls.
1698
C Counters are saved for later checks of performance.
1699
C Then the error tolerance parameters are checked, and the
1700
C work array pointers are set.
1701
C-----------------------------------------------------------------------
1705
C Save counters for use later.
1707
IWORK(LNSTL)=IWORK(LNST)
1710
NCFN0 = IWORK(LNCFN)
1711
NCFL0 = IWORK(LNCFL)
1714
C Check RTOL and ATOL.
1720
IF (INFO(2) .EQ. 1) RTOLI = RTOL(I)
1721
IF (INFO(2) .EQ. 1) ATOLI = ATOL(I)
1722
IF (RTOLI .GT. 0.0D0 .OR. ATOLI .GT. 0.0D0) NZFLG = 1
1723
IF (RTOLI .LT. 0.0D0) GO TO 706
1724
IF (ATOLI .LT. 0.0D0) GO TO 707
1726
IF (NZFLG .EQ. 0) GO TO 708
1728
C Set pointers to RWORK and IWORK segments.
1729
C For direct methods, SAVR is not used.
1731
IWORK(LLCIWP) = LID + LENID
1733
IF (INFO(12) .NE. 0) LSAVR = LDELTA + NEQ
1737
IF (INFO(16) .NE. 0) LVT = LWT + NEQ
1739
LR0 = LPHI + NCPHI*NEQ
1743
IF (INFO(1) .EQ. 1) GO TO 400
1745
C-----------------------------------------------------------------------
1746
C This block is executed on the initial call only.
1747
C Set the initial step size, the error weight vector, and PHI.
1748
C Compute unknown initial components of Y and YPRIME, if requested.
1749
C-----------------------------------------------------------------------
1755
C Set error weight array WT and altered weight array VT.
1757
CALL DDAWTS1(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
1758
CALL DINVWT(NEQ,RWORK(LWT),IER)
1759
IF (IER .NE. 0) GO TO 713
1760
IF (INFO(16) .NE. 0) THEN
1762
305 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1)
1765
C Compute unit roundoff and HMIN. >>> instead of D1MACH(4) we use
1766
c DLAMCH, because the optimized compiler affects the D1MACH.
1767
c UROUND = D1MACH(4)
1768
UROUND = DLAMCH('p')
1769
RWORK(LROUND) = UROUND
1770
c ---------------- Hmind chnage ---------------------
1771
c HMIN = 4.0D0*UROUND*MAX(ABS(T),ABS(TOUT))
1773
c ---------------- Hmind chnage ---------------------
1775
C Set/check STPTOL control for initial condition calculation.
1777
IF (INFO(11) .NE. 0) THEN
1778
IF( INFO(17) .EQ. 0) THEN
1779
RWORK(LSTOL) = UROUND**.6667D0
1781
IF (RWORK(LSTOL) .LE. 0.0D0) GO TO 725
1785
C Compute EPCON and square root of NEQ and its reciprocal, used
1786
C inside iterative solver.
1788
RWORK(LEPCON) = 0.33D0
1790
RWORK(LSQRN) = SQRT(FLOATN)
1791
RWORK(LRSQRN) = 1.D0/RWORK(LSQRN)
1793
C Check initial interval to see that it is long enough.
1795
TDIST = ABS(TOUT - T)
1796
c ---------------- Hmind chnage ---------------------
1797
cc IF(TDIST .LT. HMIN) GO TO 714
1798
c ---------------- Hmind chnage ---------------------
1800
C Check H0, if this was input.
1802
IF (INFO(8) .EQ. 0) GO TO 310
1804
IF ((TOUT - T)*H0 .LT. 0.0D0) GO TO 711
1805
IF (H0 .EQ. 0.0D0) GO TO 712
1809
C Compute initial stepsize, to be used by either
1810
C DDSTP or DDASIC, depending on INFO(11).
1813
YPNORM = DDWNRM(NEQ,YPRIME,RWORK(LVT),RPAR,IPAR)
1814
IF (YPNORM .GT. 0.5D0/H0) H0 = 0.5D0/YPNORM
1815
H0 = SIGN(H0,TOUT-T)
1817
C Adjust H0 if necessary to meet HMAX bound.
1819
320 IF (INFO(7) .EQ. 0) GO TO 330
1820
RH = ABS(H0)/RWORK(LHMAX)
1821
IF (RH .GT. 1.0D0) H0 = H0/RH
1823
C Check against TSTOP, if applicable.
1825
330 IF (INFO(4) .EQ. 0) GO TO 340
1826
TSTOP = RWORK(LTSTOP)
1827
IF ((TSTOP - T)*H0 .LT. 0.0D0) GO TO 715
1828
IF ((T + H0 - TSTOP)*H0 .GT. 0.0D0) H0 = TSTOP - T
1829
IF ((TSTOP - TOUT)*H0 .LT. 0.0D0) GO TO 709
1831
340 IF (INFO(11) .EQ. 0) GO TO 370
1833
C Compute unknown components of initial Y and YPRIME, depending
1834
C on INFO(11) and INFO(12). INFO(12) represents the nonlinear
1835
C solver type (direct/Krylov). Pass the name of the specific
1836
C nonlinear solver, depending on INFO(12). The location of the work
1837
C arrays SAVR, YIC, YPIC, PWK also differ in the two cases.
1838
C For use in stopping tests, pass TSCALE = TDIST if INDEX = 0.
1843
EPCONI = RWORK(LEPIN)*RWORK(LEPCON)
1845
IF (INDEX .EQ. 0) TSCALE = TDIST
1846
350 IF (INFO(12) .EQ. 0) THEN
1850
CALL DDASIC(TN,Y,YPRIME,NEQ,INFO(11),IWORK(LID),
1851
* RES,JAC,PSOL,H0,TSCALE,RWORK(LWT),NWT,IDID,RPAR,IPAR,
1852
* RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE),
1853
* RWORK(LYIC),RWORK(LYPIC),RWORK(LPWK),RWORK(LWM),IWORK(LIWM),
1854
* RWORK(LROUND),RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN),
1855
* EPCONI,RWORK(LSTOL),INFO(15),ICNFLG,IWORK(LICNS),DDASID)
1856
ELSE IF (INFO(12) .EQ. 1) THEN
1860
CALL DDASIC(TN,Y,YPRIME,NEQ,INFO(11),IWORK(LID),
1861
* RES,JAC,PSOL,H0,TSCALE,RWORK(LWT),NWT,IDID,RPAR,IPAR,
1862
* RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE),
1863
* RWORK(LYIC),RWORK(LYPIC),RWORK(LPWK),RWORK(LWM),IWORK(LIWM),
1864
* RWORK(LROUND),RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN),
1865
* EPCONI,RWORK(LSTOL),INFO(15),ICNFLG,IWORK(LICNS),DDASIK)
1868
IF (IDID .LT. 0) GO TO 600
1870
C DDASIC was successful. If this was the first call to DDASIC,
1871
C update the WT array (with the current Y) and call it again.
1873
IF (NWT .EQ. 2) GO TO 355
1875
CALL DDAWTS1(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
1876
CALL DINVWT(NEQ,RWORK(LWT),IER)
1877
IF (IER .NE. 0) GO TO 713
1880
C If INFO(14) = 1, return now with IDID = 4.
1883
355 IF (INFO(14) .EQ. 1) THEN
1886
IF (INFO(11) .EQ. 1) RWORK(LHOLD) = H0
1890
C Update the WT and VT arrays one more time, with the new Y.
1892
CALL DDAWTS1(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
1893
CALL DINVWT(NEQ,RWORK(LWT),IER)
1894
IF (IER .NE. 0) GO TO 713
1895
IF (INFO(16) .NE. 0) THEN
1897
357 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1)
1900
C Reset the initial stepsize to be used by DDSTP.
1901
C Use H0, if this was input. Otherwise, recompute H0,
1902
C and adjust it if necessary to meet HMAX bound.
1904
IF (INFO(8) .NE. 0) THEN
1910
YPNORM = DDWNRM(NEQ,YPRIME,RWORK(LVT),RPAR,IPAR)
1911
IF (YPNORM .GT. 0.5D0/H0) H0 = 0.5D0/YPNORM
1912
H0 = SIGN(H0,TOUT-T)
1914
360 IF (INFO(7) .NE. 0) THEN
1915
RH = ABS(H0)/RWORK(LHMAX)
1916
IF (RH .GT. 1.0D0) H0 = H0/RH
1919
C Check against TSTOP, if applicable.
1921
IF (INFO(4) .NE. 0) THEN
1922
TSTOP = RWORK(LTSTOP)
1923
IF ((T + H0 - TSTOP)*H0 .GT. 0.0D0) H0 = TSTOP - T
1926
C Load H and RWORK(LH) with H0.
1931
C Load Y and H*YPRIME into PHI(*,1) and PHI(*,2).
1935
RWORK(LPHI + I - 1) = Y(I)
1936
380 RWORK(ITEMP + I - 1) = H*YPRIME(I)
1938
C Initialize T0 in RWORK; check for a zero of R near initial T.
1943
RWORK(LPSI+1)=2.0D0*H
1945
IF (NRT .EQ. 0) GO TO 390
1946
c -------------- masking ----------------->>>
1947
CALL DRCHEK1(1,RT,NRT,NEQ,T,TOUT,Y,YPRIME,RWORK(LPHI),
1948
* RWORK(LPSI),IWORK(LKOLD),RWORK(LR0),RWORK(LR1),
1949
* RWORK(LRX),JROOT,IRT,RWORK(LROUND),INFO(3),
1950
* RWORK,IWORK,RPAR,IPAR)
1951
IF (IRT .LT. 0) GO TO 731
1952
IF (IRT .EQ. 2) IWORK(LIRFND) = 2
1956
C-----------------------------------------------------------------------
1957
C This block is for continuation calls only.
1958
C Its purpose is to check stop conditions before taking a step.
1959
C Adjust H if necessary to meet HMAX bound.
1960
C-----------------------------------------------------------------------
1963
UROUND=RWORK(LROUND)
1967
IF(NRT .EQ. 0) GO TO 405
1969
C Check for a zero of R near TN.
1971
CALL DRCHEK1(2,RT,NRT,NEQ,TN,TOUT,Y,YPRIME,RWORK(LPHI),
1972
* RWORK(LPSI),IWORK(LKOLD),RWORK(LR0),RWORK(LR1),
1973
* RWORK(LRX),JROOT,IRT,RWORK(LROUND),INFO(3),
1974
* RWORK,IWORK,RPAR,IPAR)
1975
IF (IRT .LT. 0) GO TO 731
1977
IF (IRT .EQ. 1) THEN
1987
IF(INFO(7) .EQ. 0) GO TO 410
1988
RH = ABS(H)/RWORK(LHMAX)
1989
IF(RH .GT. 1.0D0) H = H/RH
1991
IF(T .EQ. TOUT) GO TO 719
1992
IF((T - TOUT)*H .GT. 0.0D0) GO TO 711
1993
IF(INFO(4) .EQ. 1) GO TO 430
1994
IF(INFO(3) .EQ. 1) GO TO 420
1995
IF((TN-TOUT)*H.LT.0.0D0)GO TO 490
1996
CALL DDATRP1(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1997
* RWORK(LPHI),RWORK(LPSI))
2002
420 IF((TN-T)*H .LE. 0.0D0) GO TO 490
2003
IF((TN - TOUT)*H .GE. 0.0D0) GO TO 425
2004
CALL DDATRP1(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
2005
* RWORK(LPHI),RWORK(LPSI))
2011
CALL DDATRP1(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
2012
* RWORK(LPHI),RWORK(LPSI))
2017
430 IF(INFO(3) .EQ. 1) GO TO 440
2019
IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715
2020
IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709
2021
IF((TN-TOUT)*H.LT.0.0D0)GO TO 450
2022
CALL DDATRP1(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
2023
* RWORK(LPHI),RWORK(LPSI))
2028
440 TSTOP = RWORK(LTSTOP)
2029
IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715
2030
IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709
2031
IF((TN-T)*H .LE. 0.0D0) GO TO 450
2032
IF((TN - TOUT)*H .GE. 0.0D0) GO TO 445
2033
CALL DDATRP1(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
2034
* RWORK(LPHI),RWORK(LPSI))
2040
CALL DDATRP1(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
2041
* RWORK(LPHI),RWORK(LPSI))
2048
C Check whether we are within roundoff of TSTOP.
2050
IF(ABS(TN-TSTOP).GT.100.0D0*UROUND*
2051
* (ABS(TN)+ABS(H)))GO TO 460
2052
CALL DDATRP1(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD),
2053
* RWORK(LPHI),RWORK(LPSI))
2059
IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490
2063
490 IF (DONE) GO TO 590
2065
C-----------------------------------------------------------------------
2066
C The next block contains the call to the one-step integrator DDSTP.
2067
C This is a looping point for the integration steps.
2068
C Check for too many steps.
2069
C Check for poor Newton/Krylov performance.
2070
C Update WT. Check for too much accuracy requested.
2071
C Compute minimum stepsize.
2072
C-----------------------------------------------------------------------
2076
C Check for too many steps.
2078
IF((IWORK(LNST)-IWORK(LNSTL)).LT.500) GO TO 505
2082
C Check for poor Newton/Krylov performance.
2084
505 IF (INFO(12) .EQ. 0) GO TO 510
2085
NSTD = IWORK(LNST) - IWORK(LNSTL)
2086
NNID = IWORK(LNNI) - NNI0
2087
IF (NSTD .LT. 10 .OR. NNID .EQ. 0) GO TO 510
2088
AVLIN = REAL(IWORK(LNLI) - NLI0)/REAL(NNID)
2089
RCFN = REAL(IWORK(LNCFN) - NCFN0)/REAL(NSTD)
2090
RCFL = REAL(IWORK(LNCFL) - NCFL0)/REAL(NNID)
2091
FMAXL = IWORK(LMAXL)
2092
LAVL = AVLIN .GT. FMAXL
2093
LCFN = RCFN .GT. 0.9D0
2094
LCFL = RCFL .GT. 0.9D0
2095
LWARN = LAVL .OR. LCFN .OR. LCFL
2096
IF (.NOT.LWARN) GO TO 510
2098
IF (NWARN .GT. 10) GO TO 510
2100
MSG = 'DASKR-- Warning. Poor iterative algorithm performance '
2101
CALL XERRWD (MSG, 56, 501, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
2102
MSG = ' at T = R1. Average no. of linear iterations = R2 '
2103
CALL XERRWD (MSG, 56, 501, 0, 0, 0, 0, 2, TN, AVLIN)
2106
MSG = 'DASKR-- Warning. Poor iterative algorithm performance '
2107
CALL XERRWD (MSG, 56, 502, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
2108
MSG = ' at T = R1. Nonlinear convergence failure rate = R2'
2109
CALL XERRWD (MSG, 56, 502, 0, 0, 0, 0, 2, TN, RCFN)
2112
MSG = 'DASKR-- Warning. Poor iterative algorithm performance '
2113
CALL XERRWD (MSG, 56, 503, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
2114
MSG = ' at T = R1. Linear convergence failure rate = R2 '
2115
CALL XERRWD (MSG, 56, 503, 0, 0, 0, 0, 2, TN, RCFL)
2118
C Update WT and VT, if this is not the first call.
2120
510 CALL DDAWTS1(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI),RWORK(LWT),
2122
CALL DINVWT(NEQ,RWORK(LWT),IER)
2123
IF (IER .NE. 0) THEN
2127
IF (INFO(16) .NE. 0) THEN
2129
RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1)
2133
C Test for too much accuracy requested.
2135
R = DDWNRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)*100.0D0*UROUND
2136
IF (R .LE. 1.0D0) GO TO 525
2138
C Multiply RTOL and ATOL by R and return.
2140
IF(INFO(2).EQ.1)GO TO 523
2153
C Compute minimum stepsize.
2155
c ----------------------- Hmin Change----------------------
2157
c HMIN is intended to be a value slightly above the roundoff level
2158
c in the current T. As such it is appropriate that it varies with
2159
c T. In DASKR, HMIN is used in two ways:
2161
c 1. At the start, ABS(TOUT - T) is required to be at least HMIN, to
2162
c guarantee that the user has provided the direction of the
2163
c integration reliably. For this test, it would be sufficient to
2164
c ignore HMIN and simply require that TOUT - T is nonzero on the
2167
c 2. If the integration has difficulty passing the convergence or
2168
c the error test with step size H of size ABS(H) = HMIN, it stops
2169
c with an error message saying that. In all of these uses, it would
2170
c not hurt to use HMIN = 0, in my opinion. There are cases where
2171
c the appropriate step size is temporarily below the roundoff level
2172
c in T. The only negative impact of using HMIN = 0 then is that
2173
c some steps may be taken in which T + H = T on the machine. In
2174
c contrast to the DASSL family, the ODEPACK solvers have HMIN = 0 as
2175
c the default in these uses of HMIN, but they issue a warning when T
2176
c + H = T, because in most cases this is the result of a user
2177
c program bug or input error of some kind. On the other hand, the
2178
c value HMIN = 4*UROUND makes no sense, because it ignores the scale
2179
c of the T variable completely. Thus it could ause invalid error
2180
c halts, when values of ABS(H) smaller than that may well be
2181
c appropriate. If the current HMIN is bothersome, I suggest using
2182
c HMIN = 0, and removing HMIN from the initial test on TOUT - T
2185
c HMIN=4.0D0*UROUND*MAX(ABS(TN),ABS(TOUT))
2187
c ----------------------- Hmin Change----------------------
2190
IF (INFO(7) .NE. 0) THEN
2191
RH = ABS(H)/RWORK(LHMAX)
2192
IF (RH .GT. 1.0D0) H = H/RH
2195
C Call the one-step integrator.
2196
C Note that INFO(12) represents the nonlinear solver type.
2197
C Pass the required nonlinear solver, depending upon INFO(12).
2199
c info(12): 0-> dierct case 1->Krylov
2200
IF (INFO(12) .EQ. 0) THEN
2201
CALL DDSTP(TN,Y,YPRIME,NEQ,
2202
* RES,JAC,PSOL,H,RWORK(LWT),RWORK(LVT),INFO(1),IDID,RPAR,IPAR,
2203
* RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE),
2204
* RWORK(LWM),IWORK(LIWM),
2205
* RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA),
2206
* RWORK(LPSI),RWORK(LSIGMA),
2207
* RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),RWORK(LS),HMIN,
2208
* RWORK(LROUND), RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN),
2209
* RWORK(LEPCON), IWORK(LPHASE),IWORK(LJCALC),INFO(15),
2210
* IWORK(LK), IWORK(LKOLD),IWORK(LNS),NONNEG,INFO(12),
2212
ELSE IF (INFO(12) .EQ. 1) THEN
2213
CALL DDSTP(TN,Y,YPRIME,NEQ,
2214
* RES,JAC,PSOL,H,RWORK(LWT),RWORK(LVT),INFO(1),IDID,RPAR,IPAR,
2215
* RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE),
2216
* RWORK(LWM),IWORK(LIWM),
2217
* RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA),
2218
* RWORK(LPSI),RWORK(LSIGMA),
2219
* RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),RWORK(LS),HMIN,
2220
* RWORK(LROUND), RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN),
2221
* RWORK(LEPCON), IWORK(LPHASE),IWORK(LJCALC),INFO(15),
2222
* IWORK(LK), IWORK(LKOLD),IWORK(LNS),NONNEG,INFO(12),
2226
527 IF(IDID.LT.0)GO TO 600
2228
C-----------------------------------------------------------------------
2229
C This block handles the case of a successful return from DDSTP
2230
C (IDID=1). Test for stop conditions.
2231
C-----------------------------------------------------------------------
2233
IF(NRT .EQ. 0) GO TO 529
2235
C Check for a zero of R near TN.
2237
CALL DRCHEK1(3,RT,NRT,NEQ,TN,TOUT,Y,YPRIME,RWORK(LPHI),
2238
* RWORK(LPSI),IWORK(LKOLD),RWORK(LR0),RWORK(LR1),
2239
* RWORK(LRX),JROOT,IRT,RWORK(LROUND),INFO(3),
2240
* RWORK,IWORK,RPAR,IPAR)
2249
IF(IRT .NE. 1) GO TO 529
2256
IF(INFO(4).NE.0)GO TO 540
2257
IF(INFO(3).NE.0)GO TO 530
2258
IF((TN-TOUT)*H.LT.0.0D0)GO TO 500
2259
CALL DDATRP1(TN,TOUT,Y,YPRIME,NEQ,
2260
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
2264
530 IF((TN-TOUT)*H.GE.0.0D0)GO TO 535
2268
535 CALL DDATRP1(TN,TOUT,Y,YPRIME,NEQ,
2269
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
2273
540 IF(INFO(3).NE.0)GO TO 550
2274
IF((TN-TOUT)*H.LT.0.0D0)GO TO 542
2275
CALL DDATRP1(TN,TOUT,Y,YPRIME,NEQ,
2276
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
2280
542 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*
2281
* (ABS(TN)+ABS(H)))GO TO 545
2283
IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500
2286
545 CALL DDATRP1(TN,TSTOP,Y,YPRIME,NEQ,
2287
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
2291
550 IF((TN-TOUT)*H.GE.0.0D0)GO TO 555
2292
IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*(ABS(TN)+ABS(H)))GO TO 552
2296
552 CALL DDATRP1(TN,TSTOP,Y,YPRIME,NEQ,
2297
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
2301
555 CALL DDATRP1(TN,TOUT,Y,YPRIME,NEQ,
2302
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
2307
C-----------------------------------------------------------------------
2308
C All successful returns from DDASKR are made from this block.
2309
C-----------------------------------------------------------------------
2318
C-----------------------------------------------------------------------
2319
C This block handles all unsuccessful returns other than for
2321
C-----------------------------------------------------------------------
2325
GO TO (610,620,630,700,655,640,650,660,670,675,
2326
* 680,685,690,695,696), ITEMP
2328
C The maximum number of steps was taken before
2331
610 MSG = 'DASKR-- AT CURRENT T (=R1) 500 STEPS'
2332
CALL XERRWD(MSG,38,610,0,0,0,0,1,TN,0.0D0)
2333
MSG = 'DASKR-- TAKEN ON THIS CALL BEFORE REACHING TOUT'
2334
CALL XERRWD(MSG,48,611,0,0,0,0,0,0.0D0,0.0D0)
2337
C Too much accuracy for machine precision.
2339
620 MSG = 'DASKR-- AT T (=R1) TOO MUCH ACCURACY REQUESTED'
2340
CALL XERRWD(MSG,47,620,0,0,0,0,1,TN,0.0D0)
2341
MSG = 'DASKR-- FOR PRECISION OF MACHINE. RTOL AND ATOL'
2342
CALL XERRWD(MSG,48,621,0,0,0,0,0,0.0D0,0.0D0)
2343
MSG = 'DASKR-- WERE INCREASED BY A FACTOR R (=R1)'
2344
CALL XERRWD(MSG,43,622,0,0,0,0,1,R,0.0D0)
2347
C WT(I) .LE. 0.0D0 for some I (not at start of problem).
2349
630 MSG = 'DASKR-- AT T (=R1) SOME ELEMENT OF WT'
2350
CALL XERRWD(MSG,38,630,0,0,0,0,1,TN,0.0D0)
2351
MSG = 'DASKR-- HAS BECOME .LE. 0.0'
2352
CALL XERRWD(MSG,28,631,0,0,0,0,0,0.0D0,0.0D0)
2355
C Error test failed repeatedly or with H=HMIN.
2357
640 MSG = 'DASKR-- AT T (=R1) AND STEPSIZE H (=R2) THE'
2358
CALL XERRWD(MSG,44,640,0,0,0,0,2,TN,H)
2359
MSG='DASKR-- ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN'
2360
CALL XERRWD(MSG,57,641,0,0,0,0,0,0.0D0,0.0D0)
2363
C Nonlinear solver failed to converge repeatedly or with H=HMIN.
2365
650 MSG = 'DASKR-- AT T (=R1) AND STEPSIZE H (=R2) THE'
2366
CALL XERRWD(MSG,44,650,0,0,0,0,2,TN,H)
2367
MSG = 'DASKR-- NONLINEAR SOLVER FAILED TO CONVERGE'
2368
CALL XERRWD(MSG,44,651,0,0,0,0,0,0.0D0,0.0D0)
2369
MSG = 'DASKR-- REPEATEDLY OR WITH ABS(H)=HMIN'
2370
CALL XERRWD(MSG,40,652,0,0,0,0,0,0.0D0,0.0D0)
2373
C The preconditioner had repeated failures.
2375
655 MSG = 'DASKR-- AT T (=R1) AND STEPSIZE H (=R2) THE'
2376
CALL XERRWD(MSG,44,655,0,0,0,0,2,TN,H)
2377
MSG = 'DASKR-- PRECONDITIONER HAD REPEATED FAILURES.'
2378
CALL XERRWD(MSG,46,656,0,0,0,0,0,0.0D0,0.0D0)
2381
C The iteration matrix is singular.
2383
660 MSG = 'DASKR-- AT T (=R1) AND STEPSIZE H (=R2) THE'
2384
CALL XERRWD(MSG,44,660,0,0,0,0,2,TN,H)
2385
MSG = 'DASKR-- ITERATION MATRIX IS SINGULAR.'
2386
CALL XERRWD(MSG,38,661,0,0,0,0,0,0.0D0,0.0D0)
2389
C Nonlinear system failure preceded by error test failures.
2391
670 MSG = 'DASKR-- AT T (=R1) AND STEPSIZE H (=R2) THE'
2392
CALL XERRWD(MSG,44,670,0,0,0,0,2,TN,H)
2393
MSG = 'DASKR-- NONLINEAR SOLVER COULD NOT CONVERGE.'
2394
CALL XERRWD(MSG,45,671,0,0,0,0,0,0.0D0,0.0D0)
2395
MSG = 'DASKR-- ALSO, THE ERROR TEST FAILED REPEATEDLY.'
2396
CALL XERRWD(MSG,49,672,0,0,0,0,0,0.0D0,0.0D0)
2399
C Nonlinear system failure because IRES = -1.
2401
675 MSG = 'DASKR-- AT T (=R1) AND STEPSIZE H (=R2) THE'
2402
CALL XERRWD(MSG,44,675,0,0,0,0,2,TN,H)
2403
MSG = 'DASKR-- NONLINEAR SYSTEM SOLVER COULD NOT CONVERGE'
2404
CALL XERRWD(MSG,51,676,0,0,0,0,0,0.0D0,0.0D0)
2405
MSG = 'DASKR-- BECAUSE IRES WAS EQUAL TO MINUS ONE'
2406
CALL XERRWD(MSG,44,677,0,0,0,0,0,0.0D0,0.0D0)
2409
C Failure because IRES = -2.
2411
680 MSG = 'DASKR-- AT T (=R1) AND STEPSIZE H (=R2)'
2412
CALL XERRWD(MSG,40,680,0,0,0,0,2,TN,H)
2413
MSG = 'DASKR-- IRES WAS EQUAL TO MINUS TWO'
2414
CALL XERRWD(MSG,36,681,0,0,0,0,0,0.0D0,0.0D0)
2417
C Failed to compute initial YPRIME.
2419
685 MSG = 'DASKR-- AT T (=R1) AND STEPSIZE H (=R2) THE'
2420
CALL XERRWD(MSG,44,685,0,0,0,0,0,0.0D0,0.0D0)
2421
MSG = 'DASKR-- INITIAL (Y,YPRIME) COULD NOT BE COMPUTED'
2422
CALL XERRWD(MSG,49,686,0,0,0,0,2,TN,H0)
2425
C Failure because IER was negative from PSOL.
2427
690 MSG = 'DASKR-- AT T (=R1) AND STEPSIZE H (=R2)'
2428
CALL XERRWD(MSG,40,690,0,0,0,0,2,TN,H)
2429
MSG = 'DASKR-- IER WAS NEGATIVE FROM PSOL'
2430
CALL XERRWD(MSG,35,691,0,0,0,0,0,0.0D0,0.0D0)
2433
C Failure because the linear system solver could not converge.
2435
695 MSG = 'DASKR-- AT T (=R1) AND STEPSIZE H (=R2) THE'
2436
CALL XERRWD(MSG,44,695,0,0,0,0,2,TN,H)
2437
MSG = 'DASKR-- LINEAR SYSTEM SOLVER COULD NOT CONVERGE.'
2438
CALL XERRWD(MSG,50,696,0,0,0,0,0,0.0D0,0.0D0)
2441
696 MSG = 'DASKR-- AT T (=R1) AND STEPSIZE H (=R2) THE'
2442
CALL XERRWD(MSG,44,685,0,0,0,0,0,0.0D0,0.0D0)
2443
MSG = 'DASKR-- INITIAL Jacobian COULD NOT BE COMPUTED'
2444
CALL XERRWD(MSG,49,686,0,0,0,0,2,TN,H0)
2455
C-----------------------------------------------------------------------
2456
C This block handles all error returns due to illegal input,
2457
C as detected before calling DDSTP.
2458
C First the error message routine is called. If this happens
2459
C twice in succession, execution is terminated.
2460
C-----------------------------------------------------------------------
2462
701 MSG = 'DASKR-- ELEMENT (=I1) OF INFO VECTOR IS NOT VALID'
2463
CALL XERRWD(MSG,50,1,0,1,ITEMP,0,0,0.0D0,0.0D0)
2465
702 MSG = 'DASKR-- NEQ (=I1) .LE. 0'
2466
CALL XERRWD(MSG,25,2,0,1,NEQ,0,0,0.0D0,0.0D0)
2468
703 MSG = 'DASKR-- MAXORD (=I1) NOT IN RANGE'
2469
CALL XERRWD(MSG,34,3,0,1,MXORD,0,0,0.0D0,0.0D0)
2471
704 MSG='DASKR-- RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)'
2472
CALL XERRWD(MSG,60,4,0,2,LENRW,LRW,0,0.0D0,0.0D0)
2474
705 MSG='DASKR-- IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)'
2475
CALL XERRWD(MSG,60,5,0,2,LENIW,LIW,0,0.0D0,0.0D0)
2477
706 MSG = 'DASKR-- SOME ELEMENT OF RTOL IS .LT. 0'
2478
CALL XERRWD(MSG,39,6,0,0,0,0,0,0.0D0,0.0D0)
2480
707 MSG = 'DASKR-- SOME ELEMENT OF ATOL IS .LT. 0'
2481
CALL XERRWD(MSG,39,7,0,0,0,0,0,0.0D0,0.0D0)
2483
708 MSG = 'DASKR-- ALL ELEMENTS OF RTOL AND ATOL ARE ZERO'
2484
CALL XERRWD(MSG,47,8,0,0,0,0,0,0.0D0,0.0D0)
2486
709 MSG='DASKR-- INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2)'
2487
CALL XERRWD(MSG,54,9,0,0,0,0,2,TSTOP,TOUT)
2489
710 MSG = 'DASKR-- HMAX (=R1) .LT. 0.0'
2490
CALL XERRWD(MSG,28,10,0,0,0,0,1,HMAX,0.0D0)
2492
711 MSG = 'DASKR-- TOUT (=R1) BEHIND T (=R2)'
2493
CALL XERRWD(MSG,34,11,0,0,0,0,2,TOUT,T)
2495
712 MSG = 'DASKR-- INFO(8)=1 AND H0=0.0'
2496
CALL XERRWD(MSG,29,12,0,0,0,0,0,0.0D0,0.0D0)
2498
713 MSG = 'DASKR-- SOME ELEMENT OF WT IS .LE. 0.0'
2499
CALL XERRWD(MSG,39,13,0,0,0,0,0,0.0D0,0.0D0)
2501
714 MSG='DASKR-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION'
2502
CALL XERRWD(MSG,60,14,0,0,0,0,2,TOUT,T)
2504
715 MSG = 'DASKR-- INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2)'
2505
CALL XERRWD(MSG,49,15,0,0,0,0,2,TSTOP,T)
2507
717 MSG = 'DASKR-- ML (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ'
2508
CALL XERRWD(MSG,52,17,0,1,IWORK(LML),0,0,0.0D0,0.0D0)
2510
718 MSG = 'DASKR-- MU (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ'
2511
CALL XERRWD(MSG,52,18,0,1,IWORK(LMU),0,0,0.0D0,0.0D0)
2513
719 MSG = 'DASKR-- TOUT (=R1) IS EQUAL TO T (=R2)'
2514
CALL XERRWD(MSG,39,19,0,0,0,0,2,TOUT,T)
2516
720 MSG = 'DASKR-- MAXL (=I1) ILLEGAL. EITHER .LT. 1 OR .GT. NEQ'
2517
CALL XERRWD(MSG,54,20,0,1,IWORK(LMAXL),0,0,0.0D0,0.0D0)
2519
721 MSG = 'DASKR-- KMP (=I1) ILLEGAL. EITHER .LT. 1 OR .GT. MAXL'
2520
CALL XERRWD(MSG,54,21,0,1,IWORK(LKMP),0,0,0.0D0,0.0D0)
2522
722 MSG = 'DASKR-- NRMAX (=I1) ILLEGAL. .LT. 0'
2523
CALL XERRWD(MSG,36,22,0,1,IWORK(LNRMAX),0,0,0.0D0,0.0D0)
2525
723 MSG = 'DASKR-- EPLI (=R1) ILLEGAL. EITHER .LE. 0.D0 OR .GE. 1.D0'
2526
CALL XERRWD(MSG,58,23,0,0,0,0,1,RWORK(LEPLI),0.0D0)
2528
724 MSG = 'DASKR-- ILLEGAL IWORK VALUE FOR INFO(11) .NE. 0'
2529
CALL XERRWD(MSG,48,24,0,0,0,0,0,0.0D0,0.0D0)
2531
725 MSG = 'DASKR-- ONE OF THE INPUTS FOR INFO(17) = 1 IS ILLEGAL'
2532
CALL XERRWD(MSG,54,25,0,0,0,0,0,0.0D0,0.0D0)
2534
726 MSG = 'DASKR-- ILLEGAL IWORK VALUE FOR INFO(10) .NE. 0'
2535
CALL XERRWD(MSG,48,26,0,0,0,0,0,0.0D0,0.0D0)
2537
727 MSG = 'DASKR-- Y(I) AND IWORK(40+I) (I=I1) INCONSISTENT'
2538
CALL XERRWD(MSG,49,27,0,1,IRET,0,0,0.0D0,0.0D0)
2540
730 MSG = 'DASKR-- NRT (=I1) .LT. 0'
2541
CALL XERRWD(MSG,25,30,1,1,NRT,0,0,0.0D0,0.0D0)
2543
731 MSG = 'DASKR-- R IS ILL-DEFINED. ZERO VALUES WERE FOUND AT TWO'
2544
CALL XERRWD(MSG,57,31,1,0,0,0,0,0.0D0,0.0D0)
2545
MSG = ' VERY CLOSE T VALUES, AT T = R1'
2546
CALL XERRWD(MSG,39,31,1,0,0,0,1,RWORK(LT0),0.0D0)
2548
750 IF(INFO(1).EQ.-1) GO TO 760
2552
760 MSG = 'DASKR-- REPEATED OCCURRENCES OF ILLEGAL INPUT'
2553
CALL XERRWD(MSG,46,701,0,0,0,0,0,0.0D0,0.0D0)
2554
770 MSG = 'DASKR-- RUN TERMINATED. APPARENT INFINITE LOOP'
2555
CALL XERRWD(MSG,47,702,1,0,0,0,0,0.0D0,0.0D0)
2558
C------END OF SUBROUTINE DDASKR-----------------------------------------
2560
SUBROUTINE DRCHEK1 (JOB, RT, NRT, NEQ, TN, TOUT, Y, YP, PHI, PSI,
2561
* KOLD, R0, R1, RX, JROOT, IRT, UROUND, INFO3, RWORK, IWORK,
2564
C*** BEGIN PROLOGUE DRCHEK
2565
C*** REFER TO DDASKR
2566
C*** ROUTINES CALLED DDATRP, DROOTS, DCOPY, RT
2567
C*** REVISION HISTORY (YYMMDD)
2568
C 020815 DATE WRITTEN
2569
C 021217 Added test for roots close when JOB = 2.
2570
C*** END PROLOGUE DRCHEK
2572
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2573
C Pointers into IWORK:
2574
PARAMETER (LNRTE=36, LIRFND=37)
2575
C Pointers into RWORK:
2576
PARAMETER (LT0=51, LTLAST=52)
2578
INTEGER JOB, NRT, NEQ, KOLD, JROOT, IRT, INFO3, IWORK, IPAR
2579
DOUBLE PRECISION TN, TOUT, Y, YP, PHI, PSI, R0, R1, RX, UROUND,
2581
DIMENSION Y(*), YP(*), PHI(NEQ,*), PSI(*),
2582
* R0(*), R1(*), RX(*), JROOT(*), RWORK(*), IWORK(*)
2583
INTEGER I, JFLAG, LMASK
2585
DOUBLE PRECISION HMINR, T1, TEMP1, TEMP2, X, ZERO
2587
c -------------- masking -----------------
2591
c -------------- masking -----------------
2593
C-----------------------------------------------------------------------
2594
C This routine checks for the presence of a root of R(T,Y,Y') in the
2595
C vicinity of the current T, in a manner depending on the
2596
C input flag JOB. It calls subroutine DROOTS to locate the root
2597
C as precisely as possible.
2599
C In addition to variables described previously, DRCHEK
2600
C uses the following for communication..
2601
C JOB = integer flag indicating type of call..
2602
C JOB = 1 means the problem is being initialized, and DRCHEK
2603
C is to look for a root at or very near the initial T.
2604
C JOB = 2 means a continuation call to the solver was just
2605
C made, and DRCHEK is to check for a root in the
2606
C relevant part of the step last taken.
2607
C JOB = 3 means a successful step was just taken, and DRCHEK
2608
C is to look for a root in the interval of the step.
2609
C R0 = array of length NRT, containing the value of R at T = T0.
2610
C R0 is input for JOB .ge. 2 and on output in all cases.
2611
C R1,RX = arrays of length NRT for work space.
2612
C IRT = completion flag..
2613
C IRT = 0 means no root was found.
2614
C IRT = -1 means JOB = 1 and a zero was found both at T0 and
2615
C and very close to T0.
2616
C IRT = -2 means JOB = 2 and some Ri was found to have a zero
2617
C both at T0 and very close to T0.
2618
C IRT = 1 means a legitimate root was found (JOB = 2 or 3).
2619
C On return, T0 is the root location, and Y is the
2620
C corresponding solution vector.
2621
c IRT = 2 A zero-crossing surface has detached from zero
2623
C T0 = value of T at one endpoint of interval of interest. Only
2624
C roots beyond T0 in the direction of integration are sought.
2625
C T0 is input if JOB .ge. 2, and output in all cases.
2626
C T0 is updated by DRCHEK, whether a root is found or not.
2627
C Stored in the global array RWORK.
2628
C TLAST = last value of T returned by the solver (input only).
2629
C Stored in the global array RWORK.
2630
C TOUT = final output time for the solver.
2631
C IRFND = input flag showing whether the last step taken had a root.
2632
C IRFND = 1 if it did, = 0 if not.
2633
C Stored in the global array IWORK.
2634
C INFO3 = copy of INFO(3) (input only).
2635
C-----------------------------------------------------------------------
2639
LMASK=IWORK(LNIW)-NRT
2640
HMINR = (ABS(TN) + ABS(H))*UROUND*100.0D0
2642
GO TO (100, 200, 300), JOB
2644
C Evaluate R at initial T (= RWORK(LT0)); check for zero values.--------
2648
103 IWORK(LMASK+I)=0
2649
CALL DDATRP1(TN,RWORK(LT0),Y,YP,NEQ,KOLD,PHI,PSI)
2650
CALL RT (NEQ, RWORK(LT0), Y, YP, NRT, R0, RPAR, IPAR)
2654
IF (DABS(R0(I)) .EQ. ZERO) THEN
2660
IF (.NOT. ZROOT) GO TO 190
2661
C R has a zero at T. Look at R at T + (small increment). --------------
2662
TEMP1 = SIGN(HMINR,H)
2663
RWORK(LT0) = RWORK(LT0) + TEMP1
2666
120 Y(I) = Y(I) + TEMP2*PHI(I,2)
2667
CALL RT (NEQ, RWORK(LT0), Y, YP, NRT, R0, RPAR, IPAR)
2668
IWORK(LNRTE) = IWORK(LNRTE) + 1
2669
C --------- MASKING THE STUCK ZEROS IN COLD-RESTARTS
2672
IF (JROOT(I) .EQ. 1) THEN
2673
IF (ABS(R0(I)) .EQ. ZERO) THEN
2676
c . to take one step through DDSTP then then in the next arrival->exit!
2678
JROOT(I)=SIGN(2.0D0,R0(I))
2687
c in the previous call there was not a root, so this part can be ignored.
2688
IF (IWORK(LIRFND) .EQ. 0) GO TO 260
2689
c --------------- INITIALIZING THE MASKS
2692
203 IWORK(LMASK+I)=0
2693
C If a root was found on the previous step, evaluate R0 = R(T0). -------
2694
CALL DDATRP1 (TN, RWORK(LT0), Y, YP, NEQ, KOLD, PHI, PSI)
2695
CALL RT (NEQ, RWORK(LT0), Y, YP, NRT, R0, RPAR, IPAR)
2696
IWORK(LNRTE) = IWORK(LNRTE) + 1
2699
IF (ABS(R0(I)) .EQ. ZERO) THEN
2705
IF (.NOT. ZROOT) GO TO 260
2706
C R has a zero at T0. Look at R at T0+ = T0 + (small increment). ------
2707
TEMP1 = SIGN(HMINR,H)
2708
RWORK(LT0) = RWORK(LT0) + TEMP1
2709
IF ((RWORK(LT0) - TN)*H .LT. ZERO) GO TO 230
2712
Y(I) = Y(I) + TEMP2*PHI(I,2)
2715
230 CALL DDATRP1 (TN, RWORK(LT0), Y, YP, NEQ, KOLD, PHI, PSI)
2716
240 CALL RT (NEQ, RWORK(LT0), Y, YP, NRT, R0, RPAR, IPAR)
2717
IWORK(LNRTE) = IWORK(LNRTE) + 1
2719
IF (ABS(R0(I)) .GT. ZERO) GO TO 250
2720
C If Ri has a zero at both T0+ and T0, return an error flag. -----------
2721
IF (JROOT(I) .EQ. 1) THEN
2722
C . MASKING THE STUCK ZEROS IN HOT-RESTARTS
2726
C If Ri has a zero at T0+, but not at T0, return valid root. -----------
2727
JROOT(I) = -SIGN(1.0D0,R0(I))
2732
IF (IRT .EQ. 1) RETURN
2733
C R0 has no zero components. Proceed to check relevant interval. ------
2734
260 IF (TN .EQ. RWORK(LTLAST)) RETURN
2737
C AT THE BEGINING OF THE PREVIOUS STEP THERE WERE A MASK-LIFTING
2738
IF (IWORK(LIRFND) .EQ. 2) THEN
2743
C Set T1 to TN or TOUT, whichever comes first, and get R at T1. --------
2744
IF (INFO3 .EQ. 1 .OR. (TOUT - TN)*H .GE. ZERO) THEN
2749
IF ((T1 - RWORK(LT0))*H .LE. ZERO) RETURN
2750
330 CALL DDATRP1 (TN, T1, Y, YP, NEQ, KOLD, PHI, PSI)
2751
CALL RT (NEQ, T1, Y, YP, NRT, R1, RPAR, IPAR)
2752
IWORK(LNRTE) = IWORK(LNRTE) + 1
2753
C . PASSING THE MASK THROUGH DROOT2 VIA JROOT(I)
2756
331 IF(IWORK(LMASK+I).EQ.1) jroot(i)=1
2757
C Call DROOTS to search for root in interval from T0 to T1. -----------
2760
CALL DROOTS1(NRT, HMINR, JFLAG,RWORK(LT0),T1, R0,R1,RX, X, JROOT)
2761
IF (JFLAG .GT. 1) GO TO 360
2762
CALL DDATRP1 (TN, X, Y, YP, NEQ, KOLD, PHI, PSI)
2763
CALL RT (NEQ, X, Y, YP, NRT, RX, RPAR, IPAR)
2764
IWORK(LNRTE) = IWORK(LNRTE) + 1
2767
CALL DCOPY (NRT, RX, 1, R0, 1)
2768
IF (JFLAG .NE. 4) THEN
2769
CALL DDATRP1 (TN, X, Y, YP, NEQ, KOLD, PHI, PSI)
2773
IF(ABS(JROOT(I)).EQ.1) THEN
2781
IF(ABS(JROOT(I)).EQ.2) JROOT(I)=0
2789
C---------------------- END OF SUBROUTINE DRCHE -----------------------
2791
SUBROUTINE DROOTS1(NRT, HMIN, JFLAG, X0, X1, R0, R1, RX, X, JROOT)
2793
C***BEGIN PROLOGUE DROOTS
2795
C***ROUTINES CALLED DCOPY
2796
C***REVISION HISTORY (YYMMDD)
2797
C 020815 DATE WRITTEN
2798
C 021217 Added root direction information in JROOT.
2799
C***END PROLOGUE DROOTS
2801
INTEGER NRT, JFLAG, JROOT
2802
DOUBLE PRECISION HMIN, X0, X1, R0, R1, RX, X
2803
DIMENSION R0(NRT), R1(NRT), RX(NRT), JROOT(NRT)
2804
C-----------------------------------------------------------------------
2805
C This subroutine finds the leftmost root of a set of arbitrary
2806
C functions Ri(x) (i = 1,...,NRT) in an interval (X0,X1). Only roots
2807
C of odd multiplicity (i.e. changes of sign of the Ri) are found.
2808
C Here the sign of X1 - X0 is arbitrary, but is constant for a given
2809
C problem, and -leftmost- means nearest to X0.
2810
C The values of the vector-valued function R(x) = (Ri, i=1...NRT)
2811
C are communicated through the call sequence of DROOTS.
2812
C The method used is the Illinois algorithm.
2815
C Kathie L. Hiebert and Lawrence F. Shampine, Implicitly Defined
2816
C Output Points for Solutions of ODEs, Sandia Report SAND80-0180,
2819
C Description of parameters.
2821
C NRT = number of functions Ri, or the number of components of
2822
C the vector valued function R(x). Input only.
2824
C HMIN = resolution parameter in X. Input only. When a root is
2825
C found, it is located only to within an error of HMIN in X.
2826
C Typically, HMIN should be set to something on the order of
2827
C 100 * UROUND * MAX(ABS(X0),ABS(X1)),
2828
C where UROUND is the unit roundoff of the machine.
2830
C JFLAG = integer flag for input and output communication.
2832
C On input, set JFLAG = 0 on the first call for the problem,
2833
C and leave it unchanged until the problem is completed.
2834
C (The problem is completed when JFLAG .ge. 2 on return.)
2836
C On output, JFLAG has the following values and meanings:
2837
C JFLAG = 1 means DROOTS needs a value of R(x). Set RX = R(X)
2838
C and call DROOTS again.
2839
C JFLAG = 2 means a root has been found. The root is
2840
C at X, and RX contains R(X). (Actually, X is the
2841
C rightmost approximation to the root on an interval
2842
C (X0,X1) of size HMIN or less.)
2843
C JFLAG = 3 means X = X1 is a root, with one or more of the Ri
2844
C being zero at X1 and no sign changes in (X0,X1).
2845
C RX contains R(X) on output.
2846
C JFLAG = 4 means no roots (of odd multiplicity) were
2847
C found in (X0,X1) (no sign changes).
2849
C X0,X1 = endpoints of the interval where roots are sought.
2850
C X1 and X0 are input when JFLAG = 0 (first call), and
2851
C must be left unchanged between calls until the problem is
2852
C completed. X0 and X1 must be distinct, but X1 - X0 may be
2853
C of either sign. However, the notion of -left- and -right-
2854
C will be used to mean nearer to X0 or X1, respectively.
2855
C When JFLAG .ge. 2 on return, X0 and X1 are output, and
2856
C are the endpoints of the relevant interval.
2858
C R0,R1 = arrays of length NRT containing the vectors R(X0) and R(X1),
2859
C respectively. When JFLAG = 0, R0 and R1 are input and
2860
C none of the R0(i) should be zero.
2861
C When JFLAG .ge. 2 on return, R0 and R1 are output.
2863
C RX = array of length NRT containing R(X). RX is input
2864
C when JFLAG = 1, and output when JFLAG .ge. 2.
2866
C X = independent variable value. Output only.
2867
C When JFLAG = 1 on output, X is the point at which R(x)
2868
C is to be evaluated and loaded into RX.
2869
C When JFLAG = 2 or 3, X is the root.
2870
C When JFLAG = 4, X is the right endpoint of the interval, X1.
2872
C JROOT = integer array of length NRT. Output only.
2873
C When JFLAG = 2 or 3, JROOT indicates which components
2874
C of R(x) have a root at X, and the direction of the sign
2875
C change across the root in the direction of integration.
2876
C JROOT(i) = 1 if Ri has a root and changes from - to +.
2877
C JROOT(i) = -1 if Ri has a root and changes from + to -.
2878
C Otherwise JROOT(i) = 0.
2879
C-----------------------------------------------------------------------
2880
INTEGER I, IMAX, IMXOLD, LAST, NXLAST,ISTUCK,IUNSTUCK
2881
DOUBLE PRECISION ALPHA, T2, TMAX, X2, ZERO,FRACINT,FRACSUB,TENTH
2883
LOGICAL ZROOT, SGNCHG, XROOT
2884
SAVE ALPHA, X2, IMAX, LAST
2885
DATA ZERO/0.0D0/, TENTH/0.1D0/, HALF/0.5D0/, FIVE/5.0D0/
2888
IF (JFLAG .EQ. 1) GO TO 200
2889
C JFLAG .ne. 1. Check for change in sign of R or zero at X1. ----------
2896
if ((jroot(i) .eq. 1).AND.((ABS(R1(I)) .GT. ZERO))) IUNSTUCK=I
2897
IF (ABS(R1(I)) .GT. ZERO) GO TO 110
2898
if (jroot(i) .eq. 1) GOTO 120
2901
C At this point, R0(i) has been checked and cannot be zero. ------------
2902
110 IF (SIGN(1.0D0,R0(I)) .EQ. SIGN(1.0D0,R1(I))) GO TO 120
2903
T2 = ABS(R1(I)/(R1(I)-R0(I)))
2904
IF (T2 .LE. TMAX) GO TO 120
2908
IF (IMAX .GT. 0) GO TO 130
2910
IF (IMAX .GT. 0) GO TO 130
2912
IF (IMAX .GT. 0) GO TO 130
2917
140 IF (.NOT. SGNCHG) GO TO 420
2918
C There is a sign change. Find the first root in the interval. --------
2923
C Repeat until the first root in the interval is found. Loop point. ---
2925
IF (XROOT) GO TO 300
2926
IF (NXLAST .EQ. LAST) GO TO 160
2929
160 IF (LAST .EQ. 0) GO TO 170
2932
170 ALPHA = 2.0D0*ALPHA
2933
180 if((ABS(R0(IMAX)).EQ.ZERO).OR.(ABS(R1(IMAX)).EQ.ZERO)) THEN
2934
X2=(X0+ALPHA*X1)/(1+ALPHA)
2936
X2 = X1 - (X1-X0)*R1(IMAX)/(R1(IMAX) - ALPHA*R0(IMAX))
2938
c----------------------- Hindmarsh ----------------
2939
c I recently studied the rootfinding algorithm in some detail, and
2940
c found that there is a high potential for an infinite loop within
2941
c the subroutine DROOTS/SROOTS. This is caused by an adjustment to
2942
c the new computed iterate, called X2 there at statement 180. The
2943
c adjustment following 180 is faulty, and should be replaced as
2944
c follows. This logic moves X2 away from one endpoint of the current
2945
c interval bracketing the root if it is too close, but in a way that
2946
c cannot result in an infinite loop. Even if you have not
2947
c encountered any trouble at this spot in DASKR, I recommend you
2949
cc IF ((ABS(X2-X0) .LT. HMIN) .AND.
2950
cc 1 (ABS(X1-X0) .GT. 10.0D0*HMIN)) X2 = X0 + 0.1D0*(X1-X0)
2951
IF (ABS(X2 - X0) < HALF*HMIN) THEN
2952
FRACINT = ABS(X1 - X0)/HMIN
2953
IF (FRACINT .GT. FIVE) THEN
2956
FRACSUB = HALF/FRACINT
2958
X2 = X0 + FRACSUB*(X1 - X0)
2961
IF (ABS(X1 - X2) < HALF*HMIN) THEN
2962
FRACINT = ABS(X1 - X0)/HMIN
2963
IF (FRACINT .GT. FIVE) THEN
2966
FRACSUB = HALF/FRACINT
2968
X2 = X1 - FRACSUB*(X1 - X0)
2970
c----------------------- Hindmarsh ----------------
2973
C Return to the calling routine to get a value of RX = R(X). ----
2975
C Check to see in which interval R changes sign. ----------------
2983
if ((jroot(i) .eq. 1).AND.((ABS(RX(I)) .GT. ZERO))) IUNSTUCK=I
2984
IF (ABS(RX(I)) .GT. ZERO) GO TO 210
2985
if (jroot(i) .eq. 1) go to 220
2988
C Neither R0(i) nor RX(i) can be zero at this point. -------------------
2989
210 IF (SIGN(1.0D0,R0(I)) .EQ. SIGN(1.0D0,RX(I))) GO TO 220
2990
T2 = ABS(RX(I)/(RX(I) - R0(I)))
2991
IF (T2 .LE. TMAX) GO TO 220
2995
IF (IMAX .GT. 0) GO TO 230
2997
IF (IMAX .GT. 0) GO TO 230
2999
IF (IMAX .GT. 0) GO TO 230
3005
IF (.NOT. SGNCHG) GO TO 260
3006
C Sign change between X0 and X2, so replace X1 with X2. ----------------
3008
CALL DCOPY (NRT, RX, 1, R1, 1)
3014
CALL DCOPY (NRT, RX, 1, R0, 1)
3018
270 IF (ABS(X1-X0) .LE. HMIN) XROOT = .TRUE.
3022
C Return with X1 as the root. Set JROOT. Set X = X1 and RX = R1. -----
3025
CALL DCOPY (NRT, R1, 1, RX, 1)
3027
if(jroot(i).eq.1) then
3028
if(ABS(R1(i)).ne. ZERO) THEN
3029
JROOT(I)=SIGN(2.0D0,R1(I))
3034
IF (ABS(R1(I)) .EQ. ZERO) THEN
3035
JROOT(I) = -SIGN(1.0D0,R0(I))
3037
IF (SIGN(1.0D0,R0(I)) .NE. SIGN(1.0D0,R1(I))) THEN
3038
JROOT(I) = SIGN(1.0D0,R1(I) - R0(I))
3046
C No sign changes in this interval. Set X = X1, return JFLAG = 4. -----
3047
420 CALL DCOPY (NRT, R1, 1, RX, 1)
3051
C----------------------- END OF SUBROUTINE DROOTS ----------------------
3054
SUBROUTINE DRCHEK0 (JOB, RT, NRT, NEQ, TN, TOUT, Y, YP, PHI, PSI,
3055
* KOLD, R0, R1, RX, JROOT, IRT, UROUND, INFO3, RWORK, IWORK,
3058
C***BEGIN PROLOGUE DRCHEK
3060
C***ROUTINES CALLED DDATRP, DROOTS, DCOPY, RT
3061
C***REVISION HISTORY (YYMMDD)
3062
C 020815 DATE WRITTEN
3063
C 021217 Added test for roots close when JOB = 2.
3064
C***END PROLOGUE DRCHEK
3066
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3067
C Pointers into IWORK:
3068
PARAMETER (LNRTE=36, LIRFND=37)
3069
C Pointers into RWORK:
3070
PARAMETER (LT0=51, LTLAST=52)
3072
INTEGER JOB, NRT, NEQ, KOLD, JROOT, IRT, INFO3, IWORK, IPAR
3073
DOUBLE PRECISION TN, TOUT, Y, YP, PHI, PSI, R0, R1, RX, UROUND,
3075
DIMENSION Y(*), YP(*), PHI(NEQ,*), PSI(*),
3076
* R0(*), R1(*), RX(*), JROOT(*), RWORK(*), IWORK(*)
3079
DOUBLE PRECISION HMINR, T1, TEMP1, TEMP2, X, ZERO
3082
C-----------------------------------------------------------------------
3083
C This routine checks for the presence of a root of R(T,Y,Y') in the
3084
C vicinity of the current T, in a manner depending on the
3085
C input flag JOB. It calls subroutine DROOTS to locate the root
3086
C as precisely as possible.
3088
C In addition to variables described previously, DRCHEK
3089
C uses the following for communication..
3090
C JOB = integer flag indicating type of call..
3091
C JOB = 1 means the problem is being initialized, and DRCHEK
3092
C is to look for a root at or very near the initial T.
3093
C JOB = 2 means a continuation call to the solver was just
3094
C made, and DRCHEK is to check for a root in the
3095
C relevant part of the step last taken.
3096
C JOB = 3 means a successful step was just taken, and DRCHEK
3097
C is to look for a root in the interval of the step.
3098
C R0 = array of length NRT, containing the value of R at T = T0.
3099
C R0 is input for JOB .ge. 2 and on output in all cases.
3100
C R1,RX = arrays of length NRT for work space.
3101
C IRT = completion flag..
3102
C IRT = 0 means no root was found.
3103
C IRT = -1 means JOB = 1 and a zero was found both at T0 and
3104
C and very close to T0.
3105
C IRT = -2 means JOB = 2 and some Ri was found to have a zero
3106
C both at T0 and very close to T0.
3107
C IRT = 1 means a legitimate root was found (JOB = 2 or 3).
3108
C On return, T0 is the root location, and Y is the
3109
C corresponding solution vector.
3110
C T0 = value of T at one endpoint of interval of interest. Only
3111
C roots beyond T0 in the direction of integration are sought.
3112
C T0 is input if JOB .ge. 2, and output in all cases.
3113
C T0 is updated by DRCHEK, whether a root is found or not.
3114
C Stored in the global array RWORK.
3115
C TLAST = last value of T returned by the solver (input only).
3116
C Stored in the global array RWORK.
3117
C TOUT = final output time for the solver.
3118
C IRFND = input flag showing whether the last step taken had a root.
3119
C IRFND = 1 if it did, = 0 if not.
3120
C Stored in the global array IWORK.
3121
C INFO3 = copy of INFO(3) (input only).
3122
C-----------------------------------------------------------------------
3128
HMINR = (ABS(TN) + ABS(H))*UROUND*100.0D0
3130
GO TO (100, 200, 300), JOB
3132
C Evaluate R at initial T (= RWORK(LT0)); check for zero values.--------
3134
CALL DDATRP1(TN,RWORK(LT0),Y,YP,NEQ,KOLD,PHI,PSI)
3135
CALL RT (NEQ, RWORK(LT0), Y, YP, NRT, R0, RPAR, IPAR)
3139
110 IF (ABS(R0(I)) .EQ. ZERO) ZROOT = .TRUE.
3140
IF (.NOT. ZROOT) GO TO 190
3141
C R has a zero at T. Look at R at T + (small increment). --------------
3142
TEMP1 = SIGN(HMINR,H)
3143
RWORK(LT0) = RWORK(LT0) + TEMP1
3146
120 Y(I) = Y(I) + TEMP2*PHI(I,2)
3147
CALL RT (NEQ, RWORK(LT0), Y, YP, NRT, R0, RPAR, IPAR)
3148
IWORK(LNRTE) = IWORK(LNRTE) + 1
3151
130 IF (ABS(R0(I)) .EQ. ZERO) ZROOT = .TRUE.
3152
IF (.NOT. ZROOT) GO TO 190
3153
C R has a zero at T and also close to T. Take error return. -----------
3161
c if in the previous step a z-crossing or mask lifting has occured...
3162
IF (IWORK(LIRFND) .EQ. 0) GO TO 260
3163
C If a root was found on the previous step, evaluate R0 = R(T0). ----
3164
CALL DDATRP1 (TN, RWORK(LT0), Y, YP, NEQ, KOLD, PHI, PSI)
3165
CALL RT (NEQ, RWORK(LT0), Y, YP, NRT, R0, RPAR, IPAR)
3166
IWORK(LNRTE) = IWORK(LNRTE) + 1
3169
IF (ABS(R0(I)) .EQ. ZERO) THEN
3174
IF (.NOT. ZROOT) GO TO 260
3175
C R has a zero at T0. Look at R at T0+ = T0 + (small increment). ------
3176
TEMP1 = SIGN(HMINR,H)
3177
RWORK(LT0) = RWORK(LT0) + TEMP1
3178
IF ((RWORK(LT0) - TN)*H .LT. ZERO) GO TO 230
3181
220 Y(I) = Y(I) + TEMP2*PHI(I,2)
3183
230 CALL DDATRP1 (TN, RWORK(LT0), Y, YP, NEQ, KOLD, PHI, PSI)
3184
240 CALL RT (NEQ, RWORK(LT0), Y, YP, NRT, R0, RPAR, IPAR)
3185
IWORK(LNRTE) = IWORK(LNRTE) + 1
3187
IF (ABS(R0(I)) .GT. ZERO) GO TO 250
3188
C If Ri has a zero at both T0+ and T0, return an error flag. -----------
3189
IF (JROOT(I) .EQ. 1) THEN
3193
C If Ri has a zero at T0+, but not at T0, return valid root. -----------
3194
JROOT(I) = -SIGN(1.0D0,R0(I))
3198
IF (IRT .EQ. 1) RETURN
3199
C R0 has no zero components. Proceed to check relevant interval. ------
3200
260 IF (TN .EQ. RWORK(LTLAST)) RETURN
3203
C Set T1 to TN or TOUT, whichever comes first, and get R at T1. --------
3204
IF (INFO3 .EQ. 1 .OR. (TOUT - TN)*H .GE. ZERO) THEN
3209
IF ((T1 - RWORK(LT0))*H .LE. ZERO) GO TO 390
3210
330 CALL DDATRP1 (TN, T1, Y, YP, NEQ, KOLD, PHI, PSI)
3211
CALL RT (NEQ, T1, Y, YP, NRT, R1, RPAR, IPAR)
3212
IWORK(LNRTE) = IWORK(LNRTE) + 1
3213
C Call DROOTS to search for root in interval from T0 to T1. ------------
3216
CALL DROOTS0(NRT, HMINR, JFLAG,RWORK(LT0),T1, R0,R1,RX, X, JROOT)
3217
IF (JFLAG .GT. 1) GO TO 360
3218
CALL DDATRP1 (TN, X, Y, YP, NEQ, KOLD, PHI, PSI)
3219
CALL RT (NEQ, X, Y, YP, NRT, RX, RPAR, IPAR)
3220
IWORK(LNRTE) = IWORK(LNRTE) + 1
3223
CALL DCOPY (NRT, RX, 1, R0, 1)
3224
IF (JFLAG .EQ. 4) GO TO 390
3225
C Found a root. Interpolate to X and return. --------------------------
3226
CALL DDATRP1 (TN, X, Y, YP, NEQ, KOLD, PHI, PSI)
3232
C---------------------- END OF SUBROUTINE DRCHEK -----------------------
3234
SUBROUTINE DROOTS0(NRT, HMIN, JFLAG, X0, X1, R0, R1, RX, X, JROOT)
3236
C***BEGIN PROLOGUE DROOTS
3238
C***ROUTINES CALLED DCOPY
3239
C***REVISION HISTORY (YYMMDD)
3240
C 020815 DATE WRITTEN
3241
C 021217 Added root direction information in JROOT.
3242
C***END PROLOGUE DROOTS
3244
INTEGER NRT, JFLAG, JROOT
3245
DOUBLE PRECISION HMIN, X0, X1, R0, R1, RX, X
3246
DIMENSION R0(NRT), R1(NRT), RX(NRT), JROOT(NRT)
3247
C-----------------------------------------------------------------------
3248
C This subroutine finds the leftmost root of a set of arbitrary
3249
C functions Ri(x) (i = 1,...,NRT) in an interval (X0,X1). Only roots
3250
C of odd multiplicity (i.e. changes of sign of the Ri) are found.
3251
C Here the sign of X1 - X0 is arbitrary, but is constant for a given
3252
C problem, and -leftmost- means nearest to X0.
3253
C The values of the vector-valued function R(x) = (Ri, i=1...NRT)
3254
C are communicated through the call sequence of DROOTS.
3255
C The method used is the Illinois algorithm.
3258
C Kathie L. Hiebert and Lawrence F. Shampine, Implicitly Defined
3259
C Output Points for Solutions of ODEs, Sandia Report SAND80-0180,
3262
C Description of parameters.
3264
C NRT = number of functions Ri, or the number of components of
3265
C the vector valued function R(x). Input only.
3267
C HMIN = resolution parameter in X. Input only. When a root is
3268
C found, it is located only to within an error of HMIN in X.
3269
C Typically, HMIN should be set to something on the order of
3270
C 100 * UROUND * MAX(ABS(X0),ABS(X1)),
3271
C where UROUND is the unit roundoff of the machine.
3273
C JFLAG = integer flag for input and output communication.
3275
C On input, set JFLAG = 0 on the first call for the problem,
3276
C and leave it unchanged until the problem is completed.
3277
C (The problem is completed when JFLAG .ge. 2 on return.)
3279
C On output, JFLAG has the following values and meanings:
3280
C JFLAG = 1 means DROOTS needs a value of R(x). Set RX = R(X)
3281
C and call DROOTS again.
3282
C JFLAG = 2 means a root has been found. The root is
3283
C at X, and RX contains R(X). (Actually, X is the
3284
C rightmost approximation to the root on an interval
3285
C (X0,X1) of size HMIN or less.)
3286
C JFLAG = 3 means X = X1 is a root, with one or more of the Ri
3287
C being zero at X1 and no sign changes in (X0,X1).
3288
C RX contains R(X) on output.
3289
C JFLAG = 4 means no roots (of odd multiplicity) were
3290
C found in (X0,X1) (no sign changes).
3292
C X0,X1 = endpoints of the interval where roots are sought.
3293
C X1 and X0 are input when JFLAG = 0 (first call), and
3294
C must be left unchanged between calls until the problem is
3295
C completed. X0 and X1 must be distinct, but X1 - X0 may be
3296
C of either sign. However, the notion of -left- and -right-
3297
C will be used to mean nearer to X0 or X1, respectively.
3298
C When JFLAG .ge. 2 on return, X0 and X1 are output, and
3299
C are the endpoints of the relevant interval.
3301
C R0,R1 = arrays of length NRT containing the vectors R(X0) and R(X1),
3302
C respectively. When JFLAG = 0, R0 and R1 are input and
3303
C none of the R0(i) should be zero.
3304
C When JFLAG .ge. 2 on return, R0 and R1 are output.
3306
C RX = array of length NRT containing R(X). RX is input
3307
C when JFLAG = 1, and output when JFLAG .ge. 2.
3309
C X = independent variable value. Output only.
3310
C When JFLAG = 1 on output, X is the point at which R(x)
3311
C is to be evaluated and loaded into RX.
3312
C When JFLAG = 2 or 3, X is the root.
3313
C When JFLAG = 4, X is the right endpoint of the interval, X1.
3315
C JROOT = integer array of length NRT. Output only.
3316
C When JFLAG = 2 or 3, JROOT indicates which components
3317
C of R(x) have a root at X, and the direction of the sign
3318
C change across the root in the direction of integration.
3319
C JROOT(i) = 1 if Ri has a root and changes from - to +.
3320
C JROOT(i) = -1 if Ri has a root and changes from + to -.
3321
C Otherwise JROOT(i) = 0.
3322
C-----------------------------------------------------------------------
3323
INTEGER I, IMAX, IMXOLD, LAST, NXLAST
3324
DOUBLE PRECISION ALPHA, T2, TMAX, X2, ZERO
3325
LOGICAL ZROOT, SGNCHG, XROOT
3326
SAVE ALPHA, X2, IMAX, LAST
3329
IF (JFLAG .EQ. 1) GO TO 200
3330
C JFLAG .ne. 1. Check for change in sign of R or zero at X1. ----------
3335
IF (ABS(R1(I)) .GT. ZERO) GO TO 110
3338
C At this point, R0(i) has been checked and cannot be zero. ------------
3339
110 IF (SIGN(1.0D0,R0(I)) .EQ. SIGN(1.0D0,R1(I))) GO TO 120
3340
T2 = ABS(R1(I)/(R1(I)-R0(I)))
3341
IF (T2 .LE. TMAX) GO TO 120
3345
IF (IMAX .GT. 0) GO TO 130
3349
140 IF (.NOT. SGNCHG) GO TO 400
3350
C There is a sign change. Find the first root in the interval. --------
3355
C Repeat until the first root in the interval is found. Loop point. ---
3357
IF (XROOT) GO TO 300
3358
IF (NXLAST .EQ. LAST) GO TO 160
3361
160 IF (LAST .EQ. 0) GO TO 170
3364
170 ALPHA = 2.0D0*ALPHA
3365
180 X2 = X1 - (X1-X0)*R1(IMAX)/(R1(IMAX) - ALPHA*R0(IMAX))
3366
IF ((ABS(X2-X0) .LT. HMIN) .AND.
3367
1 (ABS(X1-X0) .GT. 10.0D0*HMIN)) X2 = X0 + 0.1D0*(X1-X0)
3370
C Return to the calling routine to get a value of RX = R(X). -----------
3372
C Check to see in which interval R changes sign. -----------------------
3378
IF (ABS(RX(I)) .GT. ZERO) GO TO 210
3381
C Neither R0(i) nor RX(i) can be zero at this point. -------------------
3382
210 IF (SIGN(1.0D0,R0(I)) .EQ. SIGN(1.0D0,RX(I))) GO TO 220
3383
T2 = ABS(RX(I)/(RX(I) - R0(I)))
3384
IF (T2 .LE. TMAX) GO TO 220
3388
IF (IMAX .GT. 0) GO TO 230
3394
IF (.NOT. SGNCHG) GO TO 250
3395
C Sign change between X0 and X2, so replace X1 with X2. ----------------
3397
CALL DCOPY (NRT, RX, 1, R1, 1)
3401
250 IF (.NOT. ZROOT) GO TO 260
3402
C Zero value at X2 and no sign change in (X0,X2), so X2 is a root. -----
3404
CALL DCOPY (NRT, RX, 1, R1, 1)
3407
C No sign change between X0 and X2. Replace X0 with X2. ---------------
3409
CALL DCOPY (NRT, RX, 1, R0, 1)
3413
270 IF (ABS(X1-X0) .LE. HMIN) XROOT = .TRUE.
3416
C Return with X1 as the root. Set JROOT. Set X = X1 and RX = R1. -----
3419
CALL DCOPY (NRT, R1, 1, RX, 1)
3422
IF (ABS(R1(I)) .EQ. ZERO) THEN
3423
JROOT(I) = -SIGN(1.0D0,R0(I))
3426
IF (SIGN(1.0D0,R0(I)) .NE. SIGN(1.0D0,R1(I)))
3427
1 JROOT(I) = SIGN(1.0D0,R1(I) - R0(I))
3431
C No sign change in the interval. Check for zero at right endpoint. ---
3432
400 IF (.NOT. ZROOT) GO TO 420
3434
C Zero value at X1 and no sign change in (X0,X1). Return JFLAG = 3. ---
3436
CALL DCOPY (NRT, R1, 1, RX, 1)
3439
IF (ABS(R1(I)) .EQ. ZERO) JROOT(I) = -SIGN(1.0D0,R0(I))
3444
C No sign changes in this interval. Set X = X1, return JFLAG = 4. -----
3445
420 CALL DCOPY (NRT, R1, 1, RX, 1)
3449
C----------------------- END OF SUBROUTINE DROOTS ----------------------
3452
SUBROUTINE DDASIC (X, Y, YPRIME, NEQ, ICOPT, ID, RES, JAC, PSOL,
3453
* H, TSCALE, WT, NIC, IDID, RPAR, IPAR, PHI, SAVR, DELTA, E,
3454
* YIC, YPIC, PWK, WM, IWM, UROUND, EPLI, SQRTN, RSQRTN,
3455
* EPCONI, STPTOL, JFLG, ICNFLG, ICNSTR, NLSIC)
3457
C***BEGIN PROLOGUE DDASIC
3459
C***DATE WRITTEN 940628 (YYMMDD)
3460
C***REVISION DATE 941206 (YYMMDD)
3461
C***REVISION DATE 950714 (YYMMDD)
3462
C***REVISION DATE 000628 TSCALE argument added.
3464
C-----------------------------------------------------------------------
3467
C DDASIC is a driver routine to compute consistent initial values
3468
C for Y and YPRIME. There are two different options:
3469
C Denoting the differential variables in Y by Y_d, and
3470
C the algebraic variables by Y_a, the problem solved is either:
3471
C 1. Given Y_d, calculate Y_a and Y_d', or
3472
C 2. Given Y', calculate Y.
3473
C In either case, initial values for the given components
3474
C are input, and initial guesses for the unknown components
3475
C must also be provided as input.
3477
C The external routine NLSIC solves the resulting nonlinear system.
3479
C The parameters represent
3481
C X -- Independent variable.
3482
C Y -- Solution vector at X.
3483
C YPRIME -- Derivative of solution vector.
3484
C NEQ -- Number of equations to be integrated.
3485
C ICOPT -- Flag indicating initial condition option chosen.
3486
C ICOPT = 1 for option 1 above.
3487
C ICOPT = 2 for option 2.
3488
C ID -- Array of dimension NEQ, which must be initialized
3489
C if option 1 is chosen.
3490
C ID(i) = +1 if Y_i is a differential variable,
3491
C ID(i) = -1 if Y_i is an algebraic variable.
3492
C RES -- External user-supplied subroutine to evaluate the
3493
C residual. See RES description in DDASPK prologue.
3494
C JAC -- External user-supplied routine to update Jacobian
3495
C or preconditioner information in the nonlinear solver
3496
C (optional). See JAC description in DDASPK prologue.
3497
C PSOL -- External user-supplied routine to solve
3498
C a linear system using preconditioning.
3499
C See PSOL in DDASPK prologue.
3500
C H -- Scaling factor in iteration matrix. DDASIC may
3501
C reduce H to achieve convergence.
3502
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
3503
C WT -- Vector of weights for error criterion.
3504
C NIC -- Input number of initial condition calculation call
3506
C IDID -- Completion code. See IDID in DDASPK prologue.
3507
C RPAR,IPAR -- Real and integer parameter arrays that
3508
C are used for communication between the
3509
C calling program and external user routines.
3510
C They are not altered by DNSK
3511
C PHI -- Work space for DDASIC of length at least 2*NEQ.
3512
C SAVR -- Work vector for DDASIC of length NEQ.
3513
C DELTA -- Work vector for DDASIC of length NEQ.
3514
C E -- Work vector for DDASIC of length NEQ.
3515
C YIC,YPIC -- Work vectors for DDASIC, each of length NEQ.
3516
C PWK -- Work vector for DDASIC of length NEQ.
3517
C WM,IWM -- Real and integer arrays storing
3518
C information required by the linear solver.
3519
C EPCONI -- Test constant for Newton iteration convergence.
3520
C ICNFLG -- Flag showing whether constraints on Y are to apply.
3521
C ICNSTR -- Integer array of length NEQ with constraint types.
3523
C The other parameters are for use internally by DDASIC.
3525
C-----------------------------------------------------------------------
3529
C***END PROLOGUE DDASIC
3532
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3533
DIMENSION Y(*),YPRIME(*),ID(*),WT(*),PHI(NEQ,*)
3534
DIMENSION SAVR(*),DELTA(*),E(*),YIC(*),YPIC(*),PWK(*)
3535
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*), ICNSTR(*)
3536
EXTERNAL RES, JAC, PSOL, NLSIC
3539
PARAMETER (LMXNH=34)
3541
C The following parameters are data-loaded here:
3542
C RHCUT = factor by which H is reduced on retry of Newton solve.
3543
C RATEMX = maximum convergence rate for which Newton iteration
3544
C is considered converging.
3547
DATA RHCUT/0.1D0/, RATEMX/0.8D0/
3550
C-----------------------------------------------------------------------
3553
C JSKIP is a flag set to 1 when NIC = 2 and NH = 1, to signal that
3554
C the initial call to the JAC routine is to be skipped then.
3555
C Save Y and YPRIME in PHI. Initialize IDID, NH, and CJ.
3556
C-----------------------------------------------------------------------
3562
IF (NIC .EQ. 2) JSKIP = 1
3563
CALL DCOPY (NEQ, Y, 1, PHI(1,1), 1)
3564
CALL DCOPY (NEQ, YPRIME, 1, PHI(1,2), 1)
3566
IF (ICOPT .EQ. 2) THEN
3572
C-----------------------------------------------------------------------
3574
C Call the nonlinear system solver to obtain
3575
C consistent initial values for Y and YPRIME.
3576
C-----------------------------------------------------------------------
3580
CALL NLSIC(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JAC,PSOL,H,TSCALE,WT,
3581
* JSKIP,RPAR,IPAR,SAVR,DELTA,E,YIC,YPIC,PWK,WM,IWM,CJ,UROUND,
3582
* EPLI,SQRTN,RSQRTN,EPCONI,RATEMX,STPTOL,JFLG,ICNFLG,ICNSTR,
3585
IF (IERNLS .EQ. 0) RETURN
3587
C-----------------------------------------------------------------------
3589
C The nonlinear solver was unsuccessful. Increment NCFN.
3590
C Return with IDID = -12 if either
3591
C IERNLS = -1: error is considered unrecoverable,
3592
C ICOPT = 2: we are doing initialization problem type 2, or
3593
C NH = MXNH: the maximum number of H values has been tried.
3594
C Otherwise (problem 1 with IERNLS .GE. 1), reduce H and try again.
3595
C If IERNLS > 1, restore Y and YPRIME to their original values.
3596
C-----------------------------------------------------------------------
3598
IWM(LCFN) = IWM(LCFN) + 1
3601
IF (IERNLS .EQ. -1) GO TO 350
3602
c >>>>>>>> singular Jacobian
3603
IF (IERNLS .EQ. -2) GO TO 360
3605
IF (ICOPT .EQ. 2) GO TO 350
3606
IF (NH .EQ. MXNH) GO TO 350
3612
IF (IERNLS .EQ. 1) GO TO 200
3614
CALL DCOPY (NEQ, PHI(1,1), 1, Y, 1)
3615
CALL DCOPY (NEQ, PHI(1,2), 1, YPRIME, 1)
3620
c >> singular Jacobian
3624
C------END OF SUBROUTINE DDASIC-----------------------------------------
3626
SUBROUTINE DYYPNW (NEQ, Y, YPRIME, CJ, RL, P, ICOPT, ID,
3629
C***BEGIN PROLOGUE DYYPNW
3631
C***DATE WRITTEN 940830 (YYMMDD)
3634
C-----------------------------------------------------------------------
3637
C DYYPNW calculates the new (Y,YPRIME) pair needed in the
3638
C linesearch algorithm based on the current lambda value. It is
3639
C called by DLINSK and DLINSD. Based on the ICOPT and ID values,
3640
C the corresponding entry in Y or YPRIME is updated.
3642
C In addition to the parameters described in the calling programs,
3643
C the parameters represent
3645
C P -- Array of length NEQ that contains the current
3646
C approximate Newton step.
3647
C RL -- Scalar containing the current lambda value.
3648
C YNEW -- Array of length NEQ containing the updated Y vector.
3649
C YPNEW -- Array of length NEQ containing the updated YPRIME
3651
C-----------------------------------------------------------------------
3653
C***ROUTINES CALLED (NONE)
3655
C***END PROLOGUE DYYPNW
3658
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3659
DIMENSION Y(*), YPRIME(*), YNEW(*), YPNEW(*), ID(*), P(*)
3661
IF (ICOPT .EQ. 1) THEN
3663
IF(ID(I) .LT. 0) THEN
3664
YNEW(I) = Y(I) - RL*P(I)
3665
YPNEW(I) = YPRIME(I)
3668
YPNEW(I) = YPRIME(I) - RL*CJ*P(I)
3673
YNEW(I) = Y(I) - RL*P(I)
3674
YPNEW(I) = YPRIME(I)
3678
C----------------------- END OF SUBROUTINE DYYPNW ----------------------
3680
SUBROUTINE DDSTP(X,Y,YPRIME,NEQ,RES,JAC,PSOL,H,WT,VT,
3681
* JSTART,IDID,RPAR,IPAR,PHI,SAVR,DELTA,E,WM,IWM,
3682
* ALPHA,BETA,GAMMA,PSI,SIGMA,CJ,CJOLD,HOLD,S,HMIN,UROUND,
3683
* EPLI,SQRTN,RSQRTN,EPCON,IPHASE,JCALC,JFLG,K,KOLD,NS,NONNEG,
3686
C***BEGIN PROLOGUE DDSTP2
3688
C***DATE WRITTEN 890101 (YYMMDD)
3689
C***REVISION DATE 900926 (YYMMDD)
3690
C***REVISION DATE 940909 (YYMMDD) (Reset PSI(1), PHI(*,2) at 690)
3693
C-----------------------------------------------------------------------
3696
C DDSTP solves a system of differential/algebraic equations of
3697
C the form G(X,Y,YPRIME) = 0, for one step (normally from X to X+H).
3699
C The methods used are modified divided difference, fixed leading
3700
C coefficient forms of backward differentiation formulas.
3701
C The code adjusts the stepsize and order to control the local error
3705
C The parameters represent
3706
C X -- Independent variable.
3707
C Y -- Solution vector at X.
3708
C YPRIME -- Derivative of solution vector
3709
C after successful step.
3710
C NEQ -- Number of equations to be integrated.
3711
C RES -- External user-supplied subroutine
3712
C to evaluate the residual. See RES description
3713
C in DDASPK prologue.
3714
C JAC -- External user-supplied routine to update
3715
C Jacobian or preconditioner information in the
3716
C nonlinear solver. See JAC description in DDASPK
3718
C PSOL -- External user-supplied routine to solve
3719
C a linear system using preconditioning.
3720
C (This is optional). See PSOL in DDASPK prologue.
3721
C H -- Appropriate step size for next step.
3722
C Normally determined by the code.
3723
C WT -- Vector of weights for error criterion used in Newton test.
3724
C VT -- Masked vector of weights used in error test.
3725
C JSTART -- Integer variable set 0 for
3726
C first step, 1 otherwise.
3727
C IDID -- Completion code returned from the nonlinear solver.
3728
C See IDID description in DDASPK prologue.
3729
C RPAR,IPAR -- Real and integer parameter arrays that
3730
C are used for communication between the
3731
C calling program and external user routines.
3732
C They are not altered by DNSK
3733
C PHI -- Array of divided differences used by
3734
C DDSTP. The length is NEQ*(K+1), where
3735
C K is the maximum order.
3736
C SAVR -- Work vector for DDSTP of length NEQ.
3737
C DELTA,E -- Work vectors for DDSTP of length NEQ.
3738
C WM,IWM -- Real and integer arrays storing
3739
C information required by the linear solver.
3741
C The other parameters are information
3742
C which is needed internally by DDSTP to
3743
C continue from step to step.
3745
C-----------------------------------------------------------------------
3747
C NLS, DDWNRM, DDATRP
3749
C***END PROLOGUE DDSTP
3752
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3753
DIMENSION Y(*),YPRIME(*),WT(*),VT(*)
3754
DIMENSION PHI(NEQ,*),SAVR(*),DELTA(*),E(*)
3755
DIMENSION WM(*),IWM(*)
3756
DIMENSION PSI(*),ALPHA(*),BETA(*),GAMMA(*),SIGMA(*)
3757
DIMENSION RPAR(*),IPAR(*)
3758
EXTERNAL RES, JAC, PSOL, NLS
3760
PARAMETER (LMXORD=3)
3761
PARAMETER (LNST=11, LETF=14, LCFN=15)
3764
C-----------------------------------------------------------------------
3766
C Initialize. On the first call, set
3767
C the order to 1 and initialize
3769
C-----------------------------------------------------------------------
3771
C Initializations for all calls
3776
c //cold or hot start
3777
IF(JSTART .NE. 0) GO TO 120
3779
C If this is the first step, perform
3780
C other initializations
3795
C-----------------------------------------------------------------------
3797
C Compute coefficients of formulas for
3799
C-----------------------------------------------------------------------
3804
IF(H.NE.HOLD.OR.K .NE. KOLD) NS = 0
3805
NS=MIN0(NS+1,KOLD+2)
3807
IF(KP1 .LT. NS)GO TO 230
3817
BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2
3820
SIGMA(I)=(I-1)*SIGMA(I-1)*ALPHA(I)
3821
GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H
3826
C Compute ALPHAS, ALPHA0
3831
ALPHAS = ALPHAS - 1.0D0/I
3832
ALPHA0 = ALPHA0 - ALPHA(I)
3835
C Compute leading coefficient CJ
3840
C Compute variable stepsize error coefficient CK
3842
CK = ABS(ALPHA(KP1) + ALPHAS - ALPHA0)
3843
CK = MAX(CK,ALPHA(KP1))
3845
C Change PHI to PHI STAR
3847
IF(KP1 .LT. NSP1) GO TO 280
3850
260 PHI(I,J)=BETA(J)*PHI(I,J)
3858
C Initialize IDID to 1
3860
C-----------------------------------------------------------------------
3862
C Call the nonlinear system solver to obtain the solution and
3864
C-----------------------------------------------------------------------
3866
CALL NLS(X,Y,YPRIME,NEQ,
3867
* RES,JAC,PSOL,H,WT,JSTART,IDID,RPAR,IPAR,PHI,GAMMA,
3868
* SAVR,DELTA,E,WM,IWM,CJ,CJOLD,CJLAST,S,
3869
* UROUND,EPLI,SQRTN,RSQRTN,EPCON,JCALC,JFLG,KP1,
3870
* NONNEG,NTYPE,IERNLS)
3872
IF(IERNLS .NE. 0)GO TO 600
3873
C-----------------------------------------------------------------------
3875
C Estimate the errors at orders K,K-1,K-2
3876
C as if constant stepsize was used. Estimate
3877
C the local error at order K and test
3878
C whether the current step is successful.
3879
C-----------------------------------------------------------------------
3881
C Estimate errors at orders K,K-1,K-2
3883
ENORM = DDWNRM(NEQ,E,VT,RPAR,IPAR)
3885
ERK = SIGMA(K+1)*ENORM
3889
IF(K .EQ. 1)GO TO 430
3891
405 DELTA(I) = PHI(I,KP1) + E(I)
3892
ERKM1=SIGMA(K)*DDWNRM(NEQ,DELTA,VT,RPAR,IPAR)
3894
IF(K .GT. 2)GO TO 410
3895
IF(TERKM1 .LE. 0.5*TERK)GO TO 420
3899
415 DELTA(I) = PHI(I,K) + DELTA(I)
3900
ERKM2=SIGMA(K-1)*DDWNRM(NEQ,DELTA,VT,RPAR,IPAR)
3901
TERKM2 = (K-1)*ERKM2
3902
IF(MAX(TERKM1,TERKM2).GT.TERK)GO TO 430
3911
C Calculate the local error for the current step
3912
C to see if the step was successful
3916
IF(ERR .GT. 1.0D0)GO TO 600
3922
C-----------------------------------------------------------------------
3924
C The step is successful. Determine
3925
C the best order and stepsize for
3926
C the next step. Update the differences
3927
C for the next step.
3928
C-----------------------------------------------------------------------
3930
IWM(LNST)=IWM(LNST)+1
3936
C Estimate the error at order K+1 unless
3937
C already decided to lower order, or
3938
C already using maximum order, or
3939
C stepsize not constant, or
3940
C order raised in previous step
3942
IF(KNEW.EQ.KM1.OR.K.EQ.IWM(LMXORD))IPHASE=1
3943
IF(IPHASE .EQ. 0)GO TO 545
3944
IF(KNEW.EQ.KM1)GO TO 540
3945
IF(K.EQ.IWM(LMXORD)) GO TO 550
3946
IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550
3948
510 DELTA(I)=E(I)-PHI(I,KP2)
3949
ERKP1 = (1.0D0/(K+2))*DDWNRM(NEQ,DELTA,VT,RPAR,IPAR)
3950
TERKP1 = (K+2)*ERKP1
3953
IF(TERKP1.GE.0.5D0*TERK)GO TO 550
3955
520 IF(TERKM1.LE.MIN(TERK,TERKP1))GO TO 540
3956
IF(TERKP1.GE.TERK.OR.K.EQ.IWM(LMXORD))GO TO 550
3970
C If IPHASE = 0, increase order by one and multiply stepsize by
3979
C Determine the appropriate stepsize for
3984
R=(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
3985
IF(R .LT. 2.0D0) GO TO 555
3988
555 IF(R .GT. 1.0D0) GO TO 560
3989
R = MAX(0.5D0,MIN(0.9D0,R))
3994
C Update differences for next step
3997
IF(KOLD.EQ.IWM(LMXORD))GO TO 585
4002
590 PHI(I,KP1)=PHI(I,KP1)+E(I)
4006
595 PHI(I,J)=PHI(I,J)+PHI(I,J+1)
4014
C-----------------------------------------------------------------------
4017
C The step is unsuccessful. Restore X,PSI,PHI
4018
C Determine appropriate stepsize for
4019
C continuing the integration, or exit with
4020
C an error flag if there have been many
4022
C-----------------------------------------------------------------------
4028
IF(KP1.LT.NSP1)GO TO 630
4032
610 PHI(I,J)=TEMP1*PHI(I,J)
4036
640 PSI(I-1)=PSI(I)-H
4039
C Test whether failure is due to nonlinear solver
4042
IF(IERNLS .EQ. 0)GO TO 660
4043
IWM(LCFN)=IWM(LCFN)+1
4046
C The nonlinear solver failed to converge.
4047
C Determine the cause of the failure and take appropriate action.
4048
C If IERNLS .LT. 0, then return. Otherwise, reduce the stepsize
4049
C and try again, unless too many failures have occurred.
4051
IF (IERNLS .LT. 0) GO TO 675
4055
IF (NCF .LT. 10 .AND. ABS(H) .GE. HMIN) GO TO 690
4056
IF (IDID .EQ. 1) IDID = -7
4057
IF (NEF .GE. 3) IDID = -9
4061
C The nonlinear solver converged, and the cause
4062
C of the failure was the error estimate
4063
C exceeding the tolerance.
4066
IWM(LETF)=IWM(LETF)+1
4067
IF (NEF .GT. 1) GO TO 665
4069
C On first error test failure, keep current order or lower
4070
C order by one. Compute new stepsize based on differences
4075
R = 0.90D0*(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
4076
R = MAX(0.25D0,MIN(0.9D0,R))
4078
c ------------------ HMIN chnage---------------------
4079
c IF (ABS(H) .GE. HMIN) GO TO 690
4080
if (X+H .GT. X) GO TO 690
4081
c ------------------ HMIN chnage---------------------
4085
C On second error test failure, use the current order or
4086
C decrease order by one. Reduce the stepsize by a factor of
4089
665 IF (NEF .GT. 2) GO TO 670
4093
c ------------------ HMIN chnage---------------------
4094
c IF (ABS(H) .GE. HMIN) GO TO 690
4095
if (X+H .GT. X) GO TO 690
4096
c ------------------ HMIN chnage---------------------
4100
C On third and subsequent error test failures, set the order to
4101
C one, and reduce the stepsize by a factor of one quarter.
4106
c ------------------ HMIN chnage---------------------
4107
c IF (ABS(H) .GE. HMIN) GO TO 690
4108
if (X+H .GT. X) GO TO 690
4109
c ------------------ HMIN chnage---------------------
4116
C For all crashes, restore Y to its last value,
4117
C interpolate to find YPRIME at last X, and return.
4119
C Before returning, verify that the user has not set
4120
C IDID to a nonnegative value. If the user has set IDID
4121
C to a nonnegative value, then reset IDID to be -7, indicating
4122
C a failure in the nonlinear system solver.
4125
CALL DDATRP1(X,X,Y,YPRIME,NEQ,K,PHI,PSI)
4127
IF (IDID .GE. 0) IDID = -7
4131
C Go back and try this step again.
4132
C If this is the first step, reset PSI(1) and rescale PHI(*,2).
4134
690 IF (KOLD .EQ. 0) THEN
4137
695 PHI(I,2) = R*PHI(I,2)
4141
C------END OF SUBROUTINE DDSTP------------------------------------------
4143
SUBROUTINE DCNSTR (NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR)
4145
C***BEGIN PROLOGUE DCNSTR
4146
C***DATE WRITTEN 950808 (YYMMDD)
4147
C***REVISION DATE 950814 (YYMMDD)
4150
C-----------------------------------------------------------------------
4153
C This subroutine checks for constraint violations in the proposed
4154
C new approximate solution YNEW.
4155
C If a constraint violation occurs, then a new step length, TAU,
4156
C is calculated, and this value is to be given to the linesearch routine
4157
C to calculate a new approximate solution YNEW.
4161
C NEQ -- size of the nonlinear system, and the length of arrays
4162
C Y, YNEW and ICNSTR.
4164
C Y -- real array containing the current approximate y.
4166
C YNEW -- real array containing the new approximate y.
4168
C ICNSTR -- INTEGER array of length NEQ containing flags indicating
4169
C which entries in YNEW are to be constrained.
4170
C if ICNSTR(I) = 2, then YNEW(I) must be .GT. 0,
4171
C if ICNSTR(I) = 1, then YNEW(I) must be .GE. 0,
4172
C if ICNSTR(I) = -1, then YNEW(I) must be .LE. 0, while
4173
C if ICNSTR(I) = -2, then YNEW(I) must be .LT. 0, while
4174
C if ICNSTR(I) = 0, then YNEW(I) is not constrained.
4176
C RLX -- real scalar restricting update, if ICNSTR(I) = 2 or -2,
4177
C to ABS( (YNEW-Y)/Y ) < FAC2*RLX in component I.
4179
C TAU -- the current size of the step length for the linesearch.
4183
C TAU -- the adjusted size of the step length if a constraint
4184
C violation occurred (otherwise, it is unchanged). it is
4185
C the step length to give to the linesearch routine.
4187
C IRET -- output flag.
4188
C IRET=0 means that YNEW satisfied all constraints.
4189
C IRET=1 means that YNEW failed to satisfy all the
4190
C constraints, and a new linesearch step
4193
C IVAR -- index of variable causing constraint to be violated.
4195
C-----------------------------------------------------------------------
4196
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4197
DIMENSION Y(NEQ), YNEW(NEQ), ICNSTR(NEQ)
4198
SAVE FAC, FAC2, ZERO
4199
DATA FAC /0.6D0/, FAC2 /0.9D0/, ZERO/0.0D0/
4200
C-----------------------------------------------------------------------
4201
C Check constraints for proposed new step YNEW. If a constraint has
4202
C been violated, then calculate a new step length, TAU, to be
4203
C used in the linesearch routine.
4204
C-----------------------------------------------------------------------
4210
IF (ICNSTR(I) .EQ. 2) THEN
4211
RDY = ABS( (YNEW(I)-Y(I))/Y(I) )
4212
IF (RDY .GT. RDYMX) THEN
4216
IF (YNEW(I) .LE. ZERO) THEN
4223
ELSEIF (ICNSTR(I) .EQ. 1) THEN
4224
IF (YNEW(I) .LT. ZERO) THEN
4231
ELSEIF (ICNSTR(I) .EQ. -1) THEN
4232
IF (YNEW(I) .GT. ZERO) THEN
4239
ELSEIF (ICNSTR(I) .EQ. -2) THEN
4240
RDY = ABS( (YNEW(I)-Y(I))/Y(I) )
4241
IF (RDY .GT. RDYMX) THEN
4245
IF (YNEW(I) .GE. ZERO) THEN
4255
IF(RDYMX .GE. RLX) THEN
4256
TAU = FAC2*TAU*RLX/RDYMX
4261
C----------------------- END OF SUBROUTINE DCNSTR ----------------------
4263
SUBROUTINE DCNST0 (NEQ, Y, ICNSTR, IRET)
4265
C***BEGIN PROLOGUE DCNST0
4266
C***DATE WRITTEN 950808 (YYMMDD)
4267
C***REVISION DATE 950808 (YYMMDD)
4270
C-----------------------------------------------------------------------
4273
C This subroutine checks for constraint violations in the initial
4274
C approximate solution u.
4278
C NEQ -- size of the nonlinear system, and the length of arrays
4281
C Y -- real array containing the initial approximate root.
4283
C ICNSTR -- INTEGER array of length NEQ containing flags indicating
4284
C which entries in Y are to be constrained.
4285
C if ICNSTR(I) = 2, then Y(I) must be .GT. 0,
4286
C if ICNSTR(I) = 1, then Y(I) must be .GE. 0,
4287
C if ICNSTR(I) = -1, then Y(I) must be .LE. 0, while
4288
C if ICNSTR(I) = -2, then Y(I) must be .LT. 0, while
4289
C if ICNSTR(I) = 0, then Y(I) is not constrained.
4293
C IRET -- output flag.
4294
C IRET=0 means that u satisfied all constraints.
4295
C IRET.NE.0 means that Y(IRET) failed to satisfy its
4298
C-----------------------------------------------------------------------
4299
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4300
DIMENSION Y(NEQ), ICNSTR(NEQ)
4303
C-----------------------------------------------------------------------
4304
C Check constraints for initial Y. If a constraint has been violated,
4305
C set IRET = I to signal an error return to calling routine.
4306
C-----------------------------------------------------------------------
4309
IF (ICNSTR(I) .EQ. 2) THEN
4310
IF (Y(I) .LE. ZERO) THEN
4314
ELSEIF (ICNSTR(I) .EQ. 1) THEN
4315
IF (Y(I) .LT. ZERO) THEN
4319
ELSEIF (ICNSTR(I) .EQ. -1) THEN
4320
IF (Y(I) .GT. ZERO) THEN
4324
ELSEIF (ICNSTR(I) .EQ. -2) THEN
4325
IF (Y(I) .GE. ZERO) THEN
4332
C----------------------- END OF SUBROUTINE DCNST0 ----------------------
4334
SUBROUTINE DDAWTS1(NEQ,IWT,RTOL,ATOL,Y,WT,RPAR,IPAR)
4336
C***BEGIN PROLOGUE DDAWTS
4338
C***ROUTINES CALLED (NONE)
4339
C***DATE WRITTEN 890101 (YYMMDD)
4340
C***REVISION DATE 900926 (YYMMDD)
4341
C***END PROLOGUE DDAWTS
4342
C-----------------------------------------------------------------------
4343
C This subroutine sets the error weight vector,
4344
C WT, according to WT(I)=RTOL(I)*ABS(Y(I))+ATOL(I),
4346
C RTOL and ATOL are scalars if IWT = 0,
4347
C and vectors if IWT = 1.
4348
C-----------------------------------------------------------------------
4350
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4351
DIMENSION RTOL(*),ATOL(*),Y(*),WT(*)
4352
DIMENSION RPAR(*),IPAR(*)
4356
IF (IWT .EQ.0) GO TO 10
4359
10 WT(I)=RTOLI*ABS(Y(I))+ATOLI
4363
C------END OF SUBROUTINE DDAWTS-----------------------------------------
4365
SUBROUTINE DINVWT(NEQ,WT,IER)
4367
C***BEGIN PROLOGUE DINVWT
4369
C***ROUTINES CALLED (NONE)
4370
C***DATE WRITTEN 950125 (YYMMDD)
4371
C***END PROLOGUE DINVWT
4372
C-----------------------------------------------------------------------
4373
C This subroutine checks the error weight vector WT, of length NEQ,
4374
C for components that are .le. 0, and if none are found, it
4375
C inverts the WT(I) in place. This replaces division operations
4376
C with multiplications in all norm evaluations.
4377
C IER is returned as 0 if all WT(I) were found positive,
4378
C and the first I with WT(I) .le. 0.0 otherwise.
4379
C-----------------------------------------------------------------------
4381
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4385
IF (WT(I) .LE. 0.0D0) GO TO 30
4388
20 WT(I) = 1.0D0/WT(I)
4395
C------END OF SUBROUTINE DINVWT-----------------------------------------
4397
SUBROUTINE DDATRP1(X,XOUT,YOUT,YPOUT,NEQ,KOLD,PHI,PSI)
4399
C***BEGIN PROLOGUE DDATRP
4401
C***ROUTINES CALLED (NONE)
4402
C***DATE WRITTEN 890101 (YYMMDD)
4403
C***REVISION DATE 900926 (YYMMDD)
4404
C***END PROLOGUE DDATRP
4406
C-----------------------------------------------------------------------
4407
C The methods in subroutine DDSTP use polynomials
4408
C to approximate the solution. DDATRP approximates the
4409
C solution and its derivative at time XOUT by evaluating
4410
C one of these polynomials, and its derivative, there.
4411
C Information defining this polynomial is passed from
4412
C DDSTP, so DDATRP cannot be used alone.
4414
C The parameters are
4416
C X The current time in the integration.
4417
C XOUT The time at which the solution is desired.
4418
C YOUT The interpolated approximation to Y at XOUT.
4420
C YPOUT The interpolated approximation to YPRIME at XOUT.
4422
C NEQ Number of equations.
4423
C KOLD Order used on last successful step.
4424
C PHI Array of scaled divided differences of Y.
4425
C PSI Array of past stepsize history.
4426
C-----------------------------------------------------------------------
4428
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4429
DIMENSION YOUT(*),YPOUT(*)
4430
DIMENSION PHI(NEQ,*),PSI(*)
4440
D=D*GAMMA+C/PSI(J-1)
4442
GAMMA=(TEMP1+PSI(J-1))/PSI(J)
4444
YOUT(I)=YOUT(I)+C*PHI(I,J)
4445
20 YPOUT(I)=YPOUT(I)+D*PHI(I,J)
4449
C------END OF SUBROUTINE DDATRP-----------------------------------------
4451
DOUBLE PRECISION FUNCTION DDWNRM(NEQ,V,RWT,RPAR,IPAR)
4453
C***BEGIN PROLOGUE DDWNRM
4454
C***ROUTINES CALLED (NONE)
4455
C***DATE WRITTEN 890101 (YYMMDD)
4456
C***REVISION DATE 900926 (YYMMDD)
4457
C***END PROLOGUE DDWNRM
4458
C-----------------------------------------------------------------------
4459
C This function routine computes the weighted
4460
C root-mean-square norm of the vector of length
4461
C NEQ contained in the array V, with reciprocal weights
4462
C contained in the array RWT of length NEQ.
4463
C DDWNRM=SQRT((1/NEQ)*SUM(V(I)*RWT(I))**2)
4464
C-----------------------------------------------------------------------
4466
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4467
DIMENSION V(*),RWT(*)
4468
DIMENSION RPAR(*),IPAR(*)
4472
IF(ABS(V(I)*RWT(I)) .GT. VMAX) VMAX = ABS(V(I)*RWT(I))
4474
IF(VMAX .LE. 0.0D0) GO TO 30
4477
20 SUM = SUM + ((V(I)*RWT(I))/VMAX)**2
4479
DDWNRM = VMAX*SQRT(SUM/NEQ)
4483
C------END OF FUNCTION DDWNRM-------------------------------------------
4485
SUBROUTINE DDASID(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JACD,PDUM,H,TSCALE,
4486
* WT,JSDUM,RPAR,IPAR,DUMSVR,DELTA,R,YIC,YPIC,DUMPWK,WM,IWM,CJ,
4487
* UROUND,DUME,DUMS,DUMR,EPCON,RATEMX,STPTOL,JFDUM,
4488
* ICNFLG,ICNSTR,IERNLS)
4490
C***BEGIN PROLOGUE DDASID
4492
C***DATE WRITTEN 940701 (YYMMDD)
4493
C***REVISION DATE 950808 (YYMMDD)
4494
C***REVISION DATE 951110 Removed unreachable block 390.
4495
C***REVISION DATE 000628 TSCALE argument added.
4498
C-----------------------------------------------------------------------
4502
C DDASID solves a nonlinear system of algebraic equations of the
4503
C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in
4504
C the initial conditions.
4506
C The method used is a modified Newton scheme.
4508
C The parameters represent
4510
C X -- Independent variable.
4511
C Y -- Solution vector.
4512
C YPRIME -- Derivative of solution vector.
4513
C NEQ -- Number of unknowns.
4514
C ICOPT -- Initial condition option chosen (1 or 2).
4515
C ID -- Array of dimension NEQ, which must be initialized
4516
C if ICOPT = 1. See DDASIC.
4517
C RES -- External user-supplied subroutine to evaluate the
4518
C residual. See RES description in DDASPK prologue.
4519
C JACD -- External user-supplied routine to evaluate the
4520
C Jacobian. See JAC description for the case
4521
C INFO(12) = 0 in the DDASPK prologue.
4522
C PDUM -- Dummy argument.
4523
C H -- Scaling factor for this initial condition calc.
4524
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
4525
C WT -- Vector of weights for error criterion.
4526
C JSDUM -- Dummy argument.
4527
C RPAR,IPAR -- Real and integer arrays used for communication
4528
C between the calling program and external user
4529
C routines. They are not altered within DASPK.
4530
C DUMSVR -- Dummy argument.
4531
C DELTA -- Work vector for NLS of length NEQ.
4532
C R -- Work vector for NLS of length NEQ.
4533
C YIC,YPIC -- Work vectors for NLS, each of length NEQ.
4534
C DUMPWK -- Dummy argument.
4535
C WM,IWM -- Real and integer arrays storing matrix information
4536
C such as the matrix of partial derivatives,
4537
C permutation vector, and various other information.
4538
C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).
4539
C UROUND -- Unit roundoff.
4540
C DUME -- Dummy argument.
4541
C DUMS -- Dummy argument.
4542
C DUMR -- Dummy argument.
4543
C EPCON -- Tolerance to test for convergence of the Newton
4545
C RATEMX -- Maximum convergence rate for which Newton iteration
4546
C is considered converging.
4547
C JFDUM -- Dummy argument.
4548
C STPTOL -- Tolerance used in calculating the minimum lambda
4550
C ICNFLG -- Integer scalar. If nonzero, then constraint
4551
C violations in the proposed new approximate solution
4552
C will be checked for, and the maximum step length
4553
C will be adjusted accordingly.
4554
C ICNSTR -- Integer array of length NEQ containing flags for
4555
C checking constraints.
4556
C IERNLS -- Error flag for nonlinear solver.
4557
C 0 ==> nonlinear solver converged.
4558
C 1,2 ==> recoverable error inside nonlinear solver.
4559
C 1 => retry with current Y, YPRIME
4560
C 2 => retry with original Y, YPRIME
4561
C -1 ==> unrecoverable error in nonlinear solver.
4562
C -2 ==> Singular Jacobian
4563
C All variables with "DUM" in their names are dummy variables
4564
C which are not used in this routine.
4566
C-----------------------------------------------------------------------
4571
C***END PROLOGUE DDASID
4574
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4575
DIMENSION Y(*),YPRIME(*),ID(*),WT(*),ICNSTR(*)
4576
DIMENSION DELTA(*),R(*),YIC(*),YPIC(*)
4577
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
4580
PARAMETER (LNRE=12, LNJE=13, LMXNIT=32, LMXNJ=33)
4583
C Perform initializations.
4590
C Call RES to initialize DELTA.
4593
IWM(LNRE) = IWM(LNRE) + 1
4594
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
4595
IF (IRES .LT. 0) GO TO 370
4597
C Looping point for updating the Jacobian.
4601
C Initialize all error flags to zero.
4607
C Reevaluate the iteration matrix, J = dG/dY + CJ*dG/dYPRIME,
4608
C where G(X,Y,YPRIME) = 0.
4611
IWM(LNJE)=IWM(LNJE)+1
4613
CALL DMATD(NEQ,X,Y,YPRIME,DELTA,CJ,H,IERJ,WT,R,
4614
* WM,IWM,RES,IRES,UROUND,JACD,RPAR,IPAR)
4615
c assigning two different error message for singular-Jacobian and
4617
IF (IRES .LT. 0) GO TO 370
4618
IF (IERJ .NE. 0) GO TO 375
4620
C Call the nonlinear Newton solver for up to MXNIT iterations.
4622
CALL DNSID(X,Y,YPRIME,NEQ,ICOPT,ID,RES,WT,RPAR,IPAR,DELTA,R,
4623
* YIC,YPIC,WM,IWM,CJ,TSCALE,EPCON,RATEMX,MXNIT,STPTOL,
4624
* ICNFLG,ICNSTR,IERNEW)
4626
IF (IERNEW .EQ. 1 .AND. NJ .LT. MXNJ) THEN
4628
C MXNIT iterations were done, the convergence rate is < 1,
4629
C and the number of Jacobian evaluations is less than MXNJ.
4630
C Call RES, reevaluate the Jacobian, and try again.
4632
IWM(LNRE)=IWM(LNRE)+1
4633
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
4634
IF (IRES .LT. 0) GO TO 370
4637
IF (IERNEW .NE. 0) GO TO 380
4641
C Unsuccessful exits from nonlinear solver.
4642
C Compute IERNLS accordingly.
4644
C unrecoverable error in nonlinear solver.
4647
c >> singular Jacobian
4651
380 IERNLS = MIN(IERNEW,2)
4654
C------END OF SUBROUTINE DDASID-----------------------------------------
4656
SUBROUTINE DNSID(X,Y,YPRIME,NEQ,ICOPT,ID,RES,WT,RPAR,IPAR,
4657
* DELTA,R,YIC,YPIC,WM,IWM,CJ,TSCALE,EPCON,RATEMX,MAXIT,STPTOL,
4658
* ICNFLG,ICNSTR,IERNEW)
4660
C***BEGIN PROLOGUE DNSID
4662
C***DATE WRITTEN 940701 (YYMMDD)
4663
C***REVISION DATE 950713 (YYMMDD)
4664
C***REVISION DATE 000628 TSCALE argument added.
4667
C-----------------------------------------------------------------------
4670
C DNSID solves a nonlinear system of algebraic equations of the
4671
C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME
4672
C in the initial conditions.
4674
C The method used is a modified Newton scheme.
4676
C The parameters represent
4678
C X -- Independent variable.
4679
C Y -- Solution vector.
4680
C YPRIME -- Derivative of solution vector.
4681
C NEQ -- Number of unknowns.
4682
C ICOPT -- Initial condition option chosen (1 or 2).
4683
C ID -- Array of dimension NEQ, which must be initialized
4684
C if ICOPT = 1. See DDASIC.
4685
C RES -- External user-supplied subroutine to evaluate the
4686
C residual. See RES description in DDASPK prologue.
4687
C WT -- Vector of weights for error criterion.
4688
C RPAR,IPAR -- Real and integer arrays used for communication
4689
C between the calling program and external user
4690
C routines. They are not altered within DASPK.
4691
C DELTA -- Residual vector on entry, and work vector of
4692
C length NEQ for DNSID.
4693
C WM,IWM -- Real and integer arrays storing matrix information
4694
C such as the matrix of partial derivatives,
4695
C permutation vector, and various other information.
4696
C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).
4697
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
4698
C R -- Array of length NEQ used as workspace by the
4699
C linesearch routine DLINSD.
4700
C YIC,YPIC -- Work vectors for DLINSD, each of length NEQ.
4701
C EPCON -- Tolerance to test for convergence of the Newton
4703
C RATEMX -- Maximum convergence rate for which Newton iteration
4704
C is considered converging.
4705
C MAXIT -- Maximum allowed number of Newton iterations.
4706
C STPTOL -- Tolerance used in calculating the minimum lambda
4708
C ICNFLG -- Integer scalar. If nonzero, then constraint
4709
C violations in the proposed new approximate solution
4710
C will be checked for, and the maximum step length
4711
C will be adjusted accordingly.
4712
C ICNSTR -- Integer array of length NEQ containing flags for
4713
C checking constraints.
4714
C IERNEW -- Error flag for Newton iteration.
4715
C 0 ==> Newton iteration converged.
4716
C 1 ==> failed to converge, but RATE .le. RATEMX.
4717
C 2 ==> failed to converge, RATE .gt. RATEMX.
4718
C 3 ==> other recoverable error (IRES = -1, or
4719
C linesearch failed).
4720
C -1 ==> unrecoverable error (IRES = -2).
4722
C-----------------------------------------------------------------------
4725
C DSLVD, DDWNRM, DLINSD, DCOPY
4727
C***END PROLOGUE DNSID
4730
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4731
DIMENSION Y(*),YPRIME(*),WT(*),R(*)
4732
DIMENSION ID(*),DELTA(*), YIC(*), YPIC(*)
4733
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
4737
PARAMETER (LNNI=19, LLSOFF=35)
4740
C Initializations. M is the Newton iteration counter.
4748
C Compute a new step vector DELTA by back-substitution.
4750
CALL DSLVD (NEQ, DELTA, WM, IWM)
4753
C Get norm of DELTA. Return now if norm(DELTA) .le. EPCON.
4755
DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
4757
IF (TSCALE .GT. 0.0D0) FNRM = FNRM*TSCALE*ABS(CJ)
4759
IF (FNRM .LE. EPCON) RETURN
4761
C Newton iteration loop.
4764
IWM(LNNI) = IWM(LNNI) + 1
4766
C Call linesearch routine for global strategy and set RATE
4770
CALL DLINSD (NEQ, Y, X, YPRIME, CJ, TSCALE, DELTA, DELNRM, WT,
4771
* LSOFF, STPTOL, IRET, RES, IRES, WM, IWM, FNRM, ICOPT,
4772
* ID, R, YIC, YPIC, ICNFLG, ICNSTR, RLX, RPAR, IPAR)
4776
C Check for error condition from linesearch.
4777
IF (IRET .NE. 0) GO TO 390
4779
C Test for convergence of the iteration, and return or loop.
4780
C ------------------------------------------
4784
C ------------------------------------------
4786
IF (FNRM .LE. EPCON) RETURN
4788
C The iteration has not yet converged. Update M.
4789
C Test whether the maximum number of iterations have been tried.
4792
IF (M .GE. MAXIT) GO TO 380
4794
C Copy the residual to DELTA and its norm to DELNRM, and loop for
4795
C another iteration.
4797
CALL DCOPY (NEQ, R, 1, DELTA, 1)
4801
C The maximum number of iterations was done. Set IERNEW and return.
4804
380 IF (RATE .LE. RATEMX) THEN
4811
390 IF (IRES .LE. -2) THEN
4819
C------END OF SUBROUTINE DNSID------------------------------------------
4821
SUBROUTINE DLINSD (NEQ, Y, T, YPRIME, CJ, TSCALE, P, PNRM, WT,
4822
* LSOFF, STPTOL, IRET, RES, IRES, WM, IWM,
4823
* FNRM, ICOPT, ID, R, YNEW, YPNEW, ICNFLG,
4824
* ICNSTR, RLX, RPAR, IPAR)
4826
C***BEGIN PROLOGUE DLINSD
4828
C***DATE WRITTEN 941025 (YYMMDD)
4829
C***REVISION DATE 941215 (YYMMDD)
4830
C***REVISION DATE 960129 Moved line RL = ONE to top block.
4831
C***REVISION DATE 000628 TSCALE argument added.
4834
C-----------------------------------------------------------------------
4837
C DLINSD uses a linesearch algorithm to calculate a new (Y,YPRIME)
4838
C pair (YNEW,YPNEW) such that
4840
C f(YNEW,YPNEW) .le. (1 - 2*ALPHA*RL)*f(Y,YPRIME) ,
4842
C where 0 < RL <= 1. Here, f(y,y') is defined as
4844
C f(y,y') = (1/2)*norm( (J-inverse)*G(t,y,y') )**2 ,
4846
C where norm() is the weighted RMS vector norm, G is the DAE
4847
C system residual function, and J is the system iteration matrix
4850
C In addition to the parameters defined elsewhere, we have
4852
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
4853
C P -- Approximate Newton step used in backtracking.
4854
C PNRM -- Weighted RMS norm of P.
4855
C LSOFF -- Flag showing whether the linesearch algorithm is
4856
C to be invoked. 0 means do the linesearch, and
4857
C 1 means turn off linesearch.
4858
C STPTOL -- Tolerance used in calculating the minimum lambda
4860
C ICNFLG -- Integer scalar. If nonzero, then constraint violations
4861
C in the proposed new approximate solution will be
4862
C checked for, and the maximum step length will be
4863
C adjusted accordingly.
4864
C ICNSTR -- Integer array of length NEQ containing flags for
4865
C checking constraints.
4866
C RLX -- Real scalar restricting update size in DCNSTR.
4867
C YNEW -- Array of length NEQ used to hold the new Y in
4868
C performing the linesearch.
4869
C YPNEW -- Array of length NEQ used to hold the new YPRIME in
4870
C performing the linesearch.
4871
C Y -- Array of length NEQ containing the new Y (i.e.,=YNEW).
4872
C YPRIME -- Array of length NEQ containing the new YPRIME
4874
C FNRM -- Real scalar containing SQRT(2*f(Y,YPRIME)) for the
4875
C current (Y,YPRIME) on input and output.
4876
C R -- Work array of length NEQ, containing the scaled
4877
C residual (J-inverse)*G(t,y,y') on return.
4878
C IRET -- Return flag.
4879
C IRET=0 means that a satisfactory (Y,YPRIME) was found.
4880
C IRET=1 means that the routine failed to find a new
4881
C (Y,YPRIME) that was sufficiently distinct from
4882
C the current (Y,YPRIME) pair.
4883
C IRET=2 means IRES .ne. 0 from RES.
4884
C-----------------------------------------------------------------------
4887
C DFNRMD, DYYPNW, DCNSTR, DCOPY, XERRWD
4889
C***END PROLOGUE DLINSD
4891
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4893
DIMENSION Y(*), YPRIME(*), WT(*), R(*), ID(*)
4894
DIMENSION WM(*), IWM(*)
4895
DIMENSION YNEW(*), YPNEW(*), P(*), ICNSTR(*)
4896
DIMENSION RPAR(*), IPAR(*)
4899
PARAMETER (LNRE=12, LKPRIN=31)
4901
SAVE ALPHA, ONE, TWO
4902
DATA ALPHA/1.0D-4/, ONE/1.0D0/, TWO/2.0D0/
4906
F1NRM = (FNRM*FNRM)/TWO
4908
IF (KPRIN .GE. 2) THEN
4909
MSG = '------ IN ROUTINE DLINSD-- PNRM = (R1)'
4910
CALL XERRWD(MSG, 38, 901, 0, 0, 0, 0, 1, PNRM, 0.0D0)
4915
C-----------------------------------------------------------------------
4916
C Check for violations of the constraints, if any are imposed.
4917
C If any violations are found, the step vector P is rescaled, and the
4918
C constraint check is repeated, until no violations are found.
4919
C-----------------------------------------------------------------------
4921
IF(ICNFLG .NE. 0) THEN
4923
CALL DYYPNW (NEQ,Y,YPRIME,CJ,RL,P,ICOPT,ID,YNEW,YPNEW)
4924
CALL DCNSTR (NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR)
4925
IF (IRET .EQ. 1) THEN
4927
RATIO = RATIO*RATIO1
4929
20 P(I) = P(I)*RATIO1
4931
IF(KPRIN .GE. 2) THEN
4932
MSG = '------ CONSTRAINT VIOL., PNRM = (R1), INDEX = (I1)'
4933
CALL XERRWD(MSG, 50, 902, 0, 1, IVAR, 0, 1, PNRM, 0.0D0)
4935
IF (PNRM .LE. STPTOL) THEN
4943
SLPI = (-TWO*F1NRM)*RATIO
4945
IF (LSOFF .EQ. 0 .AND. KPRIN .GE. 2) THEN
4946
MSG = '------ MIN. LAMBDA = (R1)'
4947
CALL XERRWD(MSG, 25, 903, 0, 0, 0, 0, 1, RLMIN, 0.0D0)
4949
C-----------------------------------------------------------------------
4950
C Begin iteration to find RL value satisfying alpha-condition.
4951
C If RL becomes less than RLMIN, then terminate with IRET = 1.
4952
C-----------------------------------------------------------------------
4954
c ----------------------------------------------------------
4955
CALL DYYPNW (NEQ,Y,YPRIME,CJ,RL,P,ICOPT,ID,YNEW,YPNEW)
4956
CALL DFNRMD (NEQ, YNEW, T, YPNEW, R, CJ, TSCALE, WT, RES, IRES,
4957
* FNRMP, WM, IWM, RPAR, IPAR)
4959
IWM(LNRE) = IWM(LNRE) + 1
4960
IF (IRES .NE. 0) THEN
4965
IF (LSOFF .EQ. 1) GO TO 150
4967
F1NRMP = FNRMP*FNRMP/TWO
4968
IF (KPRIN .GE. 2) THEN
4969
MSG = '------ LAMBDA = (R1)'
4970
CALL XERRWD(MSG, 20, 904, 0, 0, 0, 0, 1, RL, 0.0D0)
4971
MSG = '------ NORM(F1) = (R1), NORM(F1NEW) = (R2)'
4972
CALL XERRWD(MSG, 43, 905, 0, 0, 0, 0, 2, F1NRM, F1NRMP)
4974
IF (F1NRMP .GT. F1NRM + ALPHA*SLPI*RL) GO TO 200
4975
C-----------------------------------------------------------------------
4976
C Alpha-condition is satisfied, or linesearch is turned off.
4977
C Copy YNEW,YPNEW to Y,YPRIME and return.
4978
C-----------------------------------------------------------------------
4980
CALL DCOPY (NEQ, YNEW, 1, Y, 1)
4981
CALL DCOPY (NEQ, YPNEW, 1, YPRIME, 1)
4983
IF (KPRIN .GE. 1) THEN
4984
MSG = '------ LEAVING ROUTINE DLINSD, FNRM = (R1)'
4985
CALL XERRWD(MSG, 42, 906, 0, 0, 0, 0, 1, FNRM, 0.0D0)
4988
C-----------------------------------------------------------------------
4989
C Alpha-condition not satisfied. Perform backtrack to compute new RL
4990
C value. If no satisfactory YNEW,YPNEW can be found sufficiently
4991
C distinct from Y,YPRIME, then return IRET = 1.
4992
C-----------------------------------------------------------------------
4994
IF (RL .LT. RLMIN) THEN
5002
C-----------------------END OF SUBROUTINE DLINSD ----------------------
5004
SUBROUTINE DFNRMD (NEQ, Y, T, YPRIME, R, CJ, TSCALE, WT,
5005
* RES, IRES, FNORM, WM, IWM, RPAR, IPAR)
5007
C*** BEGIN PROLOGUE DFNRMD
5008
C*** REFER TO DLINSD
5009
C*** DATE WRITTEN 941025 (YYMMDD)
5010
C*** REVISION DATE 000628 TSCALE argument added.
5013
C-----------------------------------------------------------------------
5016
C DFNRMD calculates the scaled preconditioned norm of the nonlinear
5017
C function used in the nonlinear iteration for obtaining consistent
5018
C initial conditions. Specifically, DFNRMD calculates the weighted
5019
C root-mean-square norm of the vector (J-inverse)*G(T,Y,YPRIME),
5020
C where J is the Jacobian matrix.
5022
C In addition to the parameters described in the calling program
5023
C DLINSD, the parameters represent
5025
C R -- Array of length NEQ that contains
5026
C (J-inverse)*G(T,Y,YPRIME) on return.
5027
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
5028
C FNORM -- Scalar containing the weighted norm of R on return.
5029
C-----------------------------------------------------------------------
5032
C RES, DSLVD, DDWNRM
5034
C***END PROLOGUE DFNRMD
5037
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5039
DIMENSION Y(*), YPRIME(*), WT(*), R(*)
5040
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
5041
C-----------------------------------------------------------------------
5043
C-----------------------------------------------------------------------
5045
CALL RES(T,Y,YPRIME,CJ,R,IRES,RPAR,IPAR)
5046
IF (IRES .LT. 0) RETURN
5047
C-----------------------------------------------------------------------
5048
C Apply inverse of Jacobian to vector R.
5049
C-----------------------------------------------------------------------
5050
CALL DSLVD(NEQ,R,WM,IWM)
5051
C-----------------------------------------------------------------------
5052
C Calculate norm of R.
5053
C-----------------------------------------------------------------------
5054
FNORM = DDWNRM(NEQ,R,WT,RPAR,IPAR)
5055
IF (TSCALE .GT. 0.0D0) FNORM = FNORM*TSCALE*ABS(CJ)
5058
C-----------------------END OF SUBROUTINE DFNRMD ----------------------
5060
SUBROUTINE DNEDD(X,Y,YPRIME,NEQ,RES,JACD,PDUM,H,WT,
5061
* JSTART,IDID,RPAR,IPAR,PHI,GAMMA,DUMSVR,DELTA,E,
5062
* WM,IWM,CJ,CJOLD,CJLAST,S,UROUND,DUME,DUMS,DUMR,
5063
* EPCON,JCALC,JFDUM,KP1,NONNEG,NTYPE,IERNLS)
5065
C***BEGIN PROLOGUE DNEDD
5067
C***DATE WRITTEN 891219 (YYMMDD)
5068
C***REVISION DATE 900926 (YYMMDD)
5071
C-----------------------------------------------------------------------
5074
C DNEDD solves a nonlinear system of
5075
C algebraic equations of the form
5076
C G(X,Y,YPRIME) = 0 for the unknown Y.
5078
C The method used is a modified Newton scheme.
5080
C The parameters represent
5082
C X -- Independent variable.
5083
C Y -- Solution vector.
5084
C YPRIME -- Derivative of solution vector.
5085
C NEQ -- Number of unknowns.
5086
C RES -- External user-supplied subroutine
5087
C to evaluate the residual. See RES description
5088
C in DDASPK prologue.
5089
C JACD -- External user-supplied routine to evaluate the
5090
C Jacobian. See JAC description for the case
5091
C INFO(12) = 0 in the DDASPK prologue.
5092
C PDUM -- Dummy argument.
5093
C H -- Appropriate step size for next step.
5094
C WT -- Vector of weights for error criterion.
5095
C JSTART -- Indicates first call to this routine.
5096
C If JSTART = 0, then this is the first call,
5097
C otherwise it is not.
5098
C IDID -- Completion flag, output by DNEDD.
5099
C See IDID description in DDASPK prologue.
5100
C RPAR,IPAR -- Real and integer arrays used for communication
5101
C between the calling program and external user
5102
C routines. They are not altered within DASPK.
5103
C PHI -- Array of divided differences used by
5104
C DNEDD. The length is NEQ*(K+1),where
5105
C K is the maximum order.
5106
C GAMMA -- Array used to predict Y and YPRIME. The length
5107
C is MAXORD+1 where MAXORD is the maximum order.
5108
C DUMSVR -- Dummy argument.
5109
C DELTA -- Work vector for NLS of length NEQ.
5110
C E -- Error accumulation vector for NLS of length NEQ.
5111
C WM,IWM -- Real and integer arrays storing
5112
C matrix information such as the matrix
5113
C of partial derivatives, permutation
5114
C vector, and various other information.
5115
C CJ -- Parameter always proportional to 1/H.
5116
C CJOLD -- Saves the value of CJ as of the last call to DMATD.
5117
C Accounts for changes in CJ needed to
5118
C decide whether to call DMATD.
5119
C CJLAST -- Previous value of CJ.
5120
C S -- A scalar determined by the approximate rate
5121
C of convergence of the Newton iteration and used
5122
C in the convergence test for the Newton iteration.
5124
C If RATE is defined to be an estimate of the
5125
C rate of convergence of the Newton iteration,
5126
C then S = RATE/(1.D0-RATE).
5128
C The closer RATE is to 0., the faster the Newton
5129
C iteration is converging; the closer RATE is to 1.,
5130
C the slower the Newton iteration is converging.
5132
C On the first Newton iteration with an up-dated
5133
C preconditioner S = 100.D0, Thus the initial
5134
C RATE of convergence is approximately 1.
5136
C S is preserved from call to call so that the rate
5137
C estimate from a previous step can be applied to
5139
C UROUND -- Unit roundoff.
5140
C DUME -- Dummy argument.
5141
C DUMS -- Dummy argument.
5142
C DUMR -- Dummy argument.
5143
C EPCON -- Tolerance to test for convergence of the Newton
5145
C JCALC -- Flag used to determine when to update
5146
C the Jacobian matrix. In general:
5148
C JCALC = -1 ==> Call the DMATD routine to update
5149
C the Jacobian matrix.
5150
C JCALC = 0 ==> Jacobian matrix is up-to-date.
5151
C JCALC = 1 ==> Jacobian matrix is out-dated,
5152
C but DMATD will not be called unless
5153
C JCALC is set to -1.
5154
C JFDUM -- Dummy argument.
5155
C KP1 -- The current order(K) + 1; updated across calls.
5156
C NONNEG -- Flag to determine nonnegativity constraints.
5157
C NTYPE -- Identification code for the NLS routine.
5158
C 0 ==> modified Newton; direct solver.
5159
C IERNLS -- Error flag for nonlinear solver.
5160
C 0 ==> nonlinear solver converged.
5161
C 1 ==> recoverable error inside nonlinear solver.
5162
C -1 ==> unrecoverable error inside nonlinear solver.
5164
C All variables with "DUM" in their names are dummy variables
5165
C which are not used in this routine.
5167
C Following is a list and description of local variables which
5168
C may not have an obvious usage. They are listed in roughly the
5169
C order they occur in this subroutine.
5171
C The following group of variables are passed as arguments to
5172
C the Newton iteration solver. They are explained in greater detail
5174
C TOLNEW, MULDEL, MAXIT, IERNEW
5176
C IERTYP -- Flag which tells whether this subroutine is correct.
5177
C 0 ==> correct subroutine.
5178
C 1 ==> incorrect subroutine.
5180
C-----------------------------------------------------------------------
5182
C DDWNRM, RES, DMATD, DNSD
5184
C***END PROLOGUE DNEDD
5187
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
5188
DIMENSION Y(*),YPRIME(*),WT(*)
5189
DIMENSION DELTA(*),E(*)
5190
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
5191
DIMENSION PHI(NEQ,*),GAMMA(*)
5194
PARAMETER (LNRE=12, LNJE=13)
5196
SAVE MULDEL, MAXIT, XRATE
5197
DATA MULDEL/1/, MAXIT/4/, XRATE/0.25D0/
5199
C Verify that this is the correct subroutine.
5202
IF (NTYPE .NE. 0) THEN
5207
C If this is the first step, perform initializations.
5209
IF (JSTART .EQ. 0) THEN
5214
C Perform all other initializations.
5218
C Decide whether new Jacobian is needed.
5220
TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE)
5222
IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1
5223
IF (CJ .NE. CJLAST) S = 100.D0
5225
C-----------------------------------------------------------------------
5226
C Entry point for updating the Jacobian with current
5228
C-----------------------------------------------------------------------
5231
C Initialize all error flags to zero.
5237
C Predict the solution and derivative and compute the tolerance
5238
C for the Newton iteration.
5246
320 YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J)
5248
PNORM = DDWNRM (NEQ,Y,WT,RPAR,IPAR)
5249
TOLNEW = 100.D0*UROUND*PNORM
5251
C Call RES to initialize DELTA.
5253
IWM(LNRE)=IWM(LNRE)+1
5254
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
5255
IF (IRES .LT. 0) GO TO 380
5257
C If indicated, reevaluate the iteration matrix
5258
C J = dG/dY + CJ*dG/dYPRIME (where G(X,Y,YPRIME)=0).
5259
C Set JCALC to 0 as an indicator that this has been done.
5261
IF(JCALC .EQ. -1) THEN
5262
IWM(LNJE)=IWM(LNJE)+1
5264
CALL DMATD(NEQ,X,Y,YPRIME,DELTA,CJ,H,IERJ,WT,E,WM,IWM,
5265
* RES,IRES,UROUND,JACD,RPAR,IPAR)
5268
IF (IRES .LT. 0) GO TO 380
5269
IF(IERJ .NE. 0)GO TO 380
5272
C Call the nonlinear Newton solver.
5274
TEMP1 = 2.0D0/(1.0D0 + CJ/CJOLD)
5275
CALL DNSD(X,Y,YPRIME,NEQ,RES,PDUM,WT,RPAR,IPAR,DUMSVR,
5276
* DELTA,E,WM,IWM,CJ,DUMS,DUMR,DUME,EPCON,S,TEMP1,
5277
* TOLNEW,MULDEL,MAXIT,IRES,IDUM,IERNEW)
5279
IF (IERNEW .GT. 0 .AND. JCALC .NE. 0) THEN
5281
C The Newton iteration had a recoverable failure with an old
5282
C iteration matrix. Retry the step with a new iteration matrix.
5288
IF (IERNEW .NE. 0) GO TO 380
5290
C The Newton iteration has converged. If nonnegativity of
5291
C solution is required, set the solution nonnegative, if the
5292
C perturbation to do it is small enough. If the change is too
5293
C large, then consider the corrector iteration to have failed.
5295
375 IF(NONNEG .EQ. 0) GO TO 390
5297
377 DELTA(I) = MIN(Y(I),0.0D0)
5298
DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
5299
IF(DELNRM .GT. EPCON) GO TO 380
5301
378 E(I) = E(I) - DELTA(I)
5305
C Exits from nonlinear solver.
5306
C No convergence with current iteration
5307
C matrix, or singular iteration matrix.
5308
C Compute IERNLS and IDID accordingly.
5311
IF (IRES .LE. -2 .OR. IERTYP .NE. 0) THEN
5313
IF (IRES .LE. -2) IDID = -11
5314
IF (IERTYP .NE. 0) IDID = -15
5317
IF (IRES .LT. 0) IDID = -10
5318
IF (IERJ .NE. 0) IDID = -8
5324
C------END OF SUBROUTINE DNEDD------------------------------------------
5326
SUBROUTINE DNSD(X,Y,YPRIME,NEQ,RES,PDUM,WT,RPAR,IPAR,
5327
* DUMSVR,DELTA,E,WM,IWM,CJ,DUMS,DUMR,DUME,EPCON,
5328
* S,CONFAC,TOLNEW,MULDEL,MAXIT,IRES,IDUM,IERNEW)
5330
C***BEGIN PROLOGUE DNSD
5332
C***DATE WRITTEN 891219 (YYMMDD)
5333
C***REVISION DATE 900926 (YYMMDD)
5334
C***REVISION DATE 950126 (YYMMDD)
5335
C***REVISION DATE 000711 (YYMMDD)
5338
C-----------------------------------------------------------------------
5341
C DNSD solves a nonlinear system of
5342
C algebraic equations of the form
5343
C G(X,Y,YPRIME) = 0 for the unknown Y.
5345
C The method used is a modified Newton scheme.
5347
C The parameters represent
5349
C X -- Independent variable.
5350
C Y -- Solution vector.
5351
C YPRIME -- Derivative of solution vector.
5352
C NEQ -- Number of unknowns.
5353
C RES -- External user-supplied subroutine
5354
C to evaluate the residual. See RES description
5355
C in DDASPK prologue.
5356
C PDUM -- Dummy argument.
5357
C WT -- Vector of weights for error criterion.
5358
C RPAR,IPAR -- Real and integer arrays used for communication
5359
C between the calling program and external user
5360
C routines. They are not altered within DASPK.
5361
C DUMSVR -- Dummy argument.
5362
C DELTA -- Work vector for DNSD of length NEQ.
5363
C E -- Error accumulation vector for DNSD of length NEQ.
5364
C WM,IWM -- Real and integer arrays storing
5365
C matrix information such as the matrix
5366
C of partial derivatives, permutation
5367
C vector, and various other information.
5368
C CJ -- Parameter always proportional to 1/H (step size).
5369
C DUMS -- Dummy argument.
5370
C DUMR -- Dummy argument.
5371
C DUME -- Dummy argument.
5372
C EPCON -- Tolerance to test for convergence of the Newton
5374
C S -- Used for error convergence tests.
5375
C In the Newton iteration: S = RATE/(1 - RATE),
5376
C where RATE is the estimated rate of convergence
5377
C of the Newton iteration.
5378
C The calling routine passes the initial value
5379
C of S to the Newton iteration.
5380
C CONFAC -- A residual scale factor to improve convergence.
5381
C TOLNEW -- Tolerance on the norm of Newton correction in
5382
C alternative Newton convergence test.
5383
C MULDEL -- A flag indicating whether or not to multiply
5385
C 0 ==> do not scale DELTA by CONFAC.
5386
C 1 ==> scale DELTA by CONFAC.
5387
C MAXIT -- Maximum allowed number of Newton iterations.
5388
C IRES -- Error flag returned from RES. See RES description
5389
C in DDASPK prologue. If IRES = -1, then IERNEW
5391
C If IRES < -1, then IERNEW will be set to -1.
5392
C IDUM -- Dummy argument.
5393
C IERNEW -- Error flag for Newton iteration.
5394
C 0 ==> Newton iteration converged.
5395
C 1 ==> recoverable error inside Newton iteration.
5396
C -1 ==> unrecoverable error inside Newton iteration.
5398
C All arguments with "DUM" in their names are dummy arguments
5399
C which are not used in this routine.
5400
C-----------------------------------------------------------------------
5403
C DSLVD, DDWNRM, RES
5405
C***END PROLOGUE DNSD
5408
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
5409
DIMENSION Y(*),YPRIME(*),WT(*),DELTA(*),E(*)
5410
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
5413
PARAMETER (LNRE=12, LNNI=19)
5415
C Initialize Newton counter M and accumulation vector E.
5424
IWM(LNNI) = IWM(LNNI) + 1
5426
C If necessary, multiply residual by convergence factor.
5428
IF (MULDEL .EQ. 1) THEN
5430
320 DELTA(I) = DELTA(I) * CONFAC
5433
C Compute a new iterate (back-substitution).
5434
C Store the correction in DELTA.
5436
CALL DSLVD(NEQ,DELTA,WM,IWM)
5438
C Update Y, E, and YPRIME.
5443
YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
5447
C Test for convergence of the iteration.
5449
DELNRM=DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
5452
IF (DELNRM .LE. TOLNEW) GO TO 370
5454
RATE = (DELNRM/OLDNRM)**(1.0D0/M)
5455
IF (RATE .GT. 0.9D0) GO TO 380
5456
S = RATE/(1.0D0 - RATE)
5458
IF (S*DELNRM .LE. EPCON) GO TO 370
5460
C The corrector has not yet converged.
5461
C Update M and test whether the
5462
C maximum number of iterations have
5466
IF(M.GE.MAXIT) GO TO 380
5468
C Evaluate the residual,
5469
C and go back to do another iteration.
5471
IWM(LNRE)=IWM(LNRE)+1
5472
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
5473
IF (IRES .LT. 0) GO TO 380
5476
C The iteration has converged.
5480
C The iteration has not converged. Set IERNEW appropriately.
5483
IF (IRES .LE. -2 ) THEN
5491
C------END OF SUBROUTINE DNSD-------------------------------------------
5493
SUBROUTINE DMATD(NEQ,X,Y,YPRIME,DELTA,CJ,H,IER,EWT,E,
5494
* WM,IWM,RES,IRES,UROUND,JACD,RPAR,IPAR)
5496
C***BEGIN PROLOGUE DMATD
5498
C***DATE WRITTEN 890101 (YYMMDD)
5499
C***REVISION DATE 900926 (YYMMDD)
5500
C***REVISION DATE 940701 (YYMMDD) (new LIPVT)
5502
C-----------------------------------------------------------------------
5505
C This routine computes the iteration matrix
5506
C J = dG/dY+CJ*dG/dYPRIME (where G(X,Y,YPRIME)=0).
5507
C Here J is computed by:
5508
C the user-supplied routine JACD if IWM(MTYPE) is 1 or 4, or
5509
C by numerical difference quotients if IWM(MTYPE) is 2 or 5.
5511
C The parameters have the following meanings.
5512
C X = Independent variable.
5513
C Y = Array containing predicted values.
5514
C YPRIME = Array containing predicted derivatives.
5515
C DELTA = Residual evaluated at (X,Y,YPRIME).
5516
C (Used only if IWM(MTYPE)=2 or 5).
5517
C CJ = Scalar parameter defining iteration matrix.
5518
C H = Current stepsize in integration.
5519
C IER = Variable which is .NE. 0 if iteration matrix
5520
C is singular, and 0 otherwise.
5521
C EWT = Vector of error weights for computing norms.
5522
C E = Work space (temporary) of length NEQ.
5523
C WM = Real work space for matrices. On output
5524
C it contains the LU decomposition
5525
C of the iteration matrix.
5526
C IWM = Integer work space containing
5527
C matrix information.
5528
C RES = External user-supplied subroutine
5529
C to evaluate the residual. See RES description
5530
C in DDASPK prologue.
5531
C IRES = Flag which is equal to zero if no illegal values
5532
C in RES, and less than zero otherwise. (If IRES
5533
C is less than zero, the matrix was not completed).
5534
C In this case (if IRES .LT. 0), then IER = 0.
5535
C UROUND = The unit roundoff error of the machine being used.
5536
C JACD = Name of the external user-supplied routine
5537
C to evaluate the iteration matrix. (This routine
5538
C is only used if IWM(MTYPE) is 1 or 4)
5539
C See JAC description for the case INFO(12) = 0
5540
C in DDASPK prologue.
5541
C RPAR,IPAR= Real and integer parameter arrays that
5542
C are used for communication between the
5543
C calling program and external user routines.
5544
C They are not altered by DMATD.
5545
C-----------------------------------------------------------------------
5547
C JACD, RES, DGEFA, DGBFA
5549
C***END PROLOGUE DMATD
5552
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
5553
DIMENSION Y(*),YPRIME(*),DELTA(*),EWT(*),E(*)
5554
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
5557
PARAMETER (LML=1, LMU=2, LMTYPE=4, LNRE=12, LNPD=22, LLCIWP=30)
5562
GO TO (100,200,300,400,500),MTYPE
5565
C Dense user-supplied matrix.
5570
CALL JACD(X,Y,YPRIME,WM,CJ,RPAR,IPAR)
5574
C Dense finite-difference-generated matrix.
5580
DEL=SQUR*MAX(1.0D0,ABS(Y(I)),ABS(H*YPRIME(I)),
5582
DEL=SIGN(DEL,H*YPRIME(I))
5587
YPRIME(I)=YPRIME(I)+CJ*DEL
5588
IWM(LNRE)=IWM(LNRE)+1
5589
CALL RES(X,Y,YPRIME,CJ,E,IRES,RPAR,IPAR)
5590
IF (IRES .LT. 0) RETURN
5593
WM(NROW+L)=(E(L)-DELTA(L))*DELINV
5601
C Do dense-matrix LU decomposition on J.
5603
230 CALL DGEFA(WM,NEQ,NEQ,IWM(LIPVT),IER)
5604
IF (IER .ne. 0) THEN
5605
write(6,'('' Singular Jacobian at IER ='',i3)')IER
5610
C Dummy section for IWM(MTYPE)=3.
5615
C Banded user-supplied matrix.
5620
CALL JACD(X,Y,YPRIME,WM,CJ,RPAR,IPAR)
5621
MEBAND=2*IWM(LML)+IWM(LMU)+1
5625
C Banded finite-difference-generated matrix.
5627
500 MBAND=IWM(LML)+IWM(LMU)+1
5629
MEBAND=MBAND+IWM(LML)
5637
DO 510 N=J,NEQ,MBAND
5640
WM(IPSAVE+K)=YPRIME(N)
5641
DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),
5643
DEL=SIGN(DEL,H*YPRIME(N))
5646
510 YPRIME(N)=YPRIME(N)+CJ*DEL
5647
IWM(LNRE)=IWM(LNRE)+1
5648
CALL RES(X,Y,YPRIME,CJ,E,IRES,RPAR,IPAR)
5649
IF (IRES .LT. 0) RETURN
5650
DO 530 N=J,NEQ,MBAND
5653
YPRIME(N)=WM(IPSAVE+K)
5654
DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),
5656
DEL=SIGN(DEL,H*YPRIME(N))
5659
I1=MAX0(1,(N-IWM(LMU)))
5660
I2=MIN0(NEQ,(N+IWM(LML)))
5663
520 WM(II+I)=(E(I)-DELTA(I))*DELINV
5668
C Do LU decomposition of banded J.
5670
550 CALL DGBFA (WM,MEBAND,NEQ,IWM(LML),IWM(LMU),IWM(LIPVT),IER)
5673
C------END OF SUBROUTINE DMATD------------------------------------------
5675
SUBROUTINE DSLVD(NEQ,DELTA,WM,IWM)
5677
C***BEGIN PROLOGUE DSLVD
5679
C***DATE WRITTEN 890101 (YYMMDD)
5680
C***REVISION DATE 900926 (YYMMDD)
5681
C***REVISION DATE 940701 (YYMMDD) (new LIPVT)
5683
C-----------------------------------------------------------------------
5686
C This routine manages the solution of the linear
5687
C system arising in the Newton iteration.
5688
C Real matrix information and real temporary storage
5689
C is stored in the array WM.
5690
C Integer matrix information is stored in the array IWM.
5691
C For a dense matrix, the LINPACK routine DGESL is called.
5692
C For a banded matrix, the LINPACK routine DGBSL is called.
5693
C-----------------------------------------------------------------------
5697
C***END PROLOGUE DSLVD
5700
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
5701
DIMENSION DELTA(*),WM(*),IWM(*)
5703
PARAMETER (LML=1, LMU=2, LMTYPE=4, LLCIWP=30)
5707
GO TO(100,100,300,400,400),MTYPE
5713
c subroutine dgesl(a,lda,n,ipvt,b,job)
5716
100 CALL DGESL(WM,NEQ,NEQ,IWM(LIPVT),DELTA,0)
5720
C Dummy section for MTYPE=3.
5727
400 MEBAND=2*IWM(LML)+IWM(LMU)+1
5728
CALL DGBSL(WM,MEBAND,NEQ,IWM(LML),
5729
* IWM(LMU),IWM(LIPVT),DELTA,0)
5732
C------END OF SUBROUTINE DSLVD------------------------------------------
5734
SUBROUTINE DDASIK(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JACK,PSOL,H,TSCALE,
5735
* WT,JSKIP,RPAR,IPAR,SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,UROUND,
5736
* EPLI,SQRTN,RSQRTN,EPCON,RATEMX,STPTOL,JFLG,
5737
* ICNFLG,ICNSTR,IERNLS)
5739
C***BEGIN PROLOGUE DDASIK
5741
C***DATE WRITTEN 941026 (YYMMDD)
5742
C***REVISION DATE 950808 (YYMMDD)
5743
C***REVISION DATE 951110 Removed unreachable block 390.
5744
C***REVISION DATE 000628 TSCALE argument added.
5747
C-----------------------------------------------------------------------
5751
C DDASIK solves a nonlinear system of algebraic equations of the
5752
C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in
5753
C the initial conditions.
5755
C An initial value for Y and initial guess for YPRIME are input.
5757
C The method used is a Newton scheme with Krylov iteration and a
5758
C linesearch algorithm.
5760
C The parameters represent
5762
C X -- Independent variable.
5763
C Y -- Solution vector at x.
5764
C YPRIME -- Derivative of solution vector.
5765
C NEQ -- Number of equations to be integrated.
5766
C ICOPT -- Initial condition option chosen (1 or 2).
5767
C ID -- Array of dimension NEQ, which must be initialized
5768
C if ICOPT = 1. See DDASIC.
5769
C RES -- External user-supplied subroutine
5770
C to evaluate the residual. See RES description
5771
C in DDASPK prologue.
5772
C JACK -- External user-supplied routine to update
5773
C the preconditioner. (This is optional).
5774
C See JAC description for the case
5775
C INFO(12) = 1 in the DDASPK prologue.
5776
C PSOL -- External user-supplied routine to solve
5777
C a linear system using preconditioning.
5778
C (This is optional). See explanation inside DDASPK.
5779
C H -- Scaling factor for this initial condition calc.
5780
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
5781
C WT -- Vector of weights for error criterion.
5782
C JSKIP -- input flag to signal if initial JAC call is to be
5783
C skipped. 1 => skip the call, 0 => do not skip call.
5784
C RPAR,IPAR -- Real and integer arrays used for communication
5785
C between the calling program and external user
5786
C routines. They are not altered within DASPK.
5787
C SAVR -- Work vector for DDASIK of length NEQ.
5788
C DELTA -- Work vector for DDASIK of length NEQ.
5789
C R -- Work vector for DDASIK of length NEQ.
5790
C YIC,YPIC -- Work vectors for DDASIK, each of length NEQ.
5791
C PWK -- Work vector for DDASIK of length NEQ.
5792
C WM,IWM -- Real and integer arrays storing
5793
C matrix information for linear system
5794
C solvers, and various other information.
5795
C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).
5796
C UROUND -- Unit roundoff. Not used here.
5797
C EPLI -- convergence test constant.
5798
C See DDASPK prologue for more details.
5799
C SQRTN -- Square root of NEQ.
5800
C RSQRTN -- reciprical of square root of NEQ.
5801
C EPCON -- Tolerance to test for convergence of the Newton
5803
C RATEMX -- Maximum convergence rate for which Newton iteration
5804
C is considered converging.
5805
C JFLG -- Flag showing whether a Jacobian routine is supplied.
5806
C ICNFLG -- Integer scalar. If nonzero, then constraint
5807
C violations in the proposed new approximate solution
5808
C will be checked for, and the maximum step length
5809
C will be adjusted accordingly.
5810
C ICNSTR -- Integer array of length NEQ containing flags for
5811
C checking constraints.
5812
C IERNLS -- Error flag for nonlinear solver.
5813
C 0 ==> nonlinear solver converged.
5814
C 1,2 ==> recoverable error inside nonlinear solver.
5815
C 1 => retry with current Y, YPRIME
5816
C 2 => retry with original Y, YPRIME
5817
C -1 ==> unrecoverable error in nonlinear solver.
5819
C-----------------------------------------------------------------------
5822
C RES, JACK, DNSIK, DCOPY
5824
C***END PROLOGUE DDASIK
5827
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
5828
DIMENSION Y(*),YPRIME(*),ID(*),WT(*),ICNSTR(*)
5829
DIMENSION SAVR(*),DELTA(*),R(*),YIC(*),YPIC(*),PWK(*)
5830
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
5831
EXTERNAL RES, JACK, PSOL
5833
PARAMETER (LNRE=12, LNJE=13, LLOCWP=29, LLCIWP=30)
5834
PARAMETER (LMXNIT=32, LMXNJ=33)
5837
C Perform initializations.
5847
C Call RES to initialize DELTA.
5850
IWM(LNRE) = IWM(LNRE) + 1
5851
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
5852
IF (IRES .LT. 0) GO TO 370
5854
C Looping point for updating the preconditioner.
5858
C Initialize all error flags to zero.
5864
C If a Jacobian routine was supplied, call it.
5866
IF (JFLG .EQ. 1 .AND. JSKIP .EQ. 0) THEN
5868
IWM(LNJE)=IWM(LNJE)+1
5869
CALL JACK (RES, IRES, NEQ, X, Y, YPRIME, WT, DELTA, R, H, CJ,
5870
* WM(LWP), IWM(LIWP), IERPJ, RPAR, IPAR)
5871
IF (IRES .LT. 0 .OR. IERPJ .NE. 0) GO TO 370
5875
C Call the nonlinear Newton solver for up to MXNIT iterations.
5877
CALL DNSIK(X,Y,YPRIME,NEQ,ICOPT,ID,RES,PSOL,WT,RPAR,IPAR,
5878
* SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,TSCALE,SQRTN,RSQRTN,
5879
* EPLIN,EPCON,RATEMX,MXNIT,STPTOL,ICNFLG,ICNSTR,IERNEW)
5881
IF (IERNEW .EQ. 1 .AND. NJ .LT. MXNJ .AND. JFLG .EQ. 1) THEN
5883
C Up to MXNIT iterations were done, the convergence rate is < 1,
5884
C a Jacobian routine is supplied, and the number of JACK calls
5885
C is less than MXNJ.
5886
C Copy the residual SAVR to DELTA, call JACK, and try again.
5888
CALL DCOPY (NEQ, SAVR, 1, DELTA, 1)
5892
IF (IERNEW .NE. 0) GO TO 380
5896
C Unsuccessful exits from nonlinear solver.
5897
C Set IERNLS accordingly.
5900
IF (IRES .LE. -2) IERNLS = -1
5903
380 IERNLS = MIN(IERNEW,2)
5906
C----------------------- END OF SUBROUTINE DDASIK-----------------------
5908
SUBROUTINE DNSIK(X,Y,YPRIME,NEQ,ICOPT,ID,RES,PSOL,WT,RPAR,IPAR,
5909
* SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,TSCALE,SQRTN,RSQRTN,EPLIN,
5910
* EPCON,RATEMX,MAXIT,STPTOL,ICNFLG,ICNSTR,IERNEW)
5912
C***BEGIN PROLOGUE DNSIK
5914
C***DATE WRITTEN 940701 (YYMMDD)
5915
C***REVISION DATE 950714 (YYMMDD)
5916
C***REVISION DATE 000628 TSCALE argument added.
5917
C***REVISION DATE 000628 Added criterion for IERNEW = 1 return.
5920
C-----------------------------------------------------------------------
5923
C DNSIK solves a nonlinear system of algebraic equations of the
5924
C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in
5925
C the initial conditions.
5927
C The method used is a Newton scheme combined with a linesearch
5928
C algorithm, using Krylov iterative linear system methods.
5930
C The parameters represent
5932
C X -- Independent variable.
5933
C Y -- Solution vector.
5934
C YPRIME -- Derivative of solution vector.
5935
C NEQ -- Number of unknowns.
5936
C ICOPT -- Initial condition option chosen (1 or 2).
5937
C ID -- Array of dimension NEQ, which must be initialized
5938
C if ICOPT = 1. See DDASIC.
5939
C RES -- External user-supplied subroutine
5940
C to evaluate the residual. See RES description
5941
C in DDASPK prologue.
5942
C PSOL -- External user-supplied routine to solve
5943
C a linear system using preconditioning.
5944
C See explanation inside DDASPK.
5945
C WT -- Vector of weights for error criterion.
5946
C RPAR,IPAR -- Real and integer arrays used for communication
5947
C between the calling program and external user
5948
C routines. They are not altered within DASPK.
5949
C SAVR -- Work vector for DNSIK of length NEQ.
5950
C DELTA -- Residual vector on entry, and work vector of
5951
C length NEQ for DNSIK.
5952
C R -- Work vector for DNSIK of length NEQ.
5953
C YIC,YPIC -- Work vectors for DNSIK, each of length NEQ.
5954
C PWK -- Work vector for DNSIK of length NEQ.
5955
C WM,IWM -- Real and integer arrays storing
5956
C matrix information such as the matrix
5957
C of partial derivatives, permutation
5958
C vector, and various other information.
5959
C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).
5960
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
5961
C SQRTN -- Square root of NEQ.
5962
C RSQRTN -- reciprical of square root of NEQ.
5963
C EPLIN -- Tolerance for linear system solver.
5964
C EPCON -- Tolerance to test for convergence of the Newton
5966
C RATEMX -- Maximum convergence rate for which Newton iteration
5967
C is considered converging.
5968
C MAXIT -- Maximum allowed number of Newton iterations.
5969
C STPTOL -- Tolerance used in calculating the minimum lambda
5971
C ICNFLG -- Integer scalar. If nonzero, then constraint
5972
C violations in the proposed new approximate solution
5973
C will be checked for, and the maximum step length
5974
C will be adjusted accordingly.
5975
C ICNSTR -- Integer array of length NEQ containing flags for
5976
C checking constraints.
5977
C IERNEW -- Error flag for Newton iteration.
5978
C 0 ==> Newton iteration converged.
5979
C 1 ==> failed to converge, but RATE .lt. 1, or the
5980
C residual norm was reduced by a factor of .1.
5981
C 2 ==> failed to converge, RATE .gt. RATEMX.
5982
C 3 ==> other recoverable error.
5983
C -1 ==> unrecoverable error inside Newton iteration.
5984
C-----------------------------------------------------------------------
5987
C DFNRMK, DSLVK, DDWNRM, DLINSK, DCOPY
5989
C***END PROLOGUE DNSIK
5992
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
5993
DIMENSION Y(*),YPRIME(*),WT(*),ID(*),DELTA(*),R(*),SAVR(*)
5994
DIMENSION YIC(*),YPIC(*),PWK(*),WM(*),IWM(*), RPAR(*),IPAR(*)
5998
PARAMETER (LNNI=19, LNPS=21, LLOCWP=29, LLCIWP=30)
5999
PARAMETER (LLSOFF=35, LSTOL=14)
6002
C Initializations. M is the Newton iteration counter.
6011
C Save residual in SAVR.
6013
CALL DCOPY (NEQ, DELTA, 1, SAVR, 1)
6015
C Compute norm of (P-inverse)*(residual).
6017
CALL DFNRMK (NEQ, Y, X, YPRIME, SAVR, R, CJ, TSCALE, WT,
6018
* SQRTN, RSQRTN, RES, IRES, PSOL, 1, IER, FNRM, EPLIN,
6019
* WM(LWP), IWM(LIWP), PWK, RPAR, IPAR)
6020
IWM(LNPS) = IWM(LNPS) + 1
6021
IF (IER .NE. 0) THEN
6026
C Return now if residual norm is .le. EPCON.
6028
IF (FNRM .LE. EPCON) RETURN
6030
C Newton iteration loop.
6034
IWM(LNNI) = IWM(LNNI) + 1
6036
C Compute a new step vector DELTA.
6038
CALL DSLVK (NEQ, Y, X, YPRIME, SAVR, DELTA, WT, WM, IWM,
6039
* RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK,
6041
IF (IRES .NE. 0 .OR. IERSL .NE. 0) GO TO 390
6043
C Get norm of DELTA. Return now if DELTA is zero.
6045
DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
6046
IF (DELNRM .EQ. 0.0D0) RETURN
6048
C Call linesearch routine for global strategy and set RATE.
6052
CALL DLINSK (NEQ, Y, X, YPRIME, SAVR, CJ, TSCALE, DELTA, DELNRM,
6053
* WT, SQRTN, RSQRTN, LSOFF, STPTOL, IRET, RES, IRES, PSOL,
6054
* WM, IWM, RHOK, FNRM, ICOPT, ID, WM(LWP), IWM(LIWP), R, EPLIN,
6055
* YIC, YPIC, PWK, ICNFLG, ICNSTR, RLX, RPAR, IPAR)
6059
C Check for error condition from linesearch.
6060
IF (IRET .NE. 0) GO TO 390
6062
C Test for convergence of the iteration, and return or loop.
6064
IF (FNRM .LE. EPCON) RETURN
6066
C The iteration has not yet converged. Update M.
6067
C Test whether the maximum number of iterations have been tried.
6070
IF(M .GE. MAXIT) GO TO 380
6072
C Copy the residual SAVR to DELTA and loop for another iteration.
6074
CALL DCOPY (NEQ, SAVR, 1, DELTA, 1)
6077
C The maximum number of iterations was done. Set IERNEW and return.
6079
380 IF (RATE .LE. RATEMX .OR. FNRM .LE. 0.1D0*FNRM0) THEN
6086
390 IF (IRES .LE. -2 .OR. IERSL .LT. 0) THEN
6090
IF (IRES .EQ. 0 .AND. IERSL .EQ. 1 .AND. M .GE. 2
6091
1 .AND. RATE .LT. 1.0D0) IERNEW = 1
6096
C----------------------- END OF SUBROUTINE DNSIK------------------------
6098
SUBROUTINE DLINSK (NEQ, Y, T, YPRIME, SAVR, CJ, TSCALE, P, PNRM,
6099
* WT, SQRTN, RSQRTN, LSOFF, STPTOL, IRET, RES, IRES, PSOL,
6100
* WM, IWM, RHOK, FNRM, ICOPT, ID, WP, IWP, R, EPLIN, YNEW, YPNEW,
6101
* PWK, ICNFLG, ICNSTR, RLX, RPAR, IPAR)
6103
C***BEGIN PROLOGUE DLINSK
6105
C***DATE WRITTEN 940830 (YYMMDD)
6106
C***REVISION DATE 951006 (Arguments SQRTN, RSQRTN added.)
6107
C***REVISION DATE 960129 Moved line RL = ONE to top block.
6108
C***REVISION DATE 000628 TSCALE argument added.
6109
C***REVISION DATE 000628 RHOK*RHOK term removed in alpha test.
6112
C-----------------------------------------------------------------------
6115
C DLINSK uses a linesearch algorithm to calculate a new (Y,YPRIME)
6116
C pair (YNEW,YPNEW) such that
6118
C f(YNEW,YPNEW) .le. (1 - 2*ALPHA*RL)*f(Y,YPRIME)
6120
C where 0 < RL <= 1, and RHOK is the scaled preconditioned norm of
6121
C the final residual vector in the Krylov iteration.
6122
C Here, f(y,y') is defined as
6124
C f(y,y') = (1/2)*norm( (P-inverse)*G(t,y,y') )**2 ,
6126
C where norm() is the weighted RMS vector norm, G is the DAE
6127
C system residual function, and P is the preconditioner used
6128
C in the Krylov iteration.
6130
C In addition to the parameters defined elsewhere, we have
6132
C SAVR -- Work array of length NEQ, containing the residual
6133
C vector G(t,y,y') on return.
6134
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
6135
C P -- Approximate Newton step used in backtracking.
6136
C PNRM -- Weighted RMS norm of P.
6137
C LSOFF -- Flag showing whether the linesearch algorithm is
6138
C to be invoked. 0 means do the linesearch,
6139
C 1 means turn off linesearch.
6140
C STPTOL -- Tolerance used in calculating the minimum lambda
6142
C ICNFLG -- Integer scalar. If nonzero, then constraint violations
6143
C in the proposed new approximate solution will be
6144
C checked for, and the maximum step length will be
6145
C adjusted accordingly.
6146
C ICNSTR -- Integer array of length NEQ containing flags for
6147
C checking constraints.
6148
C RHOK -- Weighted norm of preconditioned Krylov residual.
6149
C RLX -- Real scalar restricting update size in DCNSTR.
6150
C YNEW -- Array of length NEQ used to hold the new Y in
6151
C performing the linesearch.
6152
C YPNEW -- Array of length NEQ used to hold the new YPRIME in
6153
C performing the linesearch.
6154
C PWK -- Work vector of length NEQ for use in PSOL.
6155
C Y -- Array of length NEQ containing the new Y (i.e.,=YNEW).
6156
C YPRIME -- Array of length NEQ containing the new YPRIME
6158
C FNRM -- Real scalar containing SQRT(2*f(Y,YPRIME)) for the
6159
C current (Y,YPRIME) on input and output.
6160
C R -- Work space length NEQ for residual vector.
6161
C IRET -- Return flag.
6162
C IRET=0 means that a satisfactory (Y,YPRIME) was found.
6163
C IRET=1 means that the routine failed to find a new
6164
C (Y,YPRIME) that was sufficiently distinct from
6165
C the current (Y,YPRIME) pair.
6166
C IRET=2 means a failure in RES or PSOL.
6167
C-----------------------------------------------------------------------
6170
C DFNRMK, DYYPNW, DCNSTR, DCOPY, XERRWD
6172
C***END PROLOGUE DLINSK
6174
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
6176
DIMENSION Y(*), YPRIME(*), P(*), WT(*), SAVR(*), R(*), ID(*)
6177
DIMENSION WM(*), IWM(*), YNEW(*), YPNEW(*), PWK(*), ICNSTR(*)
6178
DIMENSION WP(*), IWP(*), RPAR(*), IPAR(*)
6181
PARAMETER (LNRE=12, LNPS=21, LKPRIN=31)
6183
SAVE ALPHA, ONE, TWO
6184
DATA ALPHA/1.0D-4/, ONE/1.0D0/, TWO/2.0D0/
6187
F1NRM = (FNRM*FNRM)/TWO
6190
IF (KPRIN .GE. 2) THEN
6191
MSG = '------ IN ROUTINE DLINSK-- PNRM = (R1)'
6192
CALL XERRWD(MSG, 38, 921, 0, 0, 0, 0, 1, PNRM, 0.0D0)
6196
C-----------------------------------------------------------------------
6197
C Check for violations of the constraints, if any are imposed.
6198
C If any violations are found, the step vector P is rescaled, and the
6199
C constraint check is repeated, until no violations are found.
6200
C-----------------------------------------------------------------------
6201
IF (ICNFLG .NE. 0) THEN
6203
CALL DYYPNW (NEQ,Y,YPRIME,CJ,RL,P,ICOPT,ID,YNEW,YPNEW)
6204
CALL DCNSTR (NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR)
6205
IF (IRET .EQ. 1) THEN
6207
RATIO = RATIO*RATIO1
6209
20 P(I) = P(I)*RATIO1
6211
IF (KPRIN .GE. 2) THEN
6212
MSG = '------ CONSTRAINT VIOL., PNRM = (R1), INDEX = (I1)'
6213
CALL XERRWD (MSG, 50, 922, 0, 1, IVAR, 0, 1, PNRM, 0.0D0)
6215
IF (PNRM .LE. STPTOL) THEN
6223
SLPI = -TWO*F1NRM*RATIO
6225
IF (LSOFF .EQ. 0 .AND. KPRIN .GE. 2) THEN
6226
MSG = '------ MIN. LAMBDA = (R1)'
6227
CALL XERRWD(MSG, 25, 923, 0, 0, 0, 0, 1, RLMIN, 0.0D0)
6229
C-----------------------------------------------------------------------
6230
C Begin iteration to find RL value satisfying alpha-condition.
6231
C Update YNEW and YPNEW, then compute norm of new scaled residual and
6232
C perform alpha condition test.
6233
C-----------------------------------------------------------------------
6235
CALL DYYPNW (NEQ,Y,YPRIME,CJ,RL,P,ICOPT,ID,YNEW,YPNEW)
6236
CALL DFNRMK (NEQ, YNEW, T, YPNEW, SAVR, R, CJ, TSCALE, WT,
6237
* SQRTN, RSQRTN, RES, IRES, PSOL, 0, IER, FNRMP, EPLIN,
6238
* WP, IWP, PWK, RPAR, IPAR)
6239
IWM(LNRE) = IWM(LNRE) + 1
6240
IF (IRES .GE. 0) IWM(LNPS) = IWM(LNPS) + 1
6241
IF (IRES .NE. 0 .OR. IER .NE. 0) THEN
6245
IF (LSOFF .EQ. 1) GO TO 150
6247
F1NRMP = FNRMP*FNRMP/TWO
6248
IF (KPRIN .GE. 2) THEN
6249
MSG = '------ LAMBDA = (R1)'
6250
CALL XERRWD(MSG, 20, 924, 0, 0, 0, 0, 1, RL, 0.0D0)
6251
MSG = '------ NORM(F1) = (R1), NORM(F1NEW) = (R2)'
6252
CALL XERRWD(MSG, 43, 925, 0, 0, 0, 0, 2, F1NRM, F1NRMP)
6254
IF (F1NRMP .GT. F1NRM + ALPHA*SLPI*RL) GO TO 200
6255
C-----------------------------------------------------------------------
6256
C Alpha-condition is satisfied, or linesearch is turned off.
6257
C Copy YNEW,YPNEW to Y,YPRIME and return.
6258
C-----------------------------------------------------------------------
6260
CALL DCOPY(NEQ, YNEW, 1, Y, 1)
6261
CALL DCOPY(NEQ, YPNEW, 1, YPRIME, 1)
6263
IF (KPRIN .GE. 1) THEN
6264
MSG = '------ LEAVING ROUTINE DLINSK, FNRM = (R1)'
6265
CALL XERRWD(MSG, 42, 926, 0, 0, 0, 0, 1, FNRM, 0.0D0)
6268
C-----------------------------------------------------------------------
6269
C Alpha-condition not satisfied. Perform backtrack to compute new RL
6270
C value. If RL is less than RLMIN, i.e. no satisfactory YNEW,YPNEW can
6271
C be found sufficiently distinct from Y,YPRIME, then return IRET = 1.
6272
C-----------------------------------------------------------------------
6274
IF (RL .LT. RLMIN) THEN
6282
C----------------------- END OF SUBROUTINE DLINSK ----------------------
6284
SUBROUTINE DFNRMK (NEQ, Y, T, YPRIME, SAVR, R, CJ, TSCALE, WT,
6285
* SQRTN, RSQRTN, RES, IRES, PSOL, IRIN, IER,
6286
* FNORM, EPLIN, WP, IWP, PWK, RPAR, IPAR)
6288
C***BEGIN PROLOGUE DFNRMK
6290
C***DATE WRITTEN 940830 (YYMMDD)
6291
C***REVISION DATE 951006 (SQRTN, RSQRTN, and scaling of WT added.)
6292
C***REVISION DATE 000628 TSCALE argument added.
6295
C-----------------------------------------------------------------------
6298
C DFNRMK calculates the scaled preconditioned norm of the nonlinear
6299
C function used in the nonlinear iteration for obtaining consistent
6300
C initial conditions. Specifically, DFNRMK calculates the weighted
6301
C root-mean-square norm of the vector (P-inverse)*G(T,Y,YPRIME),
6302
C where P is the preconditioner matrix.
6304
C In addition to the parameters described in the calling program
6305
C DLINSK, the parameters represent
6307
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
6308
C IRIN -- Flag showing whether the current residual vector is
6309
C input in SAVR. 1 means it is, 0 means it is not.
6310
C R -- Array of length NEQ that contains
6311
C (P-inverse)*G(T,Y,YPRIME) on return.
6312
C FNORM -- Scalar containing the weighted norm of R on return.
6313
C-----------------------------------------------------------------------
6316
C RES, DCOPY, DSCAL, PSOL, DDWNRM
6318
C***END PROLOGUE DFNRMK
6321
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6323
DIMENSION Y(*), YPRIME(*), WT(*), SAVR(*), R(*), PWK(*)
6324
DIMENSION WP(*), IWP(*), RPAR(*), IPAR(*)
6325
C-----------------------------------------------------------------------
6326
C Call RES routine if IRIN = 0.
6327
C-----------------------------------------------------------------------
6328
IF (IRIN .EQ. 0) THEN
6330
CALL RES (T, Y, YPRIME, CJ, SAVR, IRES, RPAR, IPAR)
6331
IF (IRES .LT. 0) RETURN
6333
C-----------------------------------------------------------------------
6334
C Apply inverse of left preconditioner to vector R.
6335
C First scale WT array by 1/sqrt(N), and undo scaling afterward.
6336
C-----------------------------------------------------------------------
6337
CALL DCOPY(NEQ, SAVR, 1, R, 1)
6338
CALL DSCAL (NEQ, RSQRTN, WT, 1)
6340
CALL PSOL (NEQ, T, Y, YPRIME, SAVR, PWK, CJ, WT, WP, IWP,
6341
* R, EPLIN, IER, RPAR, IPAR)
6342
CALL DSCAL (NEQ, SQRTN, WT, 1)
6343
IF (IER .NE. 0) RETURN
6344
C-----------------------------------------------------------------------
6345
C Calculate norm of R.
6346
C-----------------------------------------------------------------------
6347
FNORM = DDWNRM (NEQ, R, WT, RPAR, IPAR)
6348
IF (TSCALE .GT. 0.0D0) FNORM = FNORM*TSCALE*ABS(CJ)
6351
C----------------------- END OF SUBROUTINE DFNRMK ----------------------
6353
SUBROUTINE DNEDK(X,Y,YPRIME,NEQ,RES,JACK,PSOL,
6354
* H,WT,JSTART,IDID,RPAR,IPAR,PHI,GAMMA,SAVR,DELTA,E,
6355
* WM,IWM,CJ,CJOLD,CJLAST,S,UROUND,EPLI,SQRTN,RSQRTN,
6356
* EPCON,JCALC,JFLG,KP1,NONNEG,NTYPE,IERNLS)
6358
C***BEGIN PROLOGUE DNEDK
6360
C***DATE WRITTEN 891219 (YYMMDD)
6361
C***REVISION DATE 900926 (YYMMDD)
6362
C***REVISION DATE 940701 (YYMMDD)
6365
C-----------------------------------------------------------------------
6368
C DNEDK solves a nonlinear system of
6369
C algebraic equations of the form
6370
C G(X,Y,YPRIME) = 0 for the unknown Y.
6372
C The method used is a matrix-free Newton scheme.
6374
C The parameters represent
6375
C X -- Independent variable.
6376
C Y -- Solution vector at x.
6377
C YPRIME -- Derivative of solution vector
6378
C after successful step.
6379
C NEQ -- Number of equations to be integrated.
6380
C RES -- External user-supplied subroutine
6381
C to evaluate the residual. See RES description
6382
C in DDASPK prologue.
6383
C JACK -- External user-supplied routine to update
6384
C the preconditioner. (This is optional).
6385
C See JAC description for the case
6386
C INFO(12) = 1 in the DDASPK prologue.
6387
C PSOL -- External user-supplied routine to solve
6388
C a linear system using preconditioning.
6389
C (This is optional). See explanation inside DDASPK.
6390
C H -- Appropriate step size for this step.
6391
C WT -- Vector of weights for error criterion.
6392
C JSTART -- Indicates first call to this routine.
6393
C If JSTART = 0, then this is the first call,
6394
C otherwise it is not.
6395
C IDID -- Completion flag, output by DNEDK.
6396
C See IDID description in DDASPK prologue.
6397
C RPAR,IPAR -- Real and integer arrays used for communication
6398
C between the calling program and external user
6399
C routines. They are not altered within DASPK.
6400
C PHI -- Array of divided differences used by
6401
C DNEDK. The length is NEQ*(K+1), where
6402
C K is the maximum order.
6403
C GAMMA -- Array used to predict Y and YPRIME. The length
6404
C is K+1, where K is the maximum order.
6405
C SAVR -- Work vector for DNEDK of length NEQ.
6406
C DELTA -- Work vector for DNEDK of length NEQ.
6407
C E -- Error accumulation vector for DNEDK of length NEQ.
6408
C WM,IWM -- Real and integer arrays storing
6409
C matrix information for linear system
6410
C solvers, and various other information.
6411
C CJ -- Parameter always proportional to 1/H.
6412
C CJOLD -- Saves the value of CJ as of the last call to DITMD.
6413
C Accounts for changes in CJ needed to
6414
C decide whether to call DITMD.
6415
C CJLAST -- Previous value of CJ.
6416
C S -- A scalar determined by the approximate rate
6417
C of convergence of the Newton iteration and used
6418
C in the convergence test for the Newton iteration.
6420
C If RATE is defined to be an estimate of the
6421
C rate of convergence of the Newton iteration,
6422
C then S = RATE/(1.D0-RATE).
6424
C The closer RATE is to 0., the faster the Newton
6425
C iteration is converging; the closer RATE is to 1.,
6426
C the slower the Newton iteration is converging.
6428
C On the first Newton iteration with an up-dated
6429
C preconditioner S = 100.D0, Thus the initial
6430
C RATE of convergence is approximately 1.
6432
C S is preserved from call to call so that the rate
6433
C estimate from a previous step can be applied to
6435
C UROUND -- Unit roundoff. Not used here.
6436
C EPLI -- convergence test constant.
6437
C See DDASPK prologue for more details.
6438
C SQRTN -- Square root of NEQ.
6439
C RSQRTN -- reciprical of square root of NEQ.
6440
C EPCON -- Tolerance to test for convergence of the Newton
6442
C JCALC -- Flag used to determine when to update
6443
C the Jacobian matrix. In general:
6445
C JCALC = -1 ==> Call the DITMD routine to update
6446
C the Jacobian matrix.
6447
C JCALC = 0 ==> Jacobian matrix is up-to-date.
6448
C JCALC = 1 ==> Jacobian matrix is out-dated,
6449
C but DITMD will not be called unless
6450
C JCALC is set to -1.
6451
C JFLG -- Flag showing whether a Jacobian routine is supplied.
6452
C KP1 -- The current order + 1; updated across calls.
6453
C NONNEG -- Flag to determine nonnegativity constraints.
6454
C NTYPE -- Identification code for the DNEDK routine.
6455
C 1 ==> modified Newton; iterative linear solver.
6456
C 2 ==> modified Newton; user-supplied linear solver.
6457
C IERNLS -- Error flag for nonlinear solver.
6458
C 0 ==> nonlinear solver converged.
6459
C 1 ==> recoverable error inside non-linear solver.
6460
C -1 ==> unrecoverable error inside non-linear solver.
6462
C The following group of variables are passed as arguments to
6463
C the Newton iteration solver. They are explained in greater detail
6465
C TOLNEW, MULDEL, MAXIT, IERNEW
6467
C IERTYP -- Flag which tells whether this subroutine is correct.
6468
C 0 ==> correct subroutine.
6469
C 1 ==> incorrect subroutine.
6471
C-----------------------------------------------------------------------
6473
C RES, JACK, DDWNRM, DNSK
6475
C***END PROLOGUE DNEDK
6478
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
6479
DIMENSION Y(*),YPRIME(*),WT(*)
6480
DIMENSION PHI(NEQ,*),SAVR(*),DELTA(*),E(*)
6481
DIMENSION WM(*),IWM(*)
6482
DIMENSION GAMMA(*),RPAR(*),IPAR(*)
6483
EXTERNAL RES, JACK, PSOL
6485
PARAMETER (LNRE=12, LNJE=13, LLOCWP=29, LLCIWP=30)
6487
SAVE MULDEL, MAXIT, XRATE
6488
DATA MULDEL/0/, MAXIT/4/, XRATE/0.25D0/
6490
C Verify that this is the correct subroutine.
6493
IF (NTYPE .NE. 1) THEN
6498
C If this is the first step, perform initializations.
6500
IF (JSTART .EQ. 0) THEN
6506
C Perform all other initializations.
6512
C Decide whether to update the preconditioner.
6514
IF (JFLG .NE. 0) THEN
6515
TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE)
6517
IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1
6518
IF (CJ .NE. CJLAST) S = 100.D0
6523
C Looping point for updating preconditioner with current stepsize.
6527
C Initialize all error flags to zero.
6534
C Predict the solution and derivative and compute the tolerance
6535
C for the Newton iteration.
6543
320 YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J)
6548
C Call RES to initialize DELTA.
6550
IWM(LNRE)=IWM(LNRE)+1
6551
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
6552
IF (IRES .LT. 0) GO TO 380
6555
C If indicated, update the preconditioner.
6556
C Set JCALC to 0 as an indicator that this has been done.
6558
IF(JCALC .EQ. -1)THEN
6559
IWM(LNJE) = IWM(LNJE) + 1
6561
CALL JACK (RES, IRES, NEQ, X, Y, YPRIME, WT, DELTA, E, H, CJ,
6562
* WM(LWP), IWM(LIWP), IERPJ, RPAR, IPAR)
6565
IF (IRES .LT. 0) GO TO 380
6566
IF (IERPJ .NE. 0) GO TO 380
6569
C Call the nonlinear Newton solver.
6571
CALL DNSK(X,Y,YPRIME,NEQ,RES,PSOL,WT,RPAR,IPAR,SAVR,
6572
* DELTA,E,WM,IWM,CJ,SQRTN,RSQRTN,EPLIN,EPCON,
6573
* S,TEMP1,TOLNEW,MULDEL,MAXIT,IRES,IERSL,IERNEW)
6575
IF (IERNEW .GT. 0 .AND. JCALC .NE. 0) THEN
6577
C The Newton iteration had a recoverable failure with an old
6578
C preconditioner. Retry the step with a new preconditioner.
6584
IF (IERNEW .NE. 0) GO TO 380
6586
C The Newton iteration has converged. If nonnegativity of
6587
C solution is required, set the solution nonnegative, if the
6588
C perturbation to do it is small enough. If the change is too
6589
C large, then consider the corrector iteration to have failed.
6591
IF(NONNEG .EQ. 0) GO TO 390
6593
360 DELTA(I) = MIN(Y(I),0.0D0)
6594
DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
6595
IF(DELNRM .GT. EPCON) GO TO 380
6597
370 E(I) = E(I) - DELTA(I)
6601
C Exits from nonlinear solver.
6602
C No convergence with current preconditioner.
6603
C Compute IERNLS and IDID accordingly.
6606
IF (IRES .LE. -2 .OR. IERSL .LT. 0 .OR. IERTYP .NE. 0) THEN
6608
IF (IRES .LE. -2) IDID = -11
6609
IF (IERSL .LT. 0) IDID = -13
6610
IF (IERTYP .NE. 0) IDID = -15
6613
IF (IRES .EQ. -1) IDID = -10
6614
IF (IERPJ .NE. 0) IDID = -5
6615
IF (IERSL .GT. 0) IDID = -14
6622
C------END OF SUBROUTINE DNEDK------------------------------------------
6624
SUBROUTINE DNSK(X,Y,YPRIME,NEQ,RES,PSOL,WT,RPAR,IPAR,
6625
* SAVR,DELTA,E,WM,IWM,CJ,SQRTN,RSQRTN,EPLIN,EPCON,
6626
* S,CONFAC,TOLNEW,MULDEL,MAXIT,IRES,IERSL,IERNEW)
6628
C***BEGIN PROLOGUE DNSK
6630
C***DATE WRITTEN 891219 (YYMMDD)
6631
C***REVISION DATE 900926 (YYMMDD)
6632
C***REVISION DATE 950126 (YYMMDD)
6633
C***REVISION DATE 000711 (YYMMDD)
6636
C-----------------------------------------------------------------------
6639
C DNSK solves a nonlinear system of
6640
C algebraic equations of the form
6641
C G(X,Y,YPRIME) = 0 for the unknown Y.
6643
C The method used is a modified Newton scheme.
6645
C The parameters represent
6647
C X -- Independent variable.
6648
C Y -- Solution vector.
6649
C YPRIME -- Derivative of solution vector.
6650
C NEQ -- Number of unknowns.
6651
C RES -- External user-supplied subroutine
6652
C to evaluate the residual. See RES description
6653
C in DDASPK prologue.
6654
C PSOL -- External user-supplied routine to solve
6655
C a linear system using preconditioning.
6656
C See explanation inside DDASPK.
6657
C WT -- Vector of weights for error criterion.
6658
C RPAR,IPAR -- Real and integer arrays used for communication
6659
C between the calling program and external user
6660
C routines. They are not altered within DASPK.
6661
C SAVR -- Work vector for DNSK of length NEQ.
6662
C DELTA -- Work vector for DNSK of length NEQ.
6663
C E -- Error accumulation vector for DNSK of length NEQ.
6664
C WM,IWM -- Real and integer arrays storing
6665
C matrix information such as the matrix
6666
C of partial derivatives, permutation
6667
C vector, and various other information.
6668
C CJ -- Parameter always proportional to 1/H (step size).
6669
C SQRTN -- Square root of NEQ.
6670
C RSQRTN -- reciprical of square root of NEQ.
6671
C EPLIN -- Tolerance for linear system solver.
6672
C EPCON -- Tolerance to test for convergence of the Newton
6674
C S -- Used for error convergence tests.
6675
C In the Newton iteration: S = RATE/(1.D0-RATE),
6676
C where RATE is the estimated rate of convergence
6677
C of the Newton iteration.
6679
C The closer RATE is to 0., the faster the Newton
6680
C iteration is converging; the closer RATE is to 1.,
6681
C the slower the Newton iteration is converging.
6683
C The calling routine sends the initial value
6684
C of S to the Newton iteration.
6685
C CONFAC -- A residual scale factor to improve convergence.
6686
C TOLNEW -- Tolerance on the norm of Newton correction in
6687
C alternative Newton convergence test.
6688
C MULDEL -- A flag indicating whether or not to multiply
6690
C 0 ==> do not scale DELTA by CONFAC.
6691
C 1 ==> scale DELTA by CONFAC.
6692
C MAXIT -- Maximum allowed number of Newton iterations.
6693
C IRES -- Error flag returned from RES. See RES description
6694
C in DDASPK prologue. If IRES = -1, then IERNEW
6696
C If IRES < -1, then IERNEW will be set to -1.
6697
C IERSL -- Error flag for linear system solver.
6698
C See IERSL description in subroutine DSLVK.
6699
C If IERSL = 1, then IERNEW will be set to 1.
6700
C If IERSL < 0, then IERNEW will be set to -1.
6701
C IERNEW -- Error flag for Newton iteration.
6702
C 0 ==> Newton iteration converged.
6703
C 1 ==> recoverable error inside Newton iteration.
6704
C -1 ==> unrecoverable error inside Newton iteration.
6705
C-----------------------------------------------------------------------
6708
C RES, DSLVK, DDWNRM
6710
C***END PROLOGUE DNSK
6713
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
6714
DIMENSION Y(*),YPRIME(*),WT(*),DELTA(*),E(*),SAVR(*)
6715
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
6718
PARAMETER (LNNI=19, LNRE=12)
6720
C Initialize Newton counter M and accumulation vector E.
6729
IWM(LNNI) = IWM(LNNI) + 1
6731
C If necessary, multiply residual by convergence factor.
6733
IF (MULDEL .EQ. 1) THEN
6735
320 DELTA(I) = DELTA(I) * CONFAC
6738
C Save residual in SAVR.
6741
340 SAVR(I) = DELTA(I)
6743
C Compute a new iterate. Store the correction in DELTA.
6745
CALL DSLVK (NEQ, Y, X, YPRIME, SAVR, DELTA, WT, WM, IWM,
6746
* RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK,
6748
IF (IRES .NE. 0 .OR. IERSL .NE. 0) GO TO 380
6750
C Update Y, E, and YPRIME.
6753
Y(I) = Y(I) - DELTA(I)
6754
E(I) = E(I) - DELTA(I)
6755
360 YPRIME(I) = YPRIME(I) - CJ*DELTA(I)
6757
C Test for convergence of the iteration.
6759
DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
6762
IF (DELNRM .LE. TOLNEW) GO TO 370
6764
RATE = (DELNRM/OLDNRM)**(1.0D0/M)
6765
IF (RATE .GT. 0.9D0) GO TO 380
6766
S = RATE/(1.0D0 - RATE)
6768
IF (S*DELNRM .LE. EPCON) GO TO 370
6770
C The corrector has not yet converged. Update M and test whether
6771
C the maximum number of iterations have been tried.
6774
IF (M .GE. MAXIT) GO TO 380
6776
C Evaluate the residual, and go back to do another iteration.
6778
IWM(LNRE) = IWM(LNRE) + 1
6779
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
6780
IF (IRES .LT. 0) GO TO 380
6783
C The iteration has converged.
6787
C The iteration has not converged. Set IERNEW appropriately.
6790
IF (IRES .LE. -2 .OR. IERSL .LT. 0) THEN
6798
C------END OF SUBROUTINE DNSK-------------------------------------------
6800
SUBROUTINE DSLVK (NEQ, Y, TN, YPRIME, SAVR, X, EWT, WM, IWM,
6801
* RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK,
6804
C***BEGIN PROLOGUE DSLVK
6806
C***DATE WRITTEN 890101 (YYMMDD)
6807
C***REVISION DATE 900926 (YYMMDD)
6808
C***REVISION DATE 940928 Removed MNEWT and added RHOK in call list.
6811
C-----------------------------------------------------------------------
6814
C DSLVK uses a restart algorithm and interfaces to DSPIGM for
6815
C the solution of the linear system arising from a Newton iteration.
6817
C In addition to variables described elsewhere,
6818
C communication with DSLVK uses the following variables..
6819
C WM = Real work space containing data for the algorithm
6820
C (Krylov basis vectors, Hessenberg matrix, etc.).
6821
C IWM = Integer work space containing data for the algorithm.
6822
C X = The right-hand side vector on input, and the solution vector
6823
C on output, of length NEQ.
6824
C IRES = Error flag from RES.
6825
C IERSL = Output flag ..
6826
C IERSL = 0 means no trouble occurred (or user RES routine
6827
C returned IRES < 0)
6828
C IERSL = 1 means the iterative method failed to converge
6829
C (DSPIGM returned IFLAG > 0.)
6830
C IERSL = -1 means there was a nonrecoverable error in the
6831
C iterative solver, and an error exit will occur.
6832
C-----------------------------------------------------------------------
6834
C DSCAL, DCOPY, DSPIGM
6836
C***END PROLOGUE DSLVK
6838
INTEGER NEQ, IWM, IRES, IERSL, IPAR
6839
DOUBLE PRECISION Y, TN, YPRIME, SAVR, X, EWT, WM, CJ, EPLIN,
6840
1 SQRTN, RSQRTN, RHOK, RPAR
6841
DIMENSION Y(*), YPRIME(*), SAVR(*), X(*), EWT(*),
6842
1 WM(*), IWM(*), RPAR(*), IPAR(*)
6844
INTEGER IFLAG, IRST, NRSTS, NRMAX, LR, LDL, LHES, LGMR, LQ, LV,
6845
1 LWK, LZ, MAXLP1, NPSL
6846
INTEGER NLI, NPS, NCFL, NRE, MAXL, KMP, MITER
6849
PARAMETER (LNRE=12, LNCFL=16, LNLI=20, LNPS=21)
6850
PARAMETER (LLOCWP=29, LLCIWP=30)
6851
PARAMETER (LMITER=23, LMAXL=24, LKMP=25, LNRMAX=26)
6853
C-----------------------------------------------------------------------
6854
C IRST is set to 1, to indicate restarting is in effect.
6855
C NRMAX is the maximum number of restarts.
6856
C-----------------------------------------------------------------------
6871
C-----------------------------------------------------------------------
6872
C Use a restarting strategy to solve the linear system
6873
C P*X = -F. Parse the work vector, and perform initializations.
6874
C Note that zero is the initial guess for X.
6875
C-----------------------------------------------------------------------
6880
LQ = LHES + MAXL*MAXLP1
6882
LDL = LWK + MIN0(1,MAXL-KMP)*NEQ
6884
CALL DSCAL (NEQ, RSQRTN, EWT, 1)
6885
CALL DCOPY (NEQ, X, 1, WM(LR), 1)
6888
C-----------------------------------------------------------------------
6889
C Top of loop for the restart algorithm. Initial pass approximates
6890
C X and sets up a transformed system to perform subsequent restarts
6891
C to update X. NRSTS is initialized to -1, because restarting
6892
C does not occur until after the first pass.
6893
C Update NRSTS; conditionally copy DL to R; call the DSPIGM
6894
C algorithm to solve A*Z = R; updated counters; update X with
6895
C the residual solution.
6896
C Note: if convergence is not achieved after NRMAX restarts,
6897
C then the linear solver is considered to have failed.
6898
C-----------------------------------------------------------------------
6902
IF (NRSTS .GT. 0) CALL DCOPY (NEQ, WM(LDL), 1, WM(LR),1)
6903
CALL DSPIGM (NEQ, TN, Y, YPRIME, SAVR, WM(LR), EWT, MAXL, MAXLP1,
6904
1 KMP, EPLIN, CJ, RES, IRES, NRES, PSOL, NPSL, WM(LZ), WM(LV),
6905
2 WM(LHES), WM(LQ), LGMR, WM(LWP), IWM(LIWP), WM(LWK),
6906
3 WM(LDL), RHOK, IFLAG, IRST, NRSTS, RPAR, IPAR)
6911
120 X(I) = X(I) + WM(LZ+I-1)
6912
IF ((IFLAG .EQ. 1) .AND. (NRSTS .LT. NRMAX) .AND. (IRES .EQ. 0))
6914
C-----------------------------------------------------------------------
6915
C The restart scheme is finished. Test IRES and IFLAG to see if
6916
C convergence was not achieved, and set flags accordingly.
6917
C-----------------------------------------------------------------------
6918
IF (IRES .LT. 0) THEN
6920
ELSE IF (IFLAG .NE. 0) THEN
6922
IF (IFLAG .GT. 0) IERSL = 1
6923
IF (IFLAG .LT. 0) IERSL = -1
6925
C-----------------------------------------------------------------------
6926
C Update IWM with counters, rescale EWT, and return.
6927
C-----------------------------------------------------------------------
6932
CALL DSCAL (NEQ, SQRTN, EWT, 1)
6935
C------END OF SUBROUTINE DSLVK------------------------------------------
6937
SUBROUTINE DSPIGM (NEQ, TN, Y, YPRIME, SAVR, R, WGHT, MAXL,
6938
* MAXLP1, KMP, EPLIN, CJ, RES, IRES, NRE, PSOL, NPSL, Z, V,
6939
* HES, Q, LGMR, WP, IWP, WK, DL, RHOK, IFLAG, IRST, NRSTS,
6942
C***BEGIN PROLOGUE DSPIGM
6943
C***DATE WRITTEN 890101 (YYMMDD)
6944
C***REVISION DATE 900926 (YYMMDD)
6945
C***REVISION DATE 940927 Removed MNEWT and added RHOK in call list.
6948
C-----------------------------------------------------------------------
6951
C This routine solves the linear system A * Z = R using a scaled
6952
C preconditioned version of the generalized minimum residual method.
6953
C An initial guess of Z = 0 is assumed.
6957
C NEQ = Problem size, passed to PSOL.
6959
C TN = Current Value of T.
6961
C Y = Array Containing current dependent variable vector.
6963
C YPRIME = Array Containing current first derivative of Y.
6965
C SAVR = Array containing current value of G(T,Y,YPRIME).
6967
C R = The right hand side of the system A*Z = R.
6968
C R is also used as work space when computing
6969
C the final approximation and will therefore be
6971
C (R is the same as V(*,MAXL+1) in the call to DSPIGM.)
6973
C WGHT = The vector of length NEQ containing the nonzero
6974
C elements of the diagonal scaling matrix.
6976
C MAXL = The maximum allowable order of the matrix H.
6978
C MAXLP1 = MAXL + 1, used for dynamic dimensioning of HES.
6980
C KMP = The number of previous vectors the new vector, VNEW,
6981
C must be made orthogonal to. (KMP .LE. MAXL.)
6983
C EPLIN = Tolerance on residuals R-A*Z in weighted rms norm.
6985
C CJ = Scalar proportional to current value of
6988
C WK = Real work array used by routine DATV and PSOL.
6990
C DL = Real work array used for calculation of the residual
6991
C norm RHO when the method is incomplete (KMP.LT.MAXL)
6992
C and/or when using restarting.
6994
C WP = Real work array used by preconditioner PSOL.
6996
C IWP = Integer work array used by preconditioner PSOL.
6998
C IRST = Method flag indicating if restarting is being
6999
C performed. IRST .GT. 0 means restarting is active,
7000
C while IRST = 0 means restarting is not being used.
7002
C NRSTS = Counter for the number of restarts on the current
7003
C call to DSPIGM. If NRSTS .GT. 0, then the residual
7004
C R is already scaled, and so scaling of R is not
7010
C Z = The final computed approximation to the solution
7011
C of the system A*Z = R.
7013
C LGMR = The number of iterations performed and
7014
C the current order of the upper Hessenberg
7017
C NRE = The number of calls to RES (i.e. DATV)
7019
C NPSL = The number of calls to PSOL.
7021
C V = The neq by (LGMR+1) array containing the LGMR
7022
C orthogonal vectors V(*,1) to V(*,LGMR).
7024
C HES = The upper triangular factor of the QR decomposition
7025
C of the (LGMR+1) by LGMR upper Hessenberg matrix whose
7026
C entries are the scaled inner-products of A*V(*,I)
7029
C Q = Real array of length 2*MAXL containing the components
7030
C of the givens rotations used in the QR decomposition
7031
C of HES. It is loaded in DHEQR and used in DHELS.
7033
C IRES = Error flag from RES.
7035
C DL = Scaled preconditioned residual,
7036
C (D-inverse)*(P-inverse)*(R-A*Z). Only loaded when
7037
C performing restarts of the Krylov iteration.
7039
C RHOK = Weighted norm of final preconditioned residual.
7041
C IFLAG = Integer error flag..
7042
C 0 Means convergence in LGMR iterations, LGMR.LE.MAXL.
7043
C 1 Means the convergence test did not pass in MAXL
7044
C iterations, but the new residual norm (RHO) is
7045
C .LT. the old residual norm (RNRM), and so Z is
7047
C 2 Means the convergence test did not pass in MAXL
7048
C iterations, new residual norm (RHO) .GE. old residual
7049
C norm (RNRM), and the initial guess, Z = 0, is
7051
C 3 Means there was a recoverable error in PSOL
7052
C caused by the preconditioner being out of date.
7053
C -1 Means there was an unrecoverable error in PSOL.
7055
C-----------------------------------------------------------------------
7057
C PSOL, DNRM2, DSCAL, DATV, DORTH, DHEQR, DCOPY, DHELS, DAXPY
7059
C***END PROLOGUE DSPIGM
7061
INTEGER NEQ,MAXL,MAXLP1,KMP,IRES,NRE,NPSL,LGMR,IWP,
7062
1 IFLAG,IRST,NRSTS,IPAR
7063
DOUBLE PRECISION TN,Y,YPRIME,SAVR,R,WGHT,EPLIN,CJ,Z,V,HES,Q,WP,WK,
7065
DIMENSION Y(*), YPRIME(*), SAVR(*), R(*), WGHT(*), Z(*),
7066
1 V(NEQ,*), HES(MAXLP1,*), Q(*), WP(*), IWP(*), WK(*), DL(*),
7068
INTEGER I, IER, INFO, IP1, I2, J, K, LL, LLP1
7069
DOUBLE PRECISION RNRM,C,DLNRM,PROD,RHO,S,SNORMW,DNRM2,TEM
7077
C-----------------------------------------------------------------------
7078
C The initial guess for Z is 0. The initial residual is therefore
7079
C the vector R. Initialize Z to 0.
7080
C-----------------------------------------------------------------------
7083
C-----------------------------------------------------------------------
7084
C Apply inverse of left preconditioner to vector R if NRSTS .EQ. 0.
7085
C Form V(*,1), the scaled preconditioned right hand side.
7086
C-----------------------------------------------------------------------
7087
IF (NRSTS .EQ. 0) THEN
7088
CALL PSOL (NEQ, TN, Y, YPRIME, SAVR, WK, CJ, WGHT, WP, IWP,
7089
1 R, EPLIN, IER, RPAR, IPAR)
7091
IF (IER .NE. 0) GO TO 300
7093
30 V(I,1) = R(I)*WGHT(I)
7098
C-----------------------------------------------------------------------
7099
C Calculate norm of scaled vector V(*,1) and normalize it
7100
C If, however, the norm of V(*,1) (i.e. the norm of the preconditioned
7101
C residual) is .le. EPLIN, then return with Z=0.
7102
C-----------------------------------------------------------------------
7103
RNRM = DNRM2 (NEQ, V, 1)
7104
IF (RNRM .LE. EPLIN) THEN
7109
CALL DSCAL (NEQ, TEM, V(1,1), 1)
7110
C-----------------------------------------------------------------------
7111
C Zero out the HES array.
7112
C-----------------------------------------------------------------------
7117
C-----------------------------------------------------------------------
7118
C Main loop to compute the vectors V(*,2) to V(*,MAXL).
7119
C The running product PROD is needed for the convergence test.
7120
C-----------------------------------------------------------------------
7124
C-----------------------------------------------------------------------
7125
C Call routine DATV to compute VNEW = ABAR*V(LL), where ABAR is
7126
C the matrix A with scaling and inverse preconditioner factors applied.
7127
C Call routine DORTH to orthogonalize the new vector VNEW = V(*,LL+1).
7128
C call routine DHEQR to update the factors of HES.
7129
C-----------------------------------------------------------------------
7130
CALL DATV (NEQ, Y, TN, YPRIME, SAVR, V(1,LL), WGHT, Z,
7131
1 RES, IRES, PSOL, V(1,LL+1), WK, WP, IWP, CJ, EPLIN,
7132
1 IER, NRE, NPSL, RPAR, IPAR)
7133
IF (IRES .LT. 0) RETURN
7134
IF (IER .NE. 0) GO TO 300
7135
CALL DORTH (V(1,LL+1), V, HES, NEQ, LL, MAXLP1, KMP, SNORMW)
7136
HES(LL+1,LL) = SNORMW
7137
CALL DHEQR (HES, MAXLP1, LL, Q, INFO, LL)
7138
IF (INFO .EQ. LL) GO TO 120
7139
C-----------------------------------------------------------------------
7140
C Update RHO, the estimate of the norm of the residual R - A*ZL.
7141
C If KMP .LT. MAXL, then the vectors V(*,1),...,V(*,LL+1) are not
7142
C necessarily orthogonal for LL .GT. KMP. The vector DL must then
7143
C be computed, and its norm used in the calculation of RHO.
7144
C-----------------------------------------------------------------------
7146
RHO = ABS(PROD*RNRM)
7147
IF ((LL.GT.KMP) .AND. (KMP.LT.MAXL)) THEN
7148
IF (LL .EQ. KMP+1) THEN
7149
CALL DCOPY (NEQ, V(1,1), 1, DL, 1)
7156
70 DL(K) = S*DL(K) + C*V(K,IP1)
7160
C = Q(2*LL-1)/SNORMW
7163
80 DL(K) = S*DL(K) + C*V(K,LLP1)
7164
DLNRM = DNRM2 (NEQ, DL, 1)
7167
C-----------------------------------------------------------------------
7168
C Test for convergence. If passed, compute approximation ZL.
7169
C If failed and LL .LT. MAXL, then continue iterating.
7170
C-----------------------------------------------------------------------
7171
IF (RHO .LE. EPLIN) GO TO 200
7172
IF (LL .EQ. MAXL) GO TO 100
7173
C-----------------------------------------------------------------------
7174
C Rescale so that the norm of V(1,LL+1) is one.
7175
C-----------------------------------------------------------------------
7177
CALL DSCAL (NEQ, TEM, V(1,LL+1), 1)
7180
IF (RHO .LT. RNRM) GO TO 150
7187
C-----------------------------------------------------------------------
7188
C The tolerance was not met, but the residual norm was reduced.
7189
C If performing restarting (IRST .gt. 0) calculate the residual vector
7190
C RL and store it in the DL array. If the incomplete version is
7191
C being used (KMP .lt. MAXL) then DL has already been calculated.
7192
C-----------------------------------------------------------------------
7193
IF (IRST .GT. 0) THEN
7194
IF (KMP .EQ. MAXL) THEN
7196
C Calculate DL from the V(I)'s.
7198
CALL DCOPY (NEQ, V(1,1), 1, DL, 1)
7206
170 DL(K) = S*DL(K) + C*V(K,IP1)
7209
C = Q(2*MAXL-1)/SNORMW
7211
180 DL(K) = S*DL(K) + C*V(K,MAXLP1)
7214
C Scale DL by RNRM*PROD to obtain the residual RL.
7217
CALL DSCAL(NEQ, TEM, DL, 1)
7219
C-----------------------------------------------------------------------
7220
C Compute the approximation ZL to the solution.
7221
C Since the vector Z was used as work space, and the initial guess
7222
C of the Newton correction is zero, Z must be reset to zero.
7223
C-----------------------------------------------------------------------
7230
CALL DHELS (HES, MAXLP1, LL, Q, R)
7234
CALL DAXPY (NEQ, R(I), V(1,I), 1, Z, 1)
7237
240 Z(I) = Z(I)/WGHT(I)
7238
C Load RHO into RHOK.
7241
C-----------------------------------------------------------------------
7242
C This block handles error returns forced by routine PSOL.
7243
C-----------------------------------------------------------------------
7245
IF (IER .LT. 0) IFLAG = -1
7246
IF (IER .GT. 0) IFLAG = 3
7250
C------END OF SUBROUTINE DSPIGM-----------------------------------------
7252
SUBROUTINE DATV (NEQ, Y, TN, YPRIME, SAVR, V, WGHT, YPTEM, RES,
7253
* IRES, PSOL, Z, VTEM, WP, IWP, CJ, EPLIN, IER, NRE, NPSL,
7256
C***BEGIN PROLOGUE DATV
7257
C***DATE WRITTEN 890101 (YYMMDD)
7258
C***REVISION DATE 900926 (YYMMDD)
7261
C-----------------------------------------------------------------------
7264
C This routine computes the product
7266
C Z = (D-inverse)*(P-inverse)*(dF/dY)*(D*V),
7268
C where F(Y) = G(T, Y, CJ*(Y-A)), CJ is a scalar proportional to 1/H,
7269
C and A involves the past history of Y. The quantity CJ*(Y-A) is
7270
C an approximation to the first derivative of Y and is stored
7271
C in the array YPRIME. Note that dF/dY = dG/dY + CJ*dG/dYPRIME.
7273
C D is a diagonal scaling matrix, and P is the left preconditioning
7274
C matrix. V is assumed to have L2 norm equal to 1.
7275
C The product is stored in Z and is computed by means of a
7276
C difference quotient, a call to RES, and one call to PSOL.
7280
C NEQ = Problem size, passed to RES and PSOL.
7282
C Y = Array containing current dependent variable vector.
7284
C YPRIME = Array containing current first derivative of y.
7286
C SAVR = Array containing current value of G(T,Y,YPRIME).
7288
C V = Real array of length NEQ (can be the same array as Z).
7290
C WGHT = Array of length NEQ containing scale factors.
7291
C 1/WGHT(I) are the diagonal elements of the matrix D.
7293
C YPTEM = Work array of length NEQ.
7295
C VTEM = Work array of length NEQ used to store the
7296
C unscaled version of V.
7298
C WP = Real work array used by preconditioner PSOL.
7300
C IWP = Integer work array used by preconditioner PSOL.
7302
C CJ = Scalar proportional to current value of
7308
C Z = Array of length NEQ containing desired scaled
7309
C matrix-vector product.
7311
C IRES = Error flag from RES.
7313
C IER = Error flag from PSOL.
7315
C NRE = The number of calls to RES.
7317
C NPSL = The number of calls to PSOL.
7319
C-----------------------------------------------------------------------
7323
C***END PROLOGUE DATV
7325
INTEGER NEQ, IRES, IWP, IER, NRE, NPSL, IPAR
7326
DOUBLE PRECISION Y, TN, YPRIME, SAVR, V, WGHT, YPTEM, Z, VTEM,
7328
DIMENSION Y(*), YPRIME(*), SAVR(*), V(*), WGHT(*), YPTEM(*),
7329
1 Z(*), VTEM(*), WP(*), IWP(*), RPAR(*), IPAR(*)
7331
DOUBLE PRECISION EPLIN
7335
C-----------------------------------------------------------------------
7337
C-----------------------------------------------------------------------
7339
10 VTEM(I) = V(I)/WGHT(I)
7341
C-----------------------------------------------------------------------
7342
C Store Y in Z and increment Z by VTEM.
7343
C Store YPRIME in YPTEM and increment YPTEM by VTEM*CJ.
7344
C-----------------------------------------------------------------------
7346
YPTEM(I) = YPRIME(I) + VTEM(I)*CJ
7347
20 Z(I) = Y(I) + VTEM(I)
7348
C-----------------------------------------------------------------------
7349
C Call RES with incremented Y, YPRIME arguments
7350
C stored in Z, YPTEM. VTEM is overwritten with new residual.
7351
C-----------------------------------------------------------------------
7353
CALL RES(TN,Z,YPTEM,CJ,VTEM,IRES,RPAR,IPAR)
7355
IF (IRES .LT. 0) RETURN
7356
C-----------------------------------------------------------------------
7357
C Set Z = (dF/dY) * VBAR using difference quotient.
7358
C (VBAR is old value of VTEM before calling RES)
7359
C-----------------------------------------------------------------------
7361
70 Z(I) = VTEM(I) - SAVR(I)
7362
C-----------------------------------------------------------------------
7363
C Apply inverse of left preconditioner to Z.
7364
C-----------------------------------------------------------------------
7365
CALL PSOL (NEQ, TN, Y, YPRIME, SAVR, YPTEM, CJ, WGHT, WP, IWP,
7366
1 Z, EPLIN, IER, RPAR, IPAR)
7368
IF (IER .NE. 0) RETURN
7369
C-----------------------------------------------------------------------
7370
C Apply D-inverse to Z and return.
7371
C-----------------------------------------------------------------------
7373
90 Z(I) = Z(I)*WGHT(I)
7376
C------END OF SUBROUTINE DATV-------------------------------------------
7378
SUBROUTINE DORTH (VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
7380
C***BEGIN PROLOGUE DORTH
7381
C***DATE WRITTEN 890101 (YYMMDD)
7382
C***REVISION DATE 900926 (YYMMDD)
7385
C-----------------------------------------------------------------------
7388
C This routine orthogonalizes the vector VNEW against the previous
7389
C KMP vectors in the V array. It uses a modified Gram-Schmidt
7390
C orthogonalization procedure with conditional reorthogonalization.
7394
C VNEW = The vector of length N containing a scaled product
7395
C OF The Jacobian and the vector V(*,LL).
7397
C V = The N x LL array containing the previous LL
7398
C orthogonal vectors V(*,1) to V(*,LL).
7400
C HES = An LL x LL upper Hessenberg matrix containing,
7401
C in HES(I,K), K.LT.LL, scaled inner products of
7402
C A*V(*,K) and V(*,I).
7404
C LDHES = The leading dimension of the HES array.
7406
C N = The order of the matrix A, and the length of VNEW.
7408
C LL = The current order of the matrix HES.
7410
C KMP = The number of previous vectors the new vector VNEW
7411
C must be made orthogonal to (KMP .LE. MAXL).
7416
C VNEW = The new vector orthogonal to V(*,I0),
7417
C where I0 = MAX(1, LL-KMP+1).
7419
C HES = Upper Hessenberg matrix with column LL filled in with
7420
C scaled inner products of A*V(*,LL) and V(*,I).
7422
C SNORMW = L-2 norm of VNEW.
7424
C-----------------------------------------------------------------------
7426
C DDOT, DNRM2, DAXPY
7428
C***END PROLOGUE DORTH
7430
INTEGER N, LL, LDHES, KMP
7431
DOUBLE PRECISION VNEW, V, HES, SNORMW
7432
DIMENSION VNEW(*), V(N,*), HES(LDHES,*)
7434
DOUBLE PRECISION ARG, DDOT, DNRM2, SUMDSQ, TEM, VNRM
7436
C-----------------------------------------------------------------------
7437
C Get norm of unaltered VNEW for later use.
7438
C-----------------------------------------------------------------------
7439
VNRM = DNRM2 (N, VNEW, 1)
7440
C-----------------------------------------------------------------------
7441
C Do Modified Gram-Schmidt on VNEW = A*V(LL).
7442
C Scaled inner products give new column of HES.
7443
C Projections of earlier vectors are subtracted from VNEW.
7444
C-----------------------------------------------------------------------
7445
I0 = MAX0(1,LL-KMP+1)
7447
HES(I,LL) = DDOT (N, V(1,I), 1, VNEW, 1)
7449
CALL DAXPY (N, TEM, V(1,I), 1, VNEW, 1)
7451
C-----------------------------------------------------------------------
7452
C Compute SNORMW = norm of VNEW.
7453
C If VNEW is small compared to its input value (in norm), then
7454
C Reorthogonalize VNEW to V(*,1) through V(*,LL).
7455
C Correct if relative correction exceeds 1000*(unit roundoff).
7456
C Finally, correct SNORMW using the dot products involved.
7457
C-----------------------------------------------------------------------
7458
SNORMW = DNRM2 (N, VNEW, 1)
7459
IF (VNRM + 0.001D0*SNORMW .NE. VNRM) RETURN
7462
TEM = -DDOT (N, V(1,I), 1, VNEW, 1)
7463
IF (HES(I,LL) + 0.001D0*TEM .EQ. HES(I,LL)) GO TO 30
7464
HES(I,LL) = HES(I,LL) - TEM
7465
CALL DAXPY (N, TEM, V(1,I), 1, VNEW, 1)
7466
SUMDSQ = SUMDSQ + TEM**2
7468
IF (SUMDSQ .EQ. 0.0D0) RETURN
7469
ARG = MAX(0.0D0,SNORMW**2 - SUMDSQ)
7473
C------END OF SUBROUTINE DORTH------------------------------------------
7475
SUBROUTINE DHEQR (A, LDA, N, Q, INFO, IJOB)
7477
C***BEGIN PROLOGUE DHEQR
7478
C***DATE WRITTEN 890101 (YYMMDD)
7479
C***REVISION DATE 900926 (YYMMDD)
7481
C-----------------------------------------------------------------------
7484
C This routine performs a QR decomposition of an upper
7485
C Hessenberg matrix A. There are two options available:
7487
C (1) performing a fresh decomposition
7488
C (2) updating the QR factors by adding a row and A
7489
C column to the matrix A.
7491
C DHEQR decomposes an upper Hessenberg matrix by using Givens
7496
C A DOUBLE PRECISION(LDA, N)
7497
C The matrix to be decomposed.
7500
C The leading dimension of the array A.
7503
C A is an (N+1) by N Hessenberg matrix.
7506
C = 1 Means that a fresh decomposition of the
7507
C matrix A is desired.
7508
C .GE. 2 Means that the current decomposition of A
7509
C will be updated by the addition of a row
7513
C A The upper triangular matrix R.
7514
C The factorization can be written Q*A = R, where
7515
C Q is a product of Givens rotations and R is upper
7518
C Q DOUBLE PRECISION(2*N)
7519
C The factors C and S of each Givens rotation used
7524
C = K If A(K,K) .EQ. 0.0. This is not an error
7525
C condition for this subroutine, but it does
7526
C indicate that DHELS will divide by zero
7529
C Modification of LINPACK.
7530
C Peter Brown, Lawrence Livermore Natl. Lab.
7532
C-----------------------------------------------------------------------
7533
C***ROUTINES CALLED (NONE)
7535
C***END PROLOGUE DHEQR
7537
INTEGER LDA, N, INFO, IJOB
7538
DOUBLE PRECISION A(LDA,*), Q(*)
7539
INTEGER I, IQ, J, K, KM1, KP1, NM1
7540
DOUBLE PRECISION C, S, T, T1, T2
7542
IF (IJOB .GT. 1) GO TO 70
7543
C-----------------------------------------------------------------------
7544
C A new factorization is desired.
7545
C-----------------------------------------------------------------------
7547
C QR decomposition without pivoting.
7554
C Compute Kth column of R.
7555
C First, multiply the Kth column of A by the previous
7556
C K-1 Givens rotations.
7558
IF (KM1 .LT. 1) GO TO 20
7565
A(J,K) = C*T1 - S*T2
7566
A(J+1,K) = S*T1 + C*T2
7569
C Compute Givens components C and S.
7575
IF (T2 .NE. 0.0D0) GO TO 30
7580
IF (ABS(T2) .LT. ABS(T1)) GO TO 40
7582
S = -1.0D0/SQRT(1.0D0+T*T)
7587
C = 1.0D0/SQRT(1.0D0+T*T)
7592
A(K,K) = C*T1 - S*T2
7593
IF (A(K,K) .EQ. 0.0D0) INFO = K
7596
C-----------------------------------------------------------------------
7597
C The old factorization of A will be updated. A row and a column
7598
C has been added to the matrix A.
7599
C N by N-1 is now the old size of the matrix.
7600
C-----------------------------------------------------------------------
7603
C-----------------------------------------------------------------------
7604
C Multiply the new column by the N previous Givens rotations.
7605
C-----------------------------------------------------------------------
7612
A(K,N) = C*T1 - S*T2
7613
A(K+1,N) = S*T1 + C*T2
7615
C-----------------------------------------------------------------------
7616
C Complete update of decomposition by forming last Givens rotation,
7617
C and multiplying it times the column vector (A(N,N),A(NP1,N)).
7618
C-----------------------------------------------------------------------
7622
IF (T2 .NE. 0.0D0) GO TO 110
7627
IF (ABS(T2) .LT. ABS(T1)) GO TO 120
7629
S = -1.0D0/SQRT(1.0D0+T*T)
7634
C = 1.0D0/SQRT(1.0D0+T*T)
7640
A(N,N) = C*T1 - S*T2
7641
IF (A(N,N) .EQ. 0.0D0) INFO = N
7644
C------END OF SUBROUTINE DHEQR------------------------------------------
7646
SUBROUTINE DHELS (A, LDA, N, Q, B)
7648
C***BEGIN PROLOGUE DHELS
7649
C***DATE WRITTEN 890101 (YYMMDD)
7650
C***REVISION DATE 900926 (YYMMDD)
7653
C-----------------------------------------------------------------------
7656
C This is similar to the LINPACK routine DGESL except that
7657
C A is an upper Hessenberg matrix.
7659
C DHELS solves the least squares problem
7663
C using the factors computed by DHEQR.
7667
C A DOUBLE PRECISION (LDA, N)
7668
C The output from DHEQR which contains the upper
7669
C triangular factor R in the QR decomposition of A.
7672
C The leading dimension of the array A .
7675
C A is originally an (N+1) by N matrix.
7677
C Q DOUBLE PRECISION(2*N)
7678
C The coefficients of the N givens rotations
7679
C used in the QR factorization of A.
7681
C B DOUBLE PRECISION(N+1)
7682
C The right hand side vector.
7687
C B The solution vector X.
7690
C Modification of LINPACK.
7691
C Peter Brown, Lawrence Livermore Natl. Lab.
7693
C-----------------------------------------------------------------------
7697
C***END PROLOGUE DHELS
7700
DOUBLE PRECISION A(LDA,*), B(*), Q(*)
7701
INTEGER IQ, K, KB, KP1
7702
DOUBLE PRECISION C, S, T, T1, T2
7704
C Minimize (B-A*X,B-A*X).
7715
B(KP1) = S*T1 + C*T2
7718
C Now solve R*X = Q*B.
7724
CALL DAXPY (K-1, T, A(1,K), 1, B(1), 1)
7728
C------END OF SUBROUTINE DHELS------------------------------------------
added
removed
routines/integ/ddaskr.f
