fix import/export of ierode common on Windows
[scilab.git] / scilab / modules / differential_equations / src / fortran / ddassl.f
1       SUBROUTINE DDAINI (X, Y, YPRIME, NEQ, RES, JAC, H, WT, IDID, RPAR,
2      +   IPAR, PHI, DELTA, E, WM, IWM, HMIN, UROUND, NONNEG, NTEMP)
3
4       include 'stack.h'
5 C***BEGIN PROLOGUE  DDAINI
6 C***SUBSIDIARY
7 C***PURPOSE  Initialization routine for DDASSL.
8 C***LIBRARY   SLATEC (DASSL)
9 C***TYPE      DOUBLE PRECISION (SDAINI-S, DDAINI-D)
10 C***AUTHOR  PETZOLD, LINDA R., (LLNL)
11 C***DESCRIPTION
12 C-----------------------------------------------------------------
13 C     DDAINI TAKES ONE STEP OF SIZE H OR SMALLER
14 C     WITH THE BACKWARD EULER METHOD, TO
15 C     FIND YPRIME.  X AND Y ARE UPDATED TO BE CONSISTENT WITH THE
16 C     NEW STEP.  A MODIFIED DAMPED NEWTON ITERATION IS USED TO
17 C     SOLVE THE CORRECTOR ITERATION.
18 C
19 C     THE INITIAL GUESS FOR YPRIME IS USED IN THE
20 C     PREDICTION, AND IN FORMING THE ITERATION
21 C     MATRIX, BUT IS NOT INVOLVED IN THE
22 C     ERROR TEST. THIS MAY HAVE TROUBLE
23 C     CONVERGING IF THE INITIAL GUESS IS NO
24 C     GOOD, OR IF G(X,Y,YPRIME) DEPENDS
25 C     NONLINEARLY ON YPRIME.
26 C
27 C     THE PARAMETERS REPRESENT:
28 C     X --         INDEPENDENT VARIABLE
29 C     Y --         SOLUTION VECTOR AT X
30 C     YPRIME --    DERIVATIVE OF SOLUTION VECTOR
31 C     NEQ --       NUMBER OF EQUATIONS
32 C     H --         STEPSIZE. IMDER MAY USE A STEPSIZE
33 C                  SMALLER THAN H.
34 C     WT --        VECTOR OF WEIGHTS FOR ERROR
35 C                  CRITERION
36 C     IDID --      COMPLETION CODE WITH THE FOLLOWING MEANINGS
37 C                  IDID= 1 -- YPRIME WAS FOUND SUCCESSFULLY
38 C                  IDID=-12 -- DDAINI FAILED TO FIND YPRIME
39 C     RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS
40 C                  THAT ARE NOT ALTERED BY DDAINI
41 C     PHI --       WORK SPACE FOR DDAINI
42 C     DELTA,E --   WORK SPACE FOR DDAINI
43 C     WM,IWM --    REAL AND INTEGER ARRAYS STORING
44 C                  MATRIX INFORMATION
45 C
46 C-----------------------------------------------------------------
47 C***ROUTINES CALLED  DDAJAC, DDANRM, DDASLV
48 C***REVISION HISTORY  (YYMMDD)
49 C   830315  DATE WRITTEN
50 C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
51 C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
52 C   901026  Added explicit declarations for all variables and minor
53 C           cosmetic changes to prologue.  (FNF)
54 C   901030  Minor corrections to declarations.  (FNF)
55 C***END PROLOGUE  DDAINI
56 C
57       INTEGER  NEQ, IDID, IPAR(*), IWM(*), NONNEG, NTEMP
58       DOUBLE PRECISION
59      *   X, Y(*), YPRIME(*), H, WT(*), RPAR(*), PHI(NEQ,*), DELTA(*),
60      *   E(*), WM(*), HMIN, UROUND
61       EXTERNAL  RES, JAC
62 C
63       EXTERNAL  DDAJAC, DDANRM, DDASLV
64       DOUBLE PRECISION  DDANRM
65 C
66       INTEGER  I, IER, IRES, JCALC, LNJE, LNRE, M, MAXIT, MJAC, NCF,
67      *   NEF, NSF
68       DOUBLE PRECISION
69      *   CJ, DAMP, DELNRM, IERR, OLDNRM, R, RATE, S, XOLD, YNORM
70       LOGICAL  CONVGD
71 C
72       PARAMETER (LNRE=12)
73       PARAMETER (LNJE=13)
74 C
75       DATA MAXIT/10/,MJAC/5/
76       DATA DAMP/0.75D0/
77 C
78 C
79 C---------------------------------------------------
80 C     BLOCK 1.
81 C     INITIALIZATIONS.
82 C---------------------------------------------------
83 C
84 C***FIRST EXECUTABLE STATEMENT  DDAINI
85       IDID=1
86       NEF=0
87       NCF=0
88       NSF=0
89       XOLD=X
90       YNORM=DDANRM(NEQ,Y,WT,RPAR,IPAR)
91 C
92 C     SAVE Y AND YPRIME IN PHI
93       DO 100 I=1,NEQ
94          PHI(I,1)=Y(I)
95 100      PHI(I,2)=YPRIME(I)
96 C
97 C
98 C----------------------------------------------------
99 C     BLOCK 2.
100 C     DO ONE BACKWARD EULER STEP.
101 C----------------------------------------------------
102 C
103 C     SET UP FOR START OF CORRECTOR ITERATION
104 200   CJ=1.0D0/H
105       X=X+H
106 C
107 C     PREDICT SOLUTION AND DERIVATIVE
108       DO 250 I=1,NEQ
109 250     Y(I)=Y(I)+H*YPRIME(I)
110 C
111       JCALC=-1
112       M=0
113       CONVGD=.TRUE.
114 C
115 C
116 C     CORRECTOR LOOP.
117 300   IWM(LNRE)=IWM(LNRE)+1
118       IRES=0
119 C
120       CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
121       if(iero.ne.0) return
122       IF (IRES.LT.0) GO TO 430
123 C
124 C
125 C     EVALUATE THE ITERATION MATRIX
126       IF (JCALC.NE.-1) GO TO 310
127       IWM(LNJE)=IWM(LNJE)+1
128       JCALC=0
129       CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H,
130      *   IER,WT,E,WM,IWM,RES,IRES,
131      *   UROUND,JAC,RPAR,IPAR,NTEMP)
132       if(iero.ne.0) return
133 C
134       S=1000000.D0
135       IF (IRES.LT.0) GO TO 430
136       IF (IER.NE.0) GO TO 430
137       NSF=0
138 C
139 C
140 C
141 C     MULTIPLY RESIDUAL BY DAMPING FACTOR
142 310   CONTINUE
143       DO 320 I=1,NEQ
144 320      DELTA(I)=DELTA(I)*DAMP
145 C
146 C     COMPUTE A NEW ITERATE (BACK SUBSTITUTION)
147 C     STORE THE CORRECTION IN DELTA
148 C
149       CALL DDASLV(NEQ,DELTA,WM,IWM)
150 C
151 C     UPDATE Y AND YPRIME
152       DO 330 I=1,NEQ
153          Y(I)=Y(I)-DELTA(I)
154 330      YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
155 C
156 C     TEST FOR CONVERGENCE OF THE ITERATION.
157 C
158       DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
159       IF (DELNRM.LE.100.D0*UROUND*YNORM)
160      *   GO TO 400
161 C
162       IF (M.GT.0) GO TO 340
163          OLDNRM=DELNRM
164          GO TO 350
165 C
166 340   RATE=(DELNRM/OLDNRM)**(1.0D0/M)
167       IF (RATE.GT.0.90D0) GO TO 430
168       S=RATE/(1.0D0-RATE)
169 C
170 350   IF (S*DELNRM .LE. 0.33D0) GO TO 400
171 C
172 C
173 C     THE CORRECTOR HAS NOT YET CONVERGED. UPDATE
174 C     M AND AND TEST WHETHER THE MAXIMUM
175 C     NUMBER OF ITERATIONS HAVE BEEN TRIED.
176 C     EVERY MJAC ITERATIONS, GET A NEW
177 C     ITERATION MATRIX.
178 C
179       M=M+1
180       IF (M.GE.MAXIT) GO TO 430
181 C
182       IF ((M/MJAC)*MJAC.EQ.M) JCALC=-1
183       GO TO 300
184 C
185 C
186 C     THE ITERATION HAS CONVERGED.
187 C     CHECK NONNEGATIVITY CONSTRAINTS
188 400   IF (NONNEG.EQ.0) GO TO 450
189       DO 410 I=1,NEQ
190 410      DELTA(I)=MIN(Y(I),0.0D0)
191 C
192       DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
193       IF (DELNRM.GT.0.33D0) GO TO 430
194 C
195       DO 420 I=1,NEQ
196          Y(I)=Y(I)-DELTA(I)
197 420      YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
198       GO TO 450
199 C
200 C
201 C     EXITS FROM CORRECTOR LOOP.
202 430   CONVGD=.FALSE.
203 450   IF (.NOT.CONVGD) GO TO 600
204 C
205 C
206 C
207 C-----------------------------------------------------
208 C     BLOCK 3.
209 C     THE CORRECTOR ITERATION CONVERGED.
210 C     DO ERROR TEST.
211 C-----------------------------------------------------
212 C
213       DO 510 I=1,NEQ
214 510      E(I)=Y(I)-PHI(I,1)
215       IERR=DDANRM(NEQ,E,WT,RPAR,IPAR)
216 C
217       IF (IERR.LE.1.0D0) RETURN
218 C
219 C
220 C
221 C--------------------------------------------------------
222 C     BLOCK 4.
223 C     THE BACKWARD EULER STEP FAILED. RESTORE X, Y
224 C     AND YPRIME TO THEIR ORIGINAL VALUES.
225 C     REDUCE STEPSIZE AND TRY AGAIN, IF
226 C     POSSIBLE.
227 C---------------------------------------------------------
228 C
229 600   CONTINUE
230       X = XOLD
231       DO 610 I=1,NEQ
232          Y(I)=PHI(I,1)
233 610      YPRIME(I)=PHI(I,2)
234 C
235       IF (CONVGD) GO TO 640
236       IF (IER.EQ.0) GO TO 620
237          NSF=NSF+1
238          H=H*0.25D0
239          IF (NSF.LT.3.AND.ABS(H).GE.HMIN) GO TO 690
240          IDID=-12
241          RETURN
242 620   IF (IRES.GT.-2) GO TO 630
243          IDID=-12
244          RETURN
245 630   NCF=NCF+1
246       H=H*0.25D0
247       IF (NCF.LT.10.AND.ABS(H).GE.HMIN) GO TO 690
248          IDID=-12
249          RETURN
250 C
251 640   NEF=NEF+1
252       R=0.90D0/(2.0D0*IERR+0.0001D0)
253       R=MAX(0.1D0,MIN(0.5D0,R))
254       H=H*R
255       IF (ABS(H).GE.HMIN.AND.NEF.LT.10) GO TO 690
256          IDID=-12
257          RETURN
258 690      GO TO 200
259 C
260 C-------------END OF SUBROUTINE DDAINI----------------------
261       END
262       SUBROUTINE DDAJAC (NEQ, X, Y, YPRIME, DELTA, CJ, H,
263      +   IER, WT, E, WM, IWM, RES, IRES, UROUND, JAC, RPAR,
264      +   IPAR, NTEMP)
265       include 'stack.h'
266
267 C***BEGIN PROLOGUE  DDAJAC
268 C***SUBSIDIARY
269 C***PURPOSE  Compute the iteration matrix for DDASSL and form the
270 C            LU-decomposition.
271 C***LIBRARY   SLATEC (DASSL)
272 C***TYPE      DOUBLE PRECISION (SDAJAC-S, DDAJAC-D)
273 C***AUTHOR  PETZOLD, LINDA R., (LLNL)
274 C***DESCRIPTION
275 C-----------------------------------------------------------------------
276 C     THIS ROUTINE COMPUTES THE ITERATION MATRIX
277 C     PD=DG/DY+CJ*DG/DYPRIME (WHERE G(X,Y,YPRIME)=0).
278 C     HERE PD IS COMPUTED BY THE USER-SUPPLIED
279 C     ROUTINE JAC IF IWM(MTYPE) IS 1 OR 4, AND
280 C     IT IS COMPUTED BY NUMERICAL FINITE DIFFERENCING
281 C     IF IWM(MTYPE)IS 2 OR 5
282 C     THE PARAMETERS HAVE THE FOLLOWING MEANINGS.
283 C     Y        = ARRAY CONTAINING PREDICTED VALUES
284 C     YPRIME   = ARRAY CONTAINING PREDICTED DERIVATIVES
285 C     DELTA    = RESIDUAL EVALUATED AT (X,Y,YPRIME)
286 C                (USED ONLY IF IWM(MTYPE)=2 OR 5)
287 C     CJ       = SCALAR PARAMETER DEFINING ITERATION MATRIX
288 C     H        = CURRENT STEPSIZE IN INTEGRATION
289 C     IER      = VARIABLE WHICH IS .NE. 0
290 C                IF ITERATION MATRIX IS SINGULAR,
291 C                AND 0 OTHERWISE.
292 C     WT       = VECTOR OF WEIGHTS FOR COMPUTING NORMS
293 C     E        = WORK SPACE (TEMPORARY) OF LENGTH NEQ
294 C     WM       = REAL WORK SPACE FOR MATRICES. ON
295 C                OUTPUT IT CONTAINS THE LU DECOMPOSITION
296 C                OF THE ITERATION MATRIX.
297 C     IWM      = INTEGER WORK SPACE CONTAINING
298 C                MATRIX INFORMATION
299 C     RES      = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE
300 C                TO EVALUATE THE RESIDUAL FUNCTION G(X,Y,YPRIME)
301 C     IRES     = FLAG WHICH IS EQUAL TO ZERO IF NO ILLEGAL VALUES
302 C                IN RES, AND LESS THAN ZERO OTHERWISE.  (IF IRES
303 C                IS LESS THAN ZERO, THE MATRIX WAS NOT COMPLETED)
304 C                IN THIS CASE (IF IRES .LT. 0), THEN IER = 0.
305 C     UROUND   = THE UNIT ROUNDOFF ERROR OF THE MACHINE BEING USED.
306 C     JAC      = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE
307 C                TO EVALUATE THE ITERATION MATRIX (THIS ROUTINE
308 C                IS ONLY USED IF IWM(MTYPE) IS 1 OR 4)
309 C-----------------------------------------------------------------------
310 C***ROUTINES CALLED  DGBFA, DGEFA
311 C***REVISION HISTORY  (YYMMDD)
312 C   830315  DATE WRITTEN
313 C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
314 C   901010  Modified three MAX calls to be all on one line.  (FNF)
315 C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
316 C   901026  Added explicit declarations for all variables and minor
317 C           cosmetic changes to prologue.  (FNF)
318 C   901101  Corrected PURPOSE.  (FNF)
319 C***END PROLOGUE  DDAJAC
320 C
321       INTEGER  NEQ, IER, IWM(*), IRES, IPAR(*), NTEMP
322       DOUBLE PRECISION
323      *   X, Y(*), YPRIME(*), DELTA(*), CJ, H, WT(*), E(*), WM(*),
324      *   UROUND, RPAR(*)
325       EXTERNAL  RES, JAC
326 C
327       EXTERNAL  DGBFA, DGEFA
328 C
329       INTEGER  I, I1, I2, II, IPSAVE, ISAVE, J, K, L, LENPD, LIPVT,
330      *   LML, LMTYPE, LMU, MBA, MBAND, MEB1, MEBAND, MSAVE, MTYPE, N,
331      *   NPD, NPDM1, NROW
332       DOUBLE PRECISION  DEL, DELINV, SQUR, YPSAVE, YSAVE
333 C
334       PARAMETER (NPD=1)
335       PARAMETER (LML=1)
336       PARAMETER (LMU=2)
337       PARAMETER (LMTYPE=4)
338       PARAMETER (LIPVT=21)
339 C
340 C***FIRST EXECUTABLE STATEMENT  DDAJAC
341       IER = 0
342       NPDM1=NPD-1
343       MTYPE=IWM(LMTYPE)
344       GO TO (100,200,300,400,500),MTYPE
345 C
346 C
347 C     DENSE USER-SUPPLIED MATRIX
348 100   LENPD=NEQ*NEQ
349       DO 110 I=1,LENPD
350 110      WM(NPDM1+I)=0.0D0
351       CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR)
352       if(iero.ne.0) return
353       GO TO 230
354 C
355 C
356 C     DENSE FINITE-DIFFERENCE-GENERATED MATRIX
357 200   IRES=0
358       NROW=NPDM1
359       SQUR = SQRT(UROUND)
360       DO 210 I=1,NEQ
361          DEL=SQUR*MAX(ABS(Y(I)),ABS(H*YPRIME(I)),ABS(WT(I)))
362          DEL=SIGN(DEL,H*YPRIME(I))
363          DEL=(Y(I)+DEL)-Y(I)
364          YSAVE=Y(I)
365          YPSAVE=YPRIME(I)
366          Y(I)=Y(I)+DEL
367          YPRIME(I)=YPRIME(I)+CJ*DEL
368          CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR)
369          if(iero.ne.0) return
370          IF (IRES .LT. 0) RETURN
371          DELINV=1.0D0/DEL
372          DO 220 L=1,NEQ
373 220      WM(NROW+L)=(E(L)-DELTA(L))*DELINV
374       NROW=NROW+NEQ
375       Y(I)=YSAVE
376       YPRIME(I)=YPSAVE
377 210   CONTINUE
378 C
379 C
380 C     DO DENSE-MATRIX LU DECOMPOSITION ON PD
381 230      CALL DGEFA(WM(NPD),NEQ,NEQ,IWM(LIPVT),IER)
382       RETURN
383 C
384 C
385 C     DUMMY SECTION FOR IWM(MTYPE)=3
386 300   RETURN
387 C
388 C
389 C     BANDED USER-SUPPLIED MATRIX
390 400   LENPD=(2*IWM(LML)+IWM(LMU)+1)*NEQ
391       DO 410 I=1,LENPD
392 410      WM(NPDM1+I)=0.0D0
393       CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR)
394       if(iero.ne.0) return
395       MEBAND=2*IWM(LML)+IWM(LMU)+1
396       GO TO 550
397 C
398 C
399 C     BANDED FINITE-DIFFERENCE-GENERATED MATRIX
400 500   MBAND=IWM(LML)+IWM(LMU)+1
401       MBA=MIN(MBAND,NEQ)
402       MEBAND=MBAND+IWM(LML)
403       MEB1=MEBAND-1
404       MSAVE=(NEQ/MBAND)+1
405       ISAVE=NTEMP-1
406       IPSAVE=ISAVE+MSAVE
407       IRES=0
408       SQUR=SQRT(UROUND)
409       DO 540 J=1,MBA
410          DO 510 N=J,NEQ,MBAND
411           K= (N-J)/MBAND + 1
412           WM(ISAVE+K)=Y(N)
413           WM(IPSAVE+K)=YPRIME(N)
414           DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),ABS(WT(N)))
415           DEL=SIGN(DEL,H*YPRIME(N))
416           DEL=(Y(N)+DEL)-Y(N)
417           Y(N)=Y(N)+DEL
418 510       YPRIME(N)=YPRIME(N)+CJ*DEL
419       CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR)
420       if(iero.ne.0) return
421       IF (IRES .LT. 0) RETURN
422       DO 530 N=J,NEQ,MBAND
423           K= (N-J)/MBAND + 1
424           Y(N)=WM(ISAVE+K)
425           YPRIME(N)=WM(IPSAVE+K)
426           DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),ABS(WT(N)))
427           DEL=SIGN(DEL,H*YPRIME(N))
428           DEL=(Y(N)+DEL)-Y(N)
429           DELINV=1.0D0/DEL
430           I1=MAX(1,(N-IWM(LMU)))
431           I2=MIN(NEQ,(N+IWM(LML)))
432           II=N*MEB1-IWM(LML)+NPDM1
433           DO 520 I=I1,I2
434 520         WM(II+I)=(E(I)-DELTA(I))*DELINV
435 530      CONTINUE
436 540   CONTINUE
437 C
438 C
439 C     DO LU DECOMPOSITION OF BANDED PD
440 550   CALL DGBFA(WM(NPD),MEBAND,NEQ,
441      *    IWM(LML),IWM(LMU),IWM(LIPVT),IER)
442       RETURN
443 C------END OF SUBROUTINE DDAJAC------
444       END
445       DOUBLE PRECISION FUNCTION DDANRM (NEQ, V, WT, RPAR, IPAR)
446 C***BEGIN PROLOGUE  DDANRM
447 C***SUBSIDIARY
448 C***PURPOSE  Compute vector norm for DDASSL.
449 C***LIBRARY   SLATEC (DASSL)
450 C***TYPE      DOUBLE PRECISION (SDANRM-S, DDANRM-D)
451 C***AUTHOR  PETZOLD, LINDA R., (LLNL)
452 C***DESCRIPTION
453 C-----------------------------------------------------------------------
454 C     THIS FUNCTION ROUTINE COMPUTES THE WEIGHTED
455 C     ROOT-MEAN-SQUARE NORM OF THE VECTOR OF LENGTH
456 C     NEQ CONTAINED IN THE ARRAY V,WITH WEIGHTS
457 C     CONTAINED IN THE ARRAY WT OF LENGTH NEQ.
458 C        DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2)
459 C-----------------------------------------------------------------------
460 C***ROUTINES CALLED  (NONE)
461 C***REVISION HISTORY  (YYMMDD)
462 C   830315  DATE WRITTEN
463 C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
464 C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
465 C   901026  Added explicit declarations for all variables and minor
466 C           cosmetic changes to prologue.  (FNF)
467 C***END PROLOGUE  DDANRM
468 C
469       INTEGER  NEQ, IPAR(*)
470       DOUBLE PRECISION  V(NEQ), WT(NEQ), RPAR(*)
471 C
472       INTEGER  I
473       DOUBLE PRECISION  SUM, VMAX
474 C
475 C***FIRST EXECUTABLE STATEMENT  DDANRM
476       DDANRM = 0.0D0
477       VMAX = 0.0D0
478       DO 10 I = 1,NEQ
479         IF(ABS(V(I)/WT(I)) .GT. VMAX) VMAX = ABS(V(I)/WT(I))
480 10      CONTINUE
481       IF(VMAX .LE. 0.0D0) GO TO 30
482       SUM = 0.0D0
483       DO 20 I = 1,NEQ
484 20      SUM = SUM + ((V(I)/WT(I))/VMAX)**2
485       DDANRM = VMAX*SQRT(SUM/NEQ)
486 30    CONTINUE
487       RETURN
488 C------END OF FUNCTION DDANRM------
489       END
490       SUBROUTINE DDASLV (NEQ, DELTA, WM, IWM)
491 C***BEGIN PROLOGUE  DDASLV
492 C***SUBSIDIARY
493 C***PURPOSE  Linear system solver for DDASSL.
494 C***LIBRARY   SLATEC (DASSL)
495 C***TYPE      DOUBLE PRECISION (SDASLV-S, DDASLV-D)
496 C***AUTHOR  PETZOLD, LINDA R., (LLNL)
497 C***DESCRIPTION
498 C-----------------------------------------------------------------------
499 C     THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR
500 C     SYSTEM ARISING IN THE NEWTON ITERATION.
501 C     MATRICES AND REAL TEMPORARY STORAGE AND
502 C     REAL INFORMATION ARE STORED IN THE ARRAY WM.
503 C     INTEGER MATRIX INFORMATION IS STORED IN
504 C     THE ARRAY IWM.
505 C     FOR A DENSE MATRIX, THE LINPACK ROUTINE
506 C     DGESL IS CALLED.
507 C     FOR A BANDED MATRIX,THE LINPACK ROUTINE
508 C     DGBSL IS CALLED.
509 C-----------------------------------------------------------------------
510 C***ROUTINES CALLED  DGBSL, DGESL
511 C***REVISION HISTORY  (YYMMDD)
512 C   830315  DATE WRITTEN
513 C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
514 C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
515 C   901026  Added explicit declarations for all variables and minor
516 C           cosmetic changes to prologue.  (FNF)
517 C***END PROLOGUE  DDASLV
518 C
519       INTEGER  NEQ, IWM(*)
520       DOUBLE PRECISION  DELTA(*), WM(*)
521 C
522       EXTERNAL  DGBSL, DGESL
523 C
524       INTEGER  LIPVT, LML, LMU, LMTYPE, MEBAND, MTYPE, NPD
525       PARAMETER (NPD=1)
526       PARAMETER (LML=1)
527       PARAMETER (LMU=2)
528       PARAMETER (LMTYPE=4)
529       PARAMETER (LIPVT=21)
530 C
531 C***FIRST EXECUTABLE STATEMENT  DDASLV
532       MTYPE=IWM(LMTYPE)
533       GO TO(100,100,300,400,400),MTYPE
534 C
535 C     DENSE MATRIX
536 100   CALL DGESL(WM(NPD),NEQ,NEQ,IWM(LIPVT),DELTA,0)
537       RETURN
538 C
539 C     DUMMY SECTION FOR MTYPE=3
540 300   CONTINUE
541       RETURN
542 C
543 C     BANDED MATRIX
544 400   MEBAND=2*IWM(LML)+IWM(LMU)+1
545       CALL DGBSL(WM(NPD),MEBAND,NEQ,IWM(LML),
546      *  IWM(LMU),IWM(LIPVT),DELTA,0)
547       RETURN
548 C------END OF SUBROUTINE DDASLV------
549       END
550       SUBROUTINE DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
551      +   IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC)
552       include 'stack.h'
553 C***BEGIN PROLOGUE  DDASSL
554 C***PURPOSE  This code solves a system of differential/algebraic
555 C            equations of the form G(T,Y,YPRIME) = 0.
556 C***LIBRARY   SLATEC (DASSL)
557 C***CATEGORY  I1A2
558 C***TYPE      DOUBLE PRECISION (SDASSL-S, DDASSL-D)
559 C***KEYWORDS  DIFFERENTIAL/ALGEBRAIC, BACKWARD DIFFERENTIATION FORMULAS,
560 C             IMPLICIT DIFFERENTIAL SYSTEMS
561 C***AUTHOR  PETZOLD, LINDA R., (LLNL)
562 C             COMPUTING AND MATHEMATICS RESEARCH DIVISION
563 C             LAWRENCE LIVERMORE NATIONAL LABORATORY
564 C             L - 316, P.O. BOX 808,
565 C             LIVERMORE, CA.    94550
566 C***DESCRIPTION
567 C
568 C *Usage:
569 C
570 C      EXTERNAL RES, JAC
571 C      INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR
572 C      DOUBLE PRECISION T, Y(NEQ), YPRIME(NEQ), TOUT, RTOL, ATOL,
573 C     *   RWORK(LRW), RPAR
574 C
575 C      CALL DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
576 C     *   IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC)
577 C
578 C
579 C *Arguments:
580 C  (In the following, all real arrays should be type DOUBLE PRECISION.)
581 C
582 C  RES:EXT     This is a subroutine which you provide to define the
583 C              differential/algebraic system.
584 C
585 C  NEQ:IN      This is the number of equations to be solved.
586 C
587 C  T:INOUT     This is the current value of the independent variable.
588 C
589 C  Y(*):INOUT  This array contains the solution components at T.
590 C
591 C  YPRIME(*):INOUT  This array contains the derivatives of the solution
592 C              components at T.
593 C
594 C  TOUT:IN     This is a point at which a solution is desired.
595 C
596 C  INFO(N):IN  The basic task of the code is to solve the system from T
597 C              to TOUT and return an answer at TOUT.  INFO is an integer
598 C              array which is used to communicate exactly how you want
599 C              this task to be carried out.  (See below for details.)
600 C              N must be greater than or equal to 15.
601 C
602 C  RTOL,ATOL:INOUT  These quantities represent relative and absolute
603 C              error tolerances which you provide to indicate how
604 C              accurately you wish the solution to be computed.  You
605 C              may choose them to be both scalars or else both vectors.
606 C              Caution:  In Fortran 77, a scalar is not the same as an
607 C                        array of length 1.  Some compilers may object
608 C                        to using scalars for RTOL,ATOL.
609 C
610 C  IDID:OUT    This scalar quantity is an indicator reporting what the
611 C              code did.  You must monitor this integer variable to
612 C              decide  what action to take next.
613 C
614 C  RWORK:WORK  A real work array of length LRW which provides the
615 C              code with needed storage space.
616 C
617 C  LRW:IN      The length of RWORK.  (See below for required length.)
618 C
619 C  IWORK:WORK  An integer work array of length LIW which probides the
620 C              code with needed storage space.
621 C
622 C  LIW:IN      The length of IWORK.  (See below for required length.)
623 C
624 C  RPAR,IPAR:IN  These are real and integer parameter arrays which
625 C              you can use for communication between your calling
626 C              program and the RES subroutine (and the JAC subroutine)
627 C
628 C  JAC:EXT     This is the name of a subroutine which you may choose
629 C              to provide for defining a matrix of partial derivatives
630 C              described below.
631 C
632 C  Quantities which may be altered by DDASSL are:
633 C     T, Y(*), YPRIME(*), INFO(1), RTOL, ATOL,
634 C     IDID, RWORK(*) AND IWORK(*)
635 C
636 C *Description
637 C
638 C  Subroutine DDASSL uses the backward differentiation formulas of
639 C  orders one through five to solve a system of the above form for Y and
640 C  YPRIME.  Values for Y and YPRIME at the initial time must be given as
641 C  input.  These values must be consistent, (that is, if T,Y,YPRIME are
642 C  the given initial values, they must satisfy G(T,Y,YPRIME) = 0.).  The
643 C  subroutine solves the system from T to TOUT.  It is easy to continue
644 C  the solution to get results at additional TOUT.  This is the interval
645 C  mode of operation.  Intermediate results can also be obtained easily
646 C  by using the intermediate-output capability.
647 C
648 C  The following detailed description is divided into subsections:
649 C    1. Input required for the first call to DDASSL.
650 C    2. Output after any return from DDASSL.
651 C    3. What to do to continue the integration.
652 C    4. Error messages.
653 C
654 C
655 C  -------- INPUT -- WHAT TO DO ON THE FIRST CALL TO DDASSL ------------
656 C
657 C  The first call of the code is defined to be the start of each new
658 C  problem. Read through the descriptions of all the following items,
659 C  provide sufficient storage space for designated arrays, set
660 C  appropriate variables for the initialization of the problem, and
661 C  give information about how you want the problem to be solved.
662 C
663 C
664 C  RES -- Provide a subroutine of the form
665 C             SUBROUTINE RES(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
666 C         to define the system of differential/algebraic
667 C         equations which is to be solved. For the given values
668 C         of T,Y and YPRIME, the subroutine should
669 C         return the residual of the defferential/algebraic
670 C         system
671 C             DELTA = G(T,Y,YPRIME)
672 C         (DELTA(*) is a vector of length NEQ which is
673 C         output for RES.)
674 C
675 C         Subroutine RES must not alter T,Y or YPRIME.
676 C         You must declare the name RES in an external
677 C         statement in your program that calls DDASSL.
678 C         You must dimension Y,YPRIME and DELTA in RES.
679 C
680 C         IRES is an integer flag which is always equal to
681 C         zero on input. Subroutine RES should alter IRES
682 C         only if it encounters an illegal value of Y or
683 C         a stop condition. Set IRES = -1 if an input value
684 C         is illegal, and DDASSL will try to solve the problem
685 C         without getting IRES = -1. If IRES = -2, DDASSL
686 C         will return control to the calling program
687 C         with IDID = -11.
688 C
689 C         RPAR and IPAR are real and integer parameter arrays which
690 C         you can use for communication between your calling program
691 C         and subroutine RES. They are not altered by DDASSL. If you
692 C         do not need RPAR or IPAR, ignore these parameters by treat-
693 C         ing them as dummy arguments. If you do choose to use them,
694 C         dimension them in your calling program and in RES as arrays
695 C         of appropriate length.
696 C
697 C  NEQ -- Set it to the number of differential equations.
698 C         (NEQ .GE. 1)
699 C
700 C  T -- Set it to the initial point of the integration.
701 C         T must be defined as a variable.
702 C
703 C  Y(*) -- Set this vector to the initial values of the NEQ solution
704 C         components at the initial point. You must dimension Y of
705 C         length at least NEQ in your calling program.
706 C
707 C  YPRIME(*) -- Set this vector to the initial values of the NEQ
708 C         first derivatives of the solution components at the initial
709 C         point.  You must dimension YPRIME at least NEQ in your
710 C         calling program. If you do not know initial values of some
711 C         of the solution components, see the explanation of INFO(11).
712 C
713 C  TOUT -- Set it to the first point at which a solution
714 C         is desired. You can not take TOUT = T.
715 C         integration either forward in T (TOUT .GT. T) or
716 C         backward in T (TOUT .LT. T) is permitted.
717 C
718 C         The code advances the solution from T to TOUT using
719 C         step sizes which are automatically selected so as to
720 C         achieve the desired accuracy. If you wish, the code will
721 C         return with the solution and its derivative at
722 C         intermediate steps (intermediate-output mode) so that
723 C         you can monitor them, but you still must provide TOUT in
724 C         accord with the basic aim of the code.
725 C
726 C         The first step taken by the code is a critical one
727 C         because it must reflect how fast the solution changes near
728 C         the initial point. The code automatically selects an
729 C         initial step size which is practically always suitable for
730 C         the problem. By using the fact that the code will not step
731 C         past TOUT in the first step, you could, if necessary,
732 C         restrict the length of the initial step size.
733 C
734 C         For some problems it may not be permissible to integrate
735 C         past a point TSTOP because a discontinuity occurs there
736 C         or the solution or its derivative is not defined beyond
737 C         TSTOP. When you have declared a TSTOP point (SEE INFO(4)
738 C         and RWORK(1)), you have told the code not to integrate
739 C         past TSTOP. In this case any TOUT beyond TSTOP is invalid
740 C         input.
741 C
742 C  INFO(*) -- Use the INFO array to give the code more details about
743 C         how you want your problem solved.  This array should be
744 C         dimensioned of length 15, though DDASSL uses only the first
745 C         eleven entries.  You must respond to all of the following
746 C         items, which are arranged as questions.  The simplest use
747 C         of the code corresponds to answering all questions as yes,
748 C         i.e. setting all entries of INFO to 0.
749 C
750 C       INFO(1) - This parameter enables the code to initialize
751 C              itself. You must set it to indicate the start of every
752 C              new problem.
753 C
754 C          **** Is this the first call for this problem ...
755 C                Yes - Set INFO(1) = 0
756 C                 No - Not applicable here.
757 C                      See below for continuation calls.  ****
758 C
759 C       INFO(2) - How much accuracy you want of your solution
760 C              is specified by the error tolerances RTOL and ATOL.
761 C              The simplest use is to take them both to be scalars.
762 C              To obtain more flexibility, they can both be vectors.
763 C              The code must be told your choice.
764 C
765 C          **** Are both error tolerances RTOL, ATOL scalars ...
766 C                Yes - Set INFO(2) = 0
767 C                      and input scalars for both RTOL and ATOL
768 C                 No - Set INFO(2) = 1
769 C                      and input arrays for both RTOL and ATOL ****
770 C
771 C       INFO(3) - The code integrates from T in the direction
772 C              of TOUT by steps. If you wish, it will return the
773 C              computed solution and derivative at the next
774 C              intermediate step (the intermediate-output mode) or
775 C              TOUT, whichever comes first. This is a good way to
776 C              proceed if you want to see the behavior of the solution.
777 C              If you must have solutions at a great many specific
778 C              TOUT points, this code will compute them efficiently.
779 C
780 C          **** Do you want the solution only at
781 C                TOUT (and not at the next intermediate step) ...
782 C                 Yes - Set INFO(3) = 0
783 C                  No - Set INFO(3) = 1 ****
784 C
785 C       INFO(4) - To handle solutions at a great many specific
786 C              values TOUT efficiently, this code may integrate past
787 C              TOUT and interpolate to obtain the result at TOUT.
788 C              Sometimes it is not possible to integrate beyond some
789 C              point TSTOP because the equation changes there or it is
790 C              not defined past TSTOP. Then you must tell the code
791 C              not to go past.
792 C
793 C           **** Can the integration be carried out without any
794 C                restrictions on the independent variable T ...
795 C                 Yes - Set INFO(4)=0
796 C                  No - Set INFO(4)=1
797 C                       and define the stopping point TSTOP by
798 C                       setting RWORK(1)=TSTOP ****
799 C
800 C       INFO(5) - To solve differential/algebraic problems it is
801 C              necessary to use a matrix of partial derivatives of the
802 C              system of differential equations. If you do not
803 C              provide a subroutine to evaluate it analytically (see
804 C              description of the item JAC in the call list), it will
805 C              be approximated by numerical differencing in this code.
806 C              although it is less trouble for you to have the code
807 C              compute partial derivatives by numerical differencing,
808 C              the solution will be more reliable if you provide the
809 C              derivatives via JAC. Sometimes numerical differencing
810 C              is cheaper than evaluating derivatives in JAC and
811 C              sometimes it is not - this depends on your problem.
812 C
813 C           **** Do you want the code to evaluate the partial
814 C                derivatives automatically by numerical differences ...
815 C                   Yes - Set INFO(5)=0
816 C                    No - Set INFO(5)=1
817 C                  and provide subroutine JAC for evaluating the
818 C                  matrix of partial derivatives ****
819 C
820 C       INFO(6) - DDASSL will perform much better if the matrix of
821 C              partial derivatives, DG/DY + CJ*DG/DYPRIME,
822 C              (here CJ is a scalar determined by DDASSL)
823 C              is banded and the code is told this. In this
824 C              case, the storage needed will be greatly reduced,
825 C              numerical differencing will be performed much cheaper,
826 C              and a number of important algorithms will execute much
827 C              faster. The differential equation is said to have
828 C              half-bandwidths ML (lower) and MU (upper) if equation i
829 C              involves only unknowns Y(J) with
830 C                             I-ML .LE. J .LE. I+MU
831 C              for all I=1,2,...,NEQ. Thus, ML and MU are the widths
832 C              of the lower and upper parts of the band, respectively,
833 C              with the main diagonal being excluded. If you do not
834 C              indicate that the equation has a banded matrix of partial
835 C              derivatives, the code works with a full matrix of NEQ**2
836 C              elements (stored in the conventional way). Computations
837 C              with banded matrices cost less time and storage than with
838 C              full matrices if 2*ML+MU .LT. NEQ. If you tell the
839 C              code that the matrix of partial derivatives has a banded
840 C              structure and you want to provide subroutine JAC to
841 C              compute the partial derivatives, then you must be careful
842 C              to store the elements of the matrix in the special form
843 C              indicated in the description of JAC.
844 C
845 C          **** Do you want to solve the problem using a full
846 C               (dense) matrix (and not a special banded
847 C               structure) ...
848 C                Yes - Set INFO(6)=0
849 C                 No - Set INFO(6)=1
850 C                       and provide the lower (ML) and upper (MU)
851 C                       bandwidths by setting
852 C                       IWORK(1)=ML
853 C                       IWORK(2)=MU ****
854 C
855 C
856 C        INFO(7) -- You can specify a maximum (absolute value of)
857 C              stepsize, so that the code
858 C              will avoid passing over very
859 C              large regions.
860 C
861 C          ****  Do you want the code to decide
862 C                on its own maximum stepsize?
863 C                Yes - Set INFO(7)=0
864 C                 No - Set INFO(7)=1
865 C                      and define HMAX by setting
866 C                      RWORK(2)=HMAX ****
867 C
868 C        INFO(8) -- Differential/algebraic problems
869 C              may occaisionally suffer from
870 C              severe scaling difficulties on the
871 C              first step. If you know a great deal
872 C              about the scaling of your problem, you can
873 C              help to alleviate this problem by
874 C              specifying an initial stepsize HO.
875 C
876 C          ****  Do you want the code to define
877 C                its own initial stepsize?
878 C                Yes - Set INFO(8)=0
879 C                 No - Set INFO(8)=1
880 C                      and define HO by setting
881 C                      RWORK(3)=HO ****
882 C
883 C        INFO(9) -- If storage is a severe problem,
884 C              you can save some locations by
885 C              restricting the maximum order MAXORD.
886 C              the default value is 5. for each
887 C              order decrease below 5, the code
888 C              requires NEQ fewer locations, however
889 C              it is likely to be slower. In any
890 C              case, you must have 1 .LE. MAXORD .LE. 5
891 C          ****  Do you want the maximum order to
892 C                default to 5?
893 C                Yes - Set INFO(9)=0
894 C                 No - Set INFO(9)=1
895 C                      and define MAXORD by setting
896 C                      IWORK(3)=MAXORD ****
897 C
898 C        INFO(10) --If you know that the solutions to your equations
899 C               will always be nonnegative, it may help to set this
900 C               parameter. However, it is probably best to
901 C               try the code without using this option first,
902 C               and only to use this option if that doesn't
903 C               work very well.
904 C           ****  Do you want the code to solve the problem without
905 C                 invoking any special nonnegativity constraints?
906 C                  Yes - Set INFO(10)=0
907 C                   No - Set INFO(10)=1
908 C
909 C        INFO(11) --DDASSL normally requires the initial T,
910 C               Y, and YPRIME to be consistent. That is,
911 C               you must have G(T,Y,YPRIME) = 0 at the initial
912 C               time. If you do not know the initial
913 C               derivative precisely, you can let DDASSL try
914 C               to compute it.
915 C          ****   Are the initialHE INITIAL T, Y, YPRIME consistent?
916 C                 Yes - Set INFO(11) = 0
917 C                  No - Set INFO(11) = 1,
918 C                       and set YPRIME to an initial approximation
919 C                       to YPRIME.  (If you have no idea what
920 C                       YPRIME should be, set it to zero. Note
921 C                       that the initial Y should be such
922 C                       that there must exist a YPRIME so that
923 C                       G(T,Y,YPRIME) = 0.)
924 C
925 C  RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL
926 C         error tolerances to tell the code how accurately you
927 C         want the solution to be computed.  They must be defined
928 C         as variables because the code may change them.  You
929 C         have two choices --
930 C               Both RTOL and ATOL are scalars. (INFO(2)=0)
931 C               Both RTOL and ATOL are vectors. (INFO(2)=1)
932 C         in either case all components must be non-negative.
933 C
934 C         The tolerances are used by the code in a local error
935 C         test at each step which requires roughly that
936 C               ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL
937 C         for each vector component.
938 C         (More specifically, a root-mean-square norm is used to
939 C         measure the size of vectors, and the error test uses the
940 C         magnitude of the solution at the beginning of the step.)
941 C
942 C         The true (global) error is the difference between the
943 C         true solution of the initial value problem and the
944 C         computed approximation.  Practically all present day
945 C         codes, including this one, control the local error at
946 C         each step and do not even attempt to control the global
947 C         error directly.
948 C         Usually, but not always, the true accuracy of the
949 C         computed Y is comparable to the error tolerances. This
950 C         code will usually, but not always, deliver a more
951 C         accurate solution if you reduce the tolerances and
952 C         integrate again.  By comparing two such solutions you
953 C         can get a fairly reliable idea of the true error in the
954 C         solution at the bigger tolerances.
955 C
956 C         Setting ATOL=0. results in a pure relative error test on
957 C         that component.  Setting RTOL=0. results in a pure
958 C         absolute error test on that component.  A mixed test
959 C         with non-zero RTOL and ATOL corresponds roughly to a
960 C         relative error test when the solution component is much
961 C         bigger than ATOL and to an absolute error test when the
962 C         solution component is smaller than the threshhold ATOL.
963 C
964 C         The code will not attempt to compute a solution at an
965 C         accuracy unreasonable for the machine being used.  It will
966 C         advise you if you ask for too much accuracy and inform
967 C         you as to the maximum accuracy it believes possible.
968 C
969 C  RWORK(*) --  Dimension this real work array of length LRW in your
970 C         calling program.
971 C
972 C  LRW -- Set it to the declared length of the RWORK array.
973 C               You must have
974 C                    LRW .GE. 40+(MAXORD+4)*NEQ+NEQ**2
975 C               for the full (dense) JACOBIAN case (when INFO(6)=0), or
976 C                    LRW .GE. 40+(MAXORfD+4)*NEQ+(2*ML+MU+1)*NEQ
977 C               for the banded user-defined JACOBIAN case
978 C               (when INFO(5)=1 and INFO(6)=1), or
979 C                     LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ
980 C                           +2*(NEQ/(ML+MU+1)+1)
981 C               for the banded finite-difference-generated JACOBIAN case
982 C               (when INFO(5)=0 and INFO(6)=1)
983 C
984 C  IWORK(*) --  Dimension this integer work array of length LIW in
985 C         your calling program.
986 C
987 C  LIW -- Set it to the declared length of the IWORK array.
988 C               You must have LIW .GE. 20+NEQ
989 C
990 C  RPAR, IPAR -- These are parameter arrays, of real and integer
991 C         type, respectively.  You can use them for communication
992 C         between your program that calls DDASSL and the
993 C         RES subroutine (and the JAC subroutine).  They are not
994 C         altered by DDASSL.  If you do not need RPAR or IPAR,
995 C         ignore these parameters by treating them as dummy
996 C         arguments.  If you do choose to use them, dimension
997 C         them in your calling program and in RES (and in JAC)
998 C         as arrays of appropriate length.
999 C
1000 C  JAC -- If you have set INFO(5)=0, you can ignore this parameter
1001 C         by treating it as a dummy argument.  Otherwise, you must
1002 C         provide a subroutine of the form
1003 C             SUBROUTINE JAC(T,Y,YPRIME,PD,CJ,RPAR,IPAR)
1004 C         to define the matrix of partial derivatives
1005 C             PD=DG/DY+CJ*DG/DYPRIME
1006 C         CJ is a scalar which is input to JAC.
1007 C         For the given values of T,Y,YPRIME, the
1008 C         subroutine must evaluate the non-zero partial
1009 C         derivatives for each equation and each solution
1010 C         component, and store these values in the
1011 C         matrix PD.  The elements of PD are set to zero
1012 C         before each call to JAC so only non-zero elements
1013 C         need to be defined.
1014 C
1015 C         Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ.
1016 C         You must declare the name JAC in an EXTERNAL statement in
1017 C         your program that calls DDASSL.  You must dimension Y,
1018 C         YPRIME and PD in JAC.
1019 C
1020 C         The way you must store the elements into the PD matrix
1021 C         depends on the structure of the matrix which you
1022 C         indicated by INFO(6).
1023 C               *** INFO(6)=0 -- Full (dense) matrix ***
1024 C                   Give PD a first dimension of NEQ.
1025 C                   When you evaluate the (non-zero) partial derivative
1026 C                   of equation I with respect to variable J, you must
1027 C                   store it in PD according to
1028 C                   PD(I,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)"
1029 C               *** INFO(6)=1 -- Banded JACOBIAN with ML lower and MU
1030 C                   upper diagonal bands (refer to INFO(6) description
1031 C                   of ML and MU) ***
1032 C                   Give PD a first dimension of 2*ML+MU+1.
1033 C                   when you evaluate the (non-zero) partial derivative
1034 C                   of equation I with respect to variable J, you must
1035 C                   store it in PD according to
1036 C                   IROW = I - J + ML + MU + 1
1037 C                   PD(IROW,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)"
1038 C
1039 C         RPAR and IPAR are real and integer parameter arrays
1040 C         which you can use for communication between your calling
1041 C         program and your JACOBIAN subroutine JAC. They are not
1042 C         altered by DDASSL. If you do not need RPAR or IPAR,
1043 C         ignore these parameters by treating them as dummy
1044 C         arguments. If you do choose to use them, dimension
1045 C         them in your calling program and in JAC as arrays of
1046 C         appropriate length.
1047 C
1048 C
1049 C  OPTIONALLY REPLACEABLE NORM ROUTINE:
1050 C
1051 C     DDASSL uses a weighted norm DDANRM to measure the size
1052 C     of vectors such as the estimated error in each step.
1053 C     A FUNCTION subprogram
1054 C       DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR)
1055 C       DIMENSION V(NEQ),WT(NEQ)
1056 C     is used to define this norm. Here, V is the vector
1057 C     whose norm is to be computed, and WT is a vector of
1058 C     weights.  A DDANRM routine has been included with DDASSL
1059 C     which computes the weighted root-mean-square norm
1060 C     given by
1061 C       DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2)
1062 C     this norm is suitable for most problems. In some
1063 C     special cases, it may be more convenient and/or
1064 C     efficient to define your own norm by writing a function
1065 C     subprogram to be called instead of DDANRM. This should,
1066 C     however, be attempted only after careful thought and
1067 C     consideration.
1068 C
1069 C
1070 C  -------- OUTPUT -- AFTER ANY RETURN FROM DDASSL ---------------------
1071 C
1072 C  The principal aim of the code is to return a computed solution at
1073 C  TOUT, although it is also possible to obtain intermediate results
1074 C  along the way. To find out whether the code achieved its goal
1075 C  or if the integration process was interrupted before the task was
1076 C  completed, you must check the IDID parameter.
1077 C
1078 C
1079 C  T -- The solution was successfully advanced to the
1080 C               output value of T.
1081 C
1082 C  Y(*) -- Contains the computed solution approximation at T.
1083 C
1084 C  YPRIME(*) -- Contains the computed derivative
1085 C               approximation at T.
1086 C
1087 C  IDID -- Reports what the code did.
1088 C
1089 C                     *** Task completed ***
1090 C                Reported by positive values of IDID
1091 C
1092 C           IDID = 1 -- A step was successfully taken in the
1093 C                   intermediate-output mode. The code has not
1094 C                   yet reached TOUT.
1095 C
1096 C           IDID = 2 -- The integration to TSTOP was successfully
1097 C                   completed (T=TSTOP) by stepping exactly to TSTOP.
1098 C
1099 C           IDID = 3 -- The integration to TOUT was successfully
1100 C                   completed (T=TOUT) by stepping past TOUT.
1101 C                   Y(*) is obtained by interpolation.
1102 C                   YPRIME(*) is obtained by interpolation.
1103 C
1104 C                    *** Task interrupted ***
1105 C                Reported by negative values of IDID
1106 C
1107 C           IDID = -1 -- A large amount of work has been expended.
1108 C                   (About 500 steps)
1109 C
1110 C           IDID = -2 -- The error tolerances are too stringent.
1111 C
1112 C           IDID = -3 -- The local error test cannot be satisfied
1113 C                   because you specified a zero component in ATOL
1114 C                   and the corresponding computed solution
1115 C                   component is zero. Thus, a pure relative error
1116 C                   test is impossible for this component.
1117 C
1118 C           IDID = -6 -- DDASSL had repeated error test
1119 C                   failures on the last attempted step.
1120 C
1121 C           IDID = -7 -- The corrector could not converge.
1122 C
1123 C           IDID = -8 -- The matrix of partial derivatives
1124 C                   is singular.
1125 C
1126 C           IDID = -9 -- The corrector could not converge.
1127 C                   there were repeated error test failures
1128 C                   in this step.
1129 C
1130 C           IDID =-10 -- The corrector could not converge
1131 C                   because IRES was equal to minus one.
1132 C
1133 C           IDID =-11 -- IRES equal to -2 was encountered
1134 C                   and control is being returned to the
1135 C                   calling program.
1136 C
1137 C           IDID =-12 -- DDASSL failed to compute the initial
1138 C                   YPRIME.
1139 C
1140 C
1141 C
1142 C           IDID = -13,..,-32 -- Not applicable for this code
1143 C
1144 C                    *** Task terminated ***
1145 C                Reported by the value of IDID=-33
1146 C
1147 C           IDID = -33 -- The code has encountered trouble from which
1148 C                   it cannot recover. A message is printed
1149 C                   explaining the trouble and control is returned
1150 C                   to the calling program. For example, this occurs
1151 C                   when invalid input is detected.
1152 C
1153 C  RTOL, ATOL -- These quantities remain unchanged except when
1154 C               IDID = -2. In this case, the error tolerances have been
1155 C               increased by the code to values which are estimated to
1156 C               be appropriate for continuing the integration. However,
1157 C               the reported solution at T was obtained using the input
1158 C               values of RTOL and ATOL.
1159 C
1160 C  RWORK, IWORK -- Contain information which is usually of no
1161 C               interest to the user but necessary for subsequent calls.
1162 C               However, you may find use for
1163 C
1164 C               RWORK(3)--Which contains the step size H to be
1165 C                       attempted on the next step.
1166 C
1167 C               RWORK(4)--Which contains the current value of the
1168 C                       independent variable, i.e., the farthest point
1169 C                       integration has reached. This will be different
1170 C                       from T only when interpolation has been
1171 C                       performed (IDID=3).
1172 C
1173 C               RWORK(7)--Which contains the stepsize used
1174 C                       on the last successful step.
1175 C
1176 C               IWORK(7)--Which contains the order of the method to
1177 C                       be attempted on the next step.
1178 C
1179 C               IWORK(8)--Which contains the order of the method used
1180 C                       on the last step.
1181 C
1182 C               IWORK(11)--Which contains the number of steps taken so
1183 C                        far.
1184 C
1185 C               IWORK(12)--Which contains the number of calls to RES
1186 C                        so far.
1187 C
1188 C               IWORK(13)--Which contains the number of evaluations of
1189 C                        the matrix of partial derivatives needed so
1190 C                        far.
1191 C
1192 C               IWORK(14)--Which contains the total number
1193 C                        of error test failures so far.
1194 C
1195 C               IWORK(15)--Which contains the total number
1196 C                        of convergence test failures so far.
1197 C                        (includes singular iteration matrix
1198 C                        failures.)
1199 C
1200 C
1201 C  -------- INPUT -- WHAT TO DO TO CONTINUE THE INTEGRATION ------------
1202 C                    (CALLS AFTER THE FIRST)
1203 C
1204 C  This code is organized so that subsequent calls to continue the
1205 C  integration involve little (if any) additional effort on your
1206 C  part. You must monitor the IDID parameter in order to determine
1207 C  what to do next.
1208 C
1209 C  Recalling that the principal task of the code is to integrate
1210 C  from T to TOUT (the interval mode), usually all you will need
1211 C  to do is specify a new TOUT upon reaching the current TOUT.
1212 C
1213 C  Do not alter any quantity not specifically permitted below,
1214 C  in particular do not alter NEQ,T,Y(*),YPRIME(*),RWORK(*),IWORK(*)
1215 C  or the differential equation in subroutine RES. Any such
1216 C  alteration constitutes a new problem and must be treated as such,
1217 C  i.e., you must start afresh.
1218 C
1219 C  You cannot change from vector to scalar error control or vice
1220 C  versa (INFO(2)), but you can change the size of the entries of
1221 C  RTOL, ATOL. Increasing a tolerance makes the equation easier
1222 C  to integrate. Decreasing a tolerance will make the equation
1223 C  harder to integrate and should generally be avoided.
1224 C
1225 C  You can switch from the intermediate-output mode to the
1226 C  interval mode (INFO(3)) or vice versa at any time.
1227 C
1228 C  If it has been necessary to prevent the integration from going
1229 C  past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the
1230 C  code will not integrate to any TOUT beyond the currently
1231 C  specified TSTOP. Once TSTOP has been reached you must change
1232 C  the value of TSTOP or set INFO(4)=0. You may change INFO(4)
1233 C  or TSTOP at any time but you must supply the value of TSTOP in
1234 C  RWORK(1) whenever you set INFO(4)=1.
1235 C
1236 C  Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2)
1237 C  unless you are going to restart the code.
1238 C
1239 C                 *** Following a completed task ***
1240 C  If
1241 C     IDID = 1, call the code again to continue the integration
1242 C                  another step in the direction of TOUT.
1243 C
1244 C     IDID = 2 or 3, define a new TOUT and call the code again.
1245 C                  TOUT must be different from T. You cannot change
1246 C                  the direction of integration without restarting.
1247 C
1248 C                 *** Following an interrupted task ***
1249 C               To show the code that you realize the task was
1250 C               interrupted and that you want to continue, you
1251 C               must take appropriate action and set INFO(1) = 1
1252 C  If
1253 C    IDID = -1, The code has taken about 500 steps.
1254 C                  If you want to continue, set INFO(1) = 1 and
1255 C                  call the code again. An additional 500 steps
1256 C                  will be allowed.
1257 C
1258 C    IDID = -2, The error tolerances RTOL, ATOL have been
1259 C                  increased to values the code estimates appropriate
1260 C                  for continuing. You may want to change them
1261 C                  yourself. If you are sure you want to continue
1262 C                  with relaxed error tolerances, set INFO(1)=1 and
1263 C                  call the code again.
1264 C
1265 C    IDID = -3, A solution component is zero and you set the
1266 C                  corresponding component of ATOL to zero. If you
1267 C                  are sure you want to continue, you must first
1268 C                  alter the error criterion to use positive values
1269 C                  for those components of ATOL corresponding to zero
1270 C                  solution components, then set INFO(1)=1 and call
1271 C                  the code again.
1272 C
1273 C    IDID = -4,-5  --- Cannot occur with this code.
1274 C
1275 C    IDID = -6, Repeated error test failures occurred on the
1276 C                  last attempted step in DDASSL. A singularity in the
1277 C                  solution may be present. If you are absolutely
1278 C                  certain you want to continue, you should restart
1279 C                  the integration. (Provide initial values of Y and
1280 C                  YPRIME which are consistent)
1281 C
1282 C    IDID = -7, Repeated convergence test failures occurred
1283 C                  on the last attempted step in DDASSL. An inaccurate
1284 C                  or ill-conditioned JACOBIAN may be the problem. If
1285 C                  you are absolutely certain you want to continue, you
1286 C                  should restart the integration.
1287 C
1288 C    IDID = -8, The matrix of partial derivatives is singular.
1289 C                  Some of your equations may be redundant.
1290 C                  DDASSL cannot solve the problem as stated.
1291 C                  It is possible that the redundant equations
1292 C                  could be removed, and then DDASSL could
1293 C                  solve the problem. It is also possible
1294 C                  that a solution to your problem either
1295 C                  does not exist or is not unique.
1296 C
1297 C    IDID = -9, DDASSL had multiple convergence test
1298 C                  failures, preceeded by multiple error
1299 C                  test failures, on the last attempted step.
1300 C                  It is possible that your problem
1301 C                  is ill-posed, and cannot be solved
1302 C                  using this code. Or, there may be a
1303 C                  discontinuity or a singularity in the
1304 C                  solution. If you are absolutely certain
1305 C                  you want to continue, you should restart
1306 C                  the integration.
1307 C
1308 C    IDID =-10, DDASSL had multiple convergence test failures
1309 C                  because IRES was equal to minus one.
1310 C                  If you are absolutely certain you want
1311 C                  to continue, you should restart the
1312 C                  integration.
1313 C
1314 C    IDID =-11, IRES=-2 was encountered, and control is being
1315 C                  returned to the calling program.
1316 C
1317 C    IDID =-12, DDASSL failed to compute the initial YPRIME.
1318 C                  This could happen because the initial
1319 C                  approximation to YPRIME was not very good, or
1320 C                  if a YPRIME consistent with the initial Y
1321 C                  does not exist. The problem could also be caused
1322 C                  by an inaccurate or singular iteration matrix.
1323 C
1324 C    IDID = -13,..,-32  --- Cannot occur with this code.
1325 C
1326 C
1327 C                 *** Following a terminated task ***
1328 C
1329 C  If IDID= -33, you cannot continue the solution of this problem.
1330 C                  An attempt to do so will result in your
1331 C                  run being terminated.
1332 C
1333 C
1334 C  -------- ERROR MESSAGES ---------------------------------------------
1335 C
1336 C      The SLATEC error print routine XERMSG is called in the event of
1337 C   unsuccessful completion of a task.  Most of these are treated as
1338 C   "recoverable errors", which means that (unless the user has directed
1339 C   otherwise) control will be returned to the calling program for
1340 C   possible action after the message has been printed.
1341 C
1342 C   In the event of a negative value of IDID other than -33, an appro-
1343 C   priate message is printed and the "error number" printed by XERMSG
1344 C   is the value of IDID.  There are quite a number of illegal input
1345 C   errors that can lead to a returned value IDID=-33.  The conditions
1346 C   and their printed "error numbers" are as follows:
1347 C
1348 C   Error number       Condition
1349 C
1350 C        1       Some element of INFO vector is not zero or one.
1351 C        2       NEQ .le. 0
1352 C        3       MAXORD not in range.
1353 C        4       LRW is less than the required length for RWORK.
1354 C        5       LIW is less than the required length for IWORK.
1355 C        6       Some element of RTOL is .lt. 0
1356 C        7       Some element of ATOL is .lt. 0
1357 C        8       All elements of RTOL and ATOL are zero.
1358 C        9       INFO(4)=1 and TSTOP is behind TOUT.
1359 C       10       HMAX .lt. 0.0
1360 C       11       TOUT is behind T.
1361 C       12       INFO(8)=1 and H0=0.0
1362 C       13       Some element of WT is .le. 0.0
1363 C       14       TOUT is too close to T to start integration.
1364 C       15       INFO(4)=1 and TSTOP is behind T.
1365 C       16       --( Not used in this version )--
1366 C       17       ML illegal.  Either .lt. 0 or .gt. NEQ
1367 C       18       MU illegal.  Either .lt. 0 or .gt. NEQ
1368 C       19       TOUT = T.
1369 C
1370 C   If DDASSL is called again without any action taken to remove the
1371 C   cause of an unsuccessful return, XERMSG will be called with a fatal
1372 C   error flag, which will cause unconditional termination of the
1373 C   program.  There are two such fatal errors:
1374 C
1375 C   Error number -998:  The last step was terminated with a negative
1376 C       value of IDID other than -33, and no appropriate action was
1377 C       taken.
1378 C
1379 C   Error number -999:  The previous call was terminated because of
1380 C       illegal input (IDID=-33) and there is illegal input in the
1381 C       present call, as well.  (Suspect infinite loop.)
1382 C
1383 C  ---------------------------------------------------------------------
1384 C
1385 C***REFERENCES  A DESCRIPTION OF DASSL: A DIFFERENTIAL/ALGEBRAIC
1386 C                 SYSTEM SOLVER, L. R. PETZOLD, SAND82-8637,
1387 C                 SANDIA NATIONAL LABORATORIES, SEPTEMBER 1982.
1388 C***ROUTINES CALLED  DLAMCH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS,
1389 C                    XERMSG
1390 C***REVISION HISTORY  (YYMMDD)
1391 C   830315  DATE WRITTEN
1392 C   880387  Code changes made.  All common statements have been
1393 C           replaced by a DATA statement, which defines pointers into
1394 C           RWORK, and PARAMETER statements which define pointers
1395 C           into IWORK.  As well the documentation has gone through
1396 C           grammatical changes.
1397 C   881005  The prologue has been changed to mixed case.
1398 C           The subordinate routines had revision dates changed to
1399 C           this date, although the documentation for these routines
1400 C           is all upper case.  No code changes.
1401 C   890511  Code changes made.  The DATA statement in the declaration
1402 C           section of DDASSL was replaced with a PARAMETER
1403 C           statement.  Also the statement S = 100.D0 was removed
1404 C           from the top of the Newton iteration in DDASTP.
1405 C           The subordinate routines had revision dates changed to
1406 C           this date.
1407 C   890517  The revision date syntax was replaced with the revision
1408 C           history syntax.  Also the "DECK" comment was added to
1409 C           the top of all subroutines.  These changes are consistent
1410 C           with new SLATEC guidelines.
1411 C           The subordinate routines had revision dates changed to
1412 C           this date.  No code changes.
1413 C   891013  Code changes made.
1414 C           Removed all occurrances of FLOAT or DBLE.  All operations
1415 C           are now performed with "mixed-mode" arithmetic.
1416 C           Also, specific function names were replaced with generic
1417 C           function names to be consistent with new SLATEC guidelines.
1418 C           In particular:
1419 C              Replaced DSQRT with SQRT everywhere.
1420 C              Replaced DABS with ABS everywhere.
1421 C              Replaced DMIN1 with MIN everywhere.
1422 C              Replaced MIN0 with MIN everywhere.
1423 C              Replaced DMAX1 with MAX everywhere.
1424 C              Replaced MAX0 with MAX everywhere.
1425 C              Replaced DSIGN with SIGN everywhere.
1426 C           Also replaced REVISION DATE with REVISION HISTORY in all
1427 C           subordinate routines.
1428 C  901004  Miscellaneous changes to prologue to complete conversion
1429 C          to SLATEC 4.0 format.  No code changes.  (F.N.Fritsch)
1430 C  901009  Corrected GAMS classification code and converted subsidiary
1431 C          routines to 4.0 format.  No code changes.  (F.N.Fritsch)
1432 C  901010  Converted XERRWV calls to XERMSG calls.  (R.Clemens,AFWL)
1433 C  901019  Code changes made.
1434 C          Merged SLATEC 4.0 changes with previous changes made
1435 C          by C. Ulrich.  Below is a history of the changes made by
1436 C          C. Ulrich. (Changes in subsidiary routines are implied
1437 C          by this history)
1438 C          891228  Bug was found and repaired inside the DDASSL
1439 C                  and DDAINI routines.  DDAINI was incorrectly
1440 C                  returning the initial T with Y and YPRIME
1441 C                  computed at T+H.  The routine now returns T+H
1442 C                  rather than the initial T.
1443 C                  Cosmetic changes made to DDASTP.
1444 C          900904  Three modifications were made to fix a bug (inside
1445 C                  DDASSL) re interpolation for continuation calls and
1446 C                  cases where TN is very close to TSTOP:
1447 C
1448 C                  1) In testing for whether H is too large, just
1449 C                     compare H to (TSTOP - TN), rather than
1450 C                     (TSTOP - TN) * (1-4*UROUND), and set H to
1451 C                     TSTOP - TN.  This will force DDASTP to step
1452 C                     exactly to TSTOP under certain situations
1453 C                     (i.e. when H returned from DDASTP would otherwise
1454 C                     take TN beyond TSTOP).
1455 C
1456 C                  2) Inside the DDASTP loop, interpolate exactly to
1457 C                     TSTOP if TN is very close to TSTOP (rather than
1458 C                     interpolating to within roundoff of TSTOP).
1459 C
1460 C                  3) Modified IDID description for IDID = 2 to say that
1461 C                     the solution is returned by stepping exactly to
1462 C                     TSTOP, rather than TOUT.  (In some cases the
1463 C                     solution is actually obtained by extrapolating
1464 C                     over a distance near unit roundoff to TSTOP,
1465 C                     but this small distance is deemed acceptable in
1466 C                     these circumstances.)
1467 C   901026  Added explicit declarations for all variables and minor
1468 C           cosmetic changes to prologue, removed unreferenced labels,
1469 C           and improved XERMSG calls.  (FNF)
1470 C   901030  Added ERROR MESSAGES section and reworked other sections to
1471 C           be of more uniform format.  (FNF)
1472 C   910624  Fixed minor bug related to HMAX (five lines ending in
1473 C           statement 526 in DDASSL).   (LRP)
1474 C
1475 C***END PROLOGUE  DDASSL
1476 C
1477 C**End
1478 C
1479 C     Declare arguments.
1480 C
1481       INTEGER  NEQ, INFO(15), IDID, LRW, IWORK(*), LIW, IPAR(*)
1482       DOUBLE PRECISION
1483      *   T, Y(*), YPRIME(*), TOUT, RTOL(*), ATOL(*), RWORK(*),
1484      *   RPAR(*)
1485       EXTERNAL  RES, JAC
1486 C
1487 C     Declare externals.
1488 C
1489       EXTERNAL  DLAMCH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS, XERMSG
1490       DOUBLE PRECISION  DLAMCH, DDANRM
1491 C
1492 C     Declare local variables.
1493 C
1494       INTEGER  I, ITEMP, LALPHA, LBETA, LCJ, LCJOLD, LCTF, LDELTA,
1495      *   LENIW, LENPD, LENRW, LE, LETF, LGAMMA, LH, LHMAX, LHOLD, LIPVT,
1496      *   LJCALC, LK, LKOLD, LIWM, LML, LMTYPE, LMU, LMXORD, LNJE, LNPD,
1497      *   LNRE, LNS, LNST, LNSTL, LPD, LPHASE, LPHI, LPSI, LROUND, LS,
1498      *   LSIGMA, LTN, LTSTOP, LWM, LWT, MBAND, MSAVE, MXORD, NPD, NTEMP,
1499      *   NZFLG
1500       DOUBLE PRECISION
1501      *   ATOLI, H, HMAX, HMIN, HO, R, RH, RTOLI, TDIST, TN, TNEXT,
1502      *   TSTOP, UROUND, YPNORM
1503       LOGICAL  DONE
1504 C       Auxiliary variables for conversion of values to be included in
1505 C       error messages.
1506       CHARACTER*8  XERN1, XERN2
1507       CHARACTER*16 XERN3, XERN4
1508 C
1509 C     SET POINTERS INTO IWORK
1510       PARAMETER (LML=1, LMU=2, LMXORD=3, LMTYPE=4, LNST=11,
1511      *  LNRE=12, LNJE=13, LETF=14, LCTF=15, LNPD=16,
1512      *  LIPVT=21, LJCALC=5, LPHASE=6, LK=7, LKOLD=8,
1513      *  LNS=9, LNSTL=10, LIWM=1)
1514 C
1515 C     SET RELATIVE OFFSET INTO RWORK
1516       PARAMETER (NPD=1)
1517 C
1518 C     SET POINTERS INTO RWORK
1519       PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4,
1520      *  LCJ=5, LCJOLD=6, LHOLD=7, LS=8, LROUND=9,
1521      *  LALPHA=11, LBETA=17, LGAMMA=23,
1522      *  LPSI=29, LSIGMA=35, LDELTA=41)
1523 C
1524 C***FIRST EXECUTABLE STATEMENT  DDASSL
1525       IF(INFO(1).NE.0)GO TO 100
1526 C
1527 C-----------------------------------------------------------------------
1528 C     THIS BLOCK IS EXECUTED FOR THE INITIAL CALL ONLY.
1529 C     IT CONTAINS CHECKING OF INPUTS AND INITIALIZATIONS.
1530 C-----------------------------------------------------------------------
1531 C
1532 C     FIRST CHECK INFO ARRAY TO MAKE SURE ALL ELEMENTS OF INFO
1533 C     ARE EITHER ZERO OR ONE.
1534       DO 10 I=2,11
1535          IF(INFO(I).NE.0.AND.INFO(I).NE.1)GO TO 701
1536 10       CONTINUE
1537 C
1538       IF(NEQ.LE.0)GO TO 702
1539 C
1540 C     CHECK AND COMPUTE MAXIMUM ORDER
1541       MXORD=5
1542       IF(INFO(9).EQ.0)GO TO 20
1543          MXORD=IWORK(LMXORD)
1544          IF(MXORD.LT.1.OR.MXORD.GT.5)GO TO 703
1545 20       IWORK(LMXORD)=MXORD
1546 C
1547 C     COMPUTE MTYPE,LENPD,LENRW.CHECK ML AND MU.
1548       IF(INFO(6).NE.0)GO TO 40
1549          LENPD=NEQ**2
1550          LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD
1551          IF(INFO(5).NE.0)GO TO 30
1552             IWORK(LMTYPE)=2
1553             GO TO 60
1554 30          IWORK(LMTYPE)=1
1555             GO TO 60
1556 40    IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717
1557       IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718
1558       LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ
1559       IF(INFO(5).NE.0)GO TO 50
1560          IWORK(LMTYPE)=5
1561          MBAND=IWORK(LML)+IWORK(LMU)+1
1562          MSAVE=(NEQ/MBAND)+1
1563          LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD+2*MSAVE
1564          GO TO 60
1565 50       IWORK(LMTYPE)=4
1566          LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD
1567 C
1568 C     CHECK LENGTHS OF RWORK AND IWORK
1569 60    LENIW=20+NEQ
1570       IWORK(LNPD)=LENPD
1571       IF(LRW.LT.LENRW)GO TO 704
1572       IF(LIW.LT.LENIW)GO TO 705
1573 C
1574 C     CHECK TO SEE THAT TOUT IS DIFFERENT FROM T
1575       IF(TOUT .EQ. T)GO TO 719
1576 C
1577 C     CHECK HMAX
1578       IF(INFO(7).EQ.0)GO TO 70
1579          HMAX=RWORK(LHMAX)
1580          IF(HMAX.LE.0.0D0)GO TO 710
1581 70    CONTINUE
1582 C
1583 C     INITIALIZE COUNTERS
1584       IWORK(LNST)=0
1585       IWORK(LNRE)=0
1586       IWORK(LNJE)=0
1587 C
1588       IWORK(LNSTL)=0
1589       IDID=1
1590       GO TO 200
1591 C
1592 C-----------------------------------------------------------------------
1593 C     THIS BLOCK IS FOR CONTINUATION CALLS
1594 C     ONLY. HERE WE CHECK INFO(1),AND IF THE
1595 C     LAST STEP WAS INTERRUPTED WE CHECK WHETHER
1596 C     APPROPRIATE ACTION WAS TAKEN.
1597 C-----------------------------------------------------------------------
1598 C
1599 100   CONTINUE
1600       IF(INFO(1).EQ.1)GO TO 110
1601       IF(INFO(1).NE.-1)GO TO 701
1602 C
1603 C     IF WE ARE HERE, THE LAST STEP WAS INTERRUPTED
1604 C     BY AN ERROR CONDITION FROM DDASTP,AND
1605 C     APPROPRIATE ACTION WAS NOT TAKEN. THIS
1606 C     IS A FATAL ERROR.
1607       WRITE (XERN1, '(I8)') IDID
1608       CALL XERMSG ('SLATEC', 'DDASSL',
1609      *   'THE LAST STEP TERMINATED WITH A NEGATIVE VALUE OF IDID = ' //
1610      *   XERN1 // ' AND NO APPROPRIATE ACTION WAS TAKEN.  ' //
1611      *   'RUN TERMINATED', -998, 2)
1612       RETURN
1613 110   CONTINUE
1614       IWORK(LNSTL)=IWORK(LNST)
1615 C
1616 C-----------------------------------------------------------------------
1617 C     THIS BLOCK IS EXECUTED ON ALL CALLS.
1618 C     THE ERROR TOLERANCE PARAMETERS ARE
1619 C     CHECKED, AND THE WORK ARRAY POINTERS
1620 C     ARE SET.
1621 C-----------------------------------------------------------------------
1622 C
1623 200   CONTINUE
1624 C     CHECK RTOL,ATOL
1625       NZFLG=0
1626       RTOLI=RTOL(1)
1627       ATOLI=ATOL(1)
1628       DO 210 I=1,NEQ
1629          IF(INFO(2).EQ.1)RTOLI=RTOL(I)
1630          IF(INFO(2).EQ.1)ATOLI=ATOL(I)
1631          IF(RTOLI.GT.0.0D0.OR.ATOLI.GT.0.0D0)NZFLG=1
1632          IF(RTOLI.LT.0.0D0)GO TO 706
1633          IF(ATOLI.LT.0.0D0)GO TO 707
1634 210      CONTINUE
1635       IF(NZFLG.EQ.0)GO TO 708
1636 C
1637 C     SET UP RWORK STORAGE.IWORK STORAGE IS FIXED
1638 C     IN DATA STATEMENT.
1639       LE=LDELTA+NEQ
1640       LWT=LE+NEQ
1641       LPHI=LWT+NEQ
1642       LPD=LPHI+(IWORK(LMXORD)+1)*NEQ
1643       LWM=LPD
1644       NTEMP=NPD+IWORK(LNPD)
1645       IF(INFO(1).EQ.1)GO TO 400
1646 C
1647 C-----------------------------------------------------------------------
1648 C     THIS BLOCK IS EXECUTED ON THE INITIAL CALL
1649 C     ONLY. SET THE INITIAL STEP SIZE, AND
1650 C     THE ERROR WEIGHT VECTOR, AND PHI.
1651 C     COMPUTE INITIAL YPRIME, IF NECESSARY.
1652 C-----------------------------------------------------------------------
1653 C
1654       TN=T
1655       IDID=1
1656 C
1657 C     SET ERROR WEIGHT VECTOR WT
1658       CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
1659       DO 305 I = 1,NEQ
1660          IF(RWORK(LWT+I-1).LE.0.0D0) GO TO 713
1661 305      CONTINUE
1662 C
1663 C     COMPUTE UNIT ROUNDOFF AND HMIN
1664       UROUND = DLAMCH('P')
1665       RWORK(LROUND) = UROUND
1666       HMIN = 4.0D0*UROUND*MAX(ABS(T),ABS(TOUT))
1667 C
1668 C     CHECK INITIAL INTERVAL TO SEE THAT IT IS LONG ENOUGH
1669       TDIST = ABS(TOUT - T)
1670       IF(TDIST .LT. HMIN) GO TO 714
1671 C
1672 C     CHECK HO, IF THIS WAS INPUT
1673       IF (INFO(8) .EQ. 0) GO TO 310
1674          HO = RWORK(LH)
1675          IF ((TOUT - T)*HO .LT. 0.0D0) GO TO 711
1676          IF (HO .EQ. 0.0D0) GO TO 712
1677          GO TO 320
1678 310    CONTINUE
1679 C
1680 C     COMPUTE INITIAL STEPSIZE, TO BE USED BY EITHER
1681 C     DDASTP OR DDAINI, DEPENDING ON INFO(11)
1682       HO = 0.001D0*TDIST
1683       YPNORM = DDANRM(NEQ,YPRIME,RWORK(LWT),RPAR,IPAR)
1684       IF (YPNORM .GT. 0.5D0/HO) HO = 0.5D0/YPNORM
1685       HO = SIGN(HO,TOUT-T)
1686 C     ADJUST HO IF NECESSARY TO MEET HMAX BOUND
1687 320   IF (INFO(7) .EQ. 0) GO TO 330
1688          RH = ABS(HO)/RWORK(LHMAX)
1689          IF (RH .GT. 1.0D0) HO = HO/RH
1690 C     COMPUTE TSTOP, IF APPLICABLE
1691 330   IF (INFO(4) .EQ. 0) GO TO 340
1692          TSTOP = RWORK(LTSTOP)
1693          IF ((TSTOP - T)*HO .LT. 0.0D0) GO TO 715
1694          IF ((T + HO - TSTOP)*HO .GT. 0.0D0) HO = TSTOP - T
1695          IF ((TSTOP - TOUT)*HO .LT. 0.0D0) GO TO 709
1696 C
1697 C     COMPUTE INITIAL DERIVATIVE, UPDATING TN AND Y, IF APPLICABLE
1698 340   IF (INFO(11) .EQ. 0) GO TO 350
1699       CALL DDAINI(TN,Y,YPRIME,NEQ,
1700      *  RES,JAC,HO,RWORK(LWT),IDID,RPAR,IPAR,
1701      *  RWORK(LPHI),RWORK(LDELTA),RWORK(LE),
1702      *  RWORK(LWM),IWORK(LIWM),HMIN,RWORK(LROUND),
1703      *  INFO(10),NTEMP)
1704       if(iero.ne.0) return
1705       IF (IDID .LT. 0) GO TO 390
1706 C
1707 C     LOAD H WITH HO.  STORE H IN RWORK(LH)
1708 350   H = HO
1709       RWORK(LH) = H
1710 C
1711 C     LOAD Y AND H*YPRIME INTO PHI(*,1) AND PHI(*,2)
1712       ITEMP = LPHI + NEQ
1713       DO 370 I = 1,NEQ
1714          RWORK(LPHI + I - 1) = Y(I)
1715 370      RWORK(ITEMP + I - 1) = H*YPRIME(I)
1716 C
1717 390   GO TO 500
1718 C
1719 C-------------------------------------------------------
1720 C     THIS BLOCK IS FOR CONTINUATION CALLS ONLY. ITS
1721 C     PURPOSE IS TO CHECK STOP CONDITIONS BEFORE
1722 C     TAKING A STEP.
1723 C     ADJUST H IF NECESSARY TO MEET HMAX BOUND
1724 C-------------------------------------------------------
1725 C
1726 400   CONTINUE
1727       UROUND=RWORK(LROUND)
1728       DONE = .FALSE.
1729       TN=RWORK(LTN)
1730       H=RWORK(LH)
1731       IF(INFO(7) .EQ. 0) GO TO 410
1732          RH = ABS(H)/RWORK(LHMAX)
1733          IF(RH .GT. 1.0D0) H = H/RH
1734 410   CONTINUE
1735       IF(T .EQ. TOUT) GO TO 719
1736       IF((T - TOUT)*H .GT. 0.0D0) GO TO 711
1737       IF(INFO(4) .EQ. 1) GO TO 430
1738       IF(INFO(3) .EQ. 1) GO TO 420
1739       IF((TN-TOUT)*H.LT.0.0D0)GO TO 490
1740       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1741      *  RWORK(LPHI),RWORK(LPSI))
1742       T=TOUT
1743       IDID = 3
1744       DONE = .TRUE.
1745       GO TO 490
1746 420   IF((TN-T)*H .LE. 0.0D0) GO TO 490
1747       IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425
1748       CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
1749      *  RWORK(LPHI),RWORK(LPSI))
1750       T = TN
1751       IDID = 1
1752       DONE = .TRUE.
1753       GO TO 490
1754 425   CONTINUE
1755       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1756      *  RWORK(LPHI),RWORK(LPSI))
1757       T = TOUT
1758       IDID = 3
1759       DONE = .TRUE.
1760       GO TO 490
1761 430   IF(INFO(3) .EQ. 1) GO TO 440
1762       TSTOP=RWORK(LTSTOP)
1763       IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715
1764       IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709
1765       IF((TN-TOUT)*H.LT.0.0D0)GO TO 450
1766       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1767      *   RWORK(LPHI),RWORK(LPSI))
1768       T=TOUT
1769       IDID = 3
1770       DONE = .TRUE.
1771       GO TO 490
1772 440   TSTOP = RWORK(LTSTOP)
1773       IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715
1774       IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709
1775       IF((TN-T)*H .LE. 0.0D0) GO TO 450
1776       IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445
1777       CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
1778      *  RWORK(LPHI),RWORK(LPSI))
1779       T = TN
1780       IDID = 1
1781       DONE = .TRUE.
1782       GO TO 490
1783 445   CONTINUE
1784       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1785      *  RWORK(LPHI),RWORK(LPSI))
1786       T = TOUT
1787       IDID = 3
1788       DONE = .TRUE.
1789       GO TO 490
1790 450   CONTINUE
1791 C     CHECK WHETHER WE ARE WITHIN ROUNDOFF OF TSTOP
1792       IF(ABS(TN-TSTOP).GT.100.0D0*UROUND*
1793      *   (ABS(TN)+ABS(H)))GO TO 460
1794       CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD),
1795      *  RWORK(LPHI),RWORK(LPSI))
1796       IDID=2
1797       T=TSTOP
1798       DONE = .TRUE.
1799       GO TO 490
1800 460   TNEXT=TN+H
1801       IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490
1802       H=TSTOP-TN
1803       RWORK(LH)=H
1804 C
1805 490   IF (DONE) GO TO 580
1806 C
1807 C-------------------------------------------------------
1808 C     THE NEXT BLOCK CONTAINS THE CALL TO THE
1809 C     ONE-STEP INTEGRATOR DDASTP.
1810 C     THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS.
1811 C     CHECK FOR TOO MANY STEPS.
1812 C     UPDATE WT.
1813 C     CHECK FOR TOO MUCH ACCURACY REQUESTED.
1814 C     COMPUTE MINIMUM STEPSIZE.
1815 C-------------------------------------------------------
1816 C
1817 500   CONTINUE
1818 C     CHECK FOR FAILURE TO COMPUTE INITIAL YPRIME
1819       IF (IDID .EQ. -12) GO TO 527
1820 C
1821 C     CHECK FOR TOO MANY STEPS
1822       IF((IWORK(LNST)-IWORK(LNSTL)).LT.500)
1823      *   GO TO 510
1824            IDID=-1
1825            GO TO 527
1826 C
1827 C     UPDATE WT
1828 510   CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI),
1829      *  RWORK(LWT),RPAR,IPAR)
1830       DO 520 I=1,NEQ
1831          IF(RWORK(I+LWT-1).GT.0.0D0)GO TO 520
1832            IDID=-3
1833            GO TO 527
1834 520   CONTINUE
1835 C
1836 C     TEST FOR TOO MUCH ACCURACY REQUESTED.
1837       R=DDANRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)*
1838      *   100.0D0*UROUND
1839       IF(R.LE.1.0D0)GO TO 525
1840 C     MULTIPLY RTOL AND ATOL BY R AND RETURN
1841       IF(INFO(2).EQ.1)GO TO 523
1842            RTOL(1)=R*RTOL(1)
1843            ATOL(1)=R*ATOL(1)
1844            IDID=-2
1845            GO TO 527
1846 523   DO 524 I=1,NEQ
1847            RTOL(I)=R*RTOL(I)
1848 524        ATOL(I)=R*ATOL(I)
1849       IDID=-2
1850       GO TO 527
1851 525   CONTINUE
1852 C
1853 C     COMPUTE MINIMUM STEPSIZE
1854       HMIN=4.0D0*UROUND*MAX(ABS(TN),ABS(TOUT))
1855 C
1856 C     TEST H VS. HMAX
1857       IF (INFO(7) .EQ. 0) GO TO 526
1858          RH = ABS(H)/RWORK(LHMAX)
1859          IF (RH .GT. 1.0D0) H = H/RH
1860 526   CONTINUE           
1861 C
1862       CALL DDASTP(TN,Y,YPRIME,NEQ,
1863      *   RES,JAC,H,RWORK(LWT),INFO(1),IDID,RPAR,IPAR,
1864      *   RWORK(LPHI),RWORK(LDELTA),RWORK(LE),
1865      *   RWORK(LWM),IWORK(LIWM),
1866      *   RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA),
1867      *   RWORK(LPSI),RWORK(LSIGMA),
1868      *   RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),
1869      *   RWORK(LS),HMIN,RWORK(LROUND),
1870      *   IWORK(LPHASE),IWORK(LJCALC),IWORK(LK),
1871      *   IWORK(LKOLD),IWORK(LNS),INFO(10),NTEMP)
1872       if(iero.ne.0) return
1873 527   IF(IDID.LT.0)GO TO 600
1874 C
1875 C--------------------------------------------------------
1876 C     THIS BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN
1877 C     FROM DDASTP (IDID=1).  TEST FOR STOP CONDITIONS.
1878 C--------------------------------------------------------
1879 C
1880       IF(INFO(4).NE.0)GO TO 540
1881            IF(INFO(3).NE.0)GO TO 530
1882              IF((TN-TOUT)*H.LT.0.0D0)GO TO 500
1883              CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1884      *         IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1885              IDID=3
1886              T=TOUT
1887              GO TO 580
1888 530          IF((TN-TOUT)*H.GE.0.0D0)GO TO 535
1889              T=TN
1890              IDID=1
1891              GO TO 580
1892 535          CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1893      *         IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1894              IDID=3
1895              T=TOUT
1896              GO TO 580
1897 540   IF(INFO(3).NE.0)GO TO 550
1898       IF((TN-TOUT)*H.LT.0.0D0)GO TO 542
1899          CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1900      *     IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1901          T=TOUT
1902          IDID=3
1903          GO TO 580
1904 542   IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*
1905      *   (ABS(TN)+ABS(H)))GO TO 545
1906       TNEXT=TN+H
1907       IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500
1908       H=TSTOP-TN
1909       GO TO 500
1910 545   CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
1911      *  IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1912       IDID=2
1913       T=TSTOP
1914       GO TO 580
1915 550   IF((TN-TOUT)*H.GE.0.0D0)GO TO 555
1916       IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*(ABS(TN)+ABS(H)))GO TO 552
1917       T=TN
1918       IDID=1
1919       GO TO 580
1920 552   CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
1921      *  IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1922       IDID=2
1923       T=TSTOP
1924       GO TO 580
1925 555   CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1926      *   IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1927       T=TOUT
1928       IDID=3
1929       GO TO 580
1930 C
1931 C--------------------------------------------------------
1932 C     ALL SUCCESSFUL RETURNS FROM DDASSL ARE MADE FROM
1933 C     THIS BLOCK.
1934 C--------------------------------------------------------
1935 C
1936 580   CONTINUE
1937       RWORK(LTN)=TN
1938       RWORK(LH)=H
1939       RETURN
1940 C
1941 C-----------------------------------------------------------------------
1942 C     THIS BLOCK HANDLES ALL UNSUCCESSFUL
1943 C     RETURNS OTHER THAN FOR ILLEGAL INPUT.
1944 C-----------------------------------------------------------------------
1945 C
1946 600   CONTINUE
1947       ITEMP=-IDID
1948       GO TO (610,620,630,690,690,640,650,660,670,675,
1949      *  680,685), ITEMP
1950 C
1951 C     THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE
1952 C     REACHING TOUT
1953 610   WRITE (XERN3, '(1P,D15.6)') TN
1954       CALL XERMSG ('SLATEC', 'DDASSL',
1955      *   'AT CURRENT T = ' // XERN3 // ' 500 STEPS TAKEN ON THIS ' //
1956      *   'CALL BEFORE REACHING TOUT', IDID, 1)
1957       GO TO 690
1958 C
1959 C     TOO MUCH ACCURACY FOR MACHINE PRECISION
1960 620   WRITE (XERN3, '(1P,D15.6)') TN
1961       CALL XERMSG ('SLATEC', 'DDASSL',
1962      *   'AT T = ' // XERN3 // ' TOO MUCH ACCURACY REQUESTED FOR ' //
1963      *   'PRECISION OF MACHINE. RTOL AND ATOL WERE INCREASED TO ' //
1964      *   'APPROPRIATE VALUES', IDID, 1)
1965       GO TO 690
1966 C
1967 C     WT(I) .LE. 0.0 FOR SOME I (NOT AT START OF PROBLEM)
1968 630   WRITE (XERN3, '(1P,D15.6)') TN
1969       CALL XERMSG ('SLATEC', 'DDASSL',
1970      *   'AT T = ' // XERN3 // ' SOME ELEMENT OF WT HAS BECOME .LE. ' //
1971      *   '0.0', IDID, 1)
1972       GO TO 690
1973 C
1974 C     ERROR TEST FAILED REPEATEDLY OR WITH H=HMIN
1975 640   WRITE (XERN3, '(1P,D15.6)') TN
1976       WRITE (XERN4, '(1P,D15.6)') H
1977       CALL XERMSG ('SLATEC', 'DDASSL',
1978      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
1979      *   ' THE ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN',
1980      *   IDID, 1)
1981       GO TO 690
1982 C
1983 C     CORRECTOR CONVERGENCE FAILED REPEATEDLY OR WITH H=HMIN
1984 650   WRITE (XERN3, '(1P,D15.6)') TN
1985       WRITE (XERN4, '(1P,D15.6)') H
1986       CALL XERMSG ('SLATEC', 'DDASSL',
1987      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
1988      *   ' THE CORRECTOR FAILED TO CONVERGE REPEATEDLY OR WITH ' //
1989      *   'ABS(H)=HMIN', IDID, 1)
1990       GO TO 690
1991 C
1992 C     THE ITERATION MATRIX IS SINGULAR
1993 660   WRITE (XERN3, '(1P,D15.6)') TN
1994       WRITE (XERN4, '(1P,D15.6)') H
1995       CALL XERMSG ('SLATEC', 'DDASSL',
1996      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
1997      *   ' THE ITERATION MATRIX IS SINGULAR', IDID, 1)
1998       GO TO 690
1999 C
2000 C     CORRECTOR FAILURE PRECEEDED BY ERROR TEST FAILURES.
2001 670   WRITE (XERN3, '(1P,D15.6)') TN
2002       WRITE (XERN4, '(1P,D15.6)') H
2003       CALL XERMSG ('SLATEC', 'DDASSL',
2004      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
2005      *   ' THE CORRECTOR COULD NOT CONVERGE.  ALSO, THE ERROR TEST ' //
2006      *   'FAILED REPEATEDLY.', IDID, 1)
2007       GO TO 690
2008 C
2009 C     CORRECTOR FAILURE BECAUSE IRES = -1
2010 675   WRITE (XERN3, '(1P,D15.6)') TN
2011       WRITE (XERN4, '(1P,D15.6)') H
2012       CALL XERMSG ('SLATEC', 'DDASSL',
2013      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
2014      *   ' THE CORRECTOR COULD NOT CONVERGE BECAUSE IRES WAS EQUAL ' //
2015      *   'TO MINUS ONE', IDID, 1)
2016       GO TO 690
2017 C
2018 C     FAILURE BECAUSE IRES = -2
2019 680   WRITE (XERN3, '(1P,D15.6)') TN
2020       WRITE (XERN4, '(1P,D15.6)') H
2021       CALL XERMSG ('SLATEC', 'DDASSL',
2022      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
2023      *   ' IRES WAS EQUAL TO MINUS TWO', IDID, 1)
2024       GO TO 690
2025 C
2026 C     FAILED TO COMPUTE INITIAL YPRIME
2027 685   WRITE (XERN3, '(1P,D15.6)') TN
2028       WRITE (XERN4, '(1P,D15.6)') HO
2029       CALL XERMSG ('SLATEC', 'DDASSL',
2030      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
2031      *   ' THE INITIAL YPRIME COULD NOT BE COMPUTED', IDID, 1)
2032       GO TO 690
2033 C
2034 690   CONTINUE
2035       INFO(1)=-1
2036       T=TN
2037       RWORK(LTN)=TN
2038       RWORK(LH)=H
2039       RETURN
2040 C
2041 C-----------------------------------------------------------------------
2042 C     THIS BLOCK HANDLES ALL ERROR RETURNS DUE
2043 C     TO ILLEGAL INPUT, AS DETECTED BEFORE CALLING
2044 C     DDASTP. FIRST THE ERROR MESSAGE ROUTINE IS
2045 C     CALLED. IF THIS HAPPENS TWICE IN
2046 C     SUCCESSION, EXECUTION IS TERMINATED
2047 C
2048 C-----------------------------------------------------------------------
2049 701   CALL XERMSG ('SLATEC', 'DDASSL',
2050      *   'SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE', 1, 1)
2051       GO TO 750
2052 C
2053 702   WRITE (XERN1, '(I8)') NEQ
2054       CALL XERMSG ('SLATEC', 'DDASSL',
2055      *   'NEQ = ' // XERN1 // ' .LE. 0', 2, 1)
2056       GO TO 750
2057 C
2058 703   WRITE (XERN1, '(I8)') MXORD
2059       CALL XERMSG ('SLATEC', 'DDASSL',
2060      *   'MAXORD = ' // XERN1 // ' NOT IN RANGE', 3, 1)
2061       GO TO 750
2062 C
2063 704   WRITE (XERN1, '(I8)') LENRW
2064       WRITE (XERN2, '(I8)') LRW
2065       CALL XERMSG ('SLATEC', 'DDASSL',
2066      *   'RWORK LENGTH NEEDED, LENRW = ' // XERN1 //
2067      *   ', EXCEEDS LRW = ' // XERN2, 4, 1)
2068       GO TO 750
2069 C
2070 705   WRITE (XERN1, '(I8)') LENIW
2071       WRITE (XERN2, '(I8)') LIW
2072       CALL XERMSG ('SLATEC', 'DDASSL',
2073      *   'IWORK LENGTH NEEDED, LENIW = ' // XERN1 //
2074      *   ', EXCEEDS LIW = ' // XERN2, 5, 1)
2075       GO TO 750
2076 C
2077 706   CALL XERMSG ('SLATEC', 'DDASSL',
2078      *   'SOME ELEMENT OF RTOL IS .LT. 0', 6, 1)
2079       GO TO 750
2080 C
2081 707   CALL XERMSG ('SLATEC', 'DDASSL',
2082      *   'SOME ELEMENT OF ATOL IS .LT. 0', 7, 1)
2083       GO TO 750
2084 C
2085 708   CALL XERMSG ('SLATEC', 'DDASSL',
2086      *   'ALL ELEMENTS OF RTOL AND ATOL ARE ZERO', 8, 1)
2087       GO TO 750
2088 C
2089 709   WRITE (XERN3, '(1P,D15.6)') TSTOP
2090       WRITE (XERN4, '(1P,D15.6)') TOUT
2091       CALL XERMSG ('SLATEC', 'DDASSL',
2092      *   'INFO(4) = 1 AND TSTOP = ' // XERN3 // ' BEHIND TOUT = ' //
2093      *   XERN4, 9, 1)
2094       GO TO 750
2095 C
2096 710   WRITE (XERN3, '(1P,D15.6)') HMAX
2097       CALL XERMSG ('SLATEC', 'DDASSL',
2098      *   'HMAX = ' // XERN3 // ' .LT. 0.0', 10, 1)
2099       GO TO 750
2100 C
2101 711   WRITE (XERN3, '(1P,D15.6)') TOUT
2102       WRITE (XERN4, '(1P,D15.6)') T
2103       CALL XERMSG ('SLATEC', 'DDASSL',
2104      *   'TOUT = ' // XERN3 // ' BEHIND T = ' // XERN4, 11, 1)
2105       GO TO 750
2106 C
2107 712   CALL XERMSG ('SLATEC', 'DDASSL',
2108      *   'INFO(8)=1 AND H0=0.0', 12, 1)
2109       GO TO 750
2110 C
2111 713   CALL XERMSG ('SLATEC', 'DDASSL',
2112      *   'SOME ELEMENT OF WT IS .LE. 0.0', 13, 1)
2113       GO TO 750
2114 C
2115 714   WRITE (XERN3, '(1P,D15.6)') TOUT
2116       WRITE (XERN4, '(1P,D15.6)') T
2117       CALL XERMSG ('SLATEC', 'DDASSL',
2118      *   'TOUT = ' // XERN3 // ' TOO CLOSE TO T = ' // XERN4 //
2119      *   ' TO START INTEGRATION', 14, 1)
2120       GO TO 750
2121 C
2122 715   WRITE (XERN3, '(1P,D15.6)') TSTOP
2123       WRITE (XERN4, '(1P,D15.6)') T
2124       CALL XERMSG ('SLATEC', 'DDASSL',
2125      *   'INFO(4)=1 AND TSTOP = ' // XERN3 // ' BEHIND T = ' // XERN4,
2126      *   15, 1)
2127       GO TO 750
2128 C
2129 717   WRITE (XERN1, '(I8)') IWORK(LML)
2130       CALL XERMSG ('SLATEC', 'DDASSL',
2131      *   'ML = ' // XERN1 // ' ILLEGAL.  EITHER .LT. 0 OR .GT. NEQ',
2132      *   17, 1)
2133       GO TO 750
2134 C
2135 718   WRITE (XERN1, '(I8)') IWORK(LMU)
2136       CALL XERMSG ('SLATEC', 'DDASSL',
2137      *   'MU = ' // XERN1 // ' ILLEGAL.  EITHER .LT. 0 OR .GT. NEQ',
2138      *   18, 1)
2139       GO TO 750
2140 C
2141 719   WRITE (XERN3, '(1P,D15.6)') TOUT
2142       CALL XERMSG ('SLATEC', 'DDASSL',
2143      *  'TOUT = T = ' // XERN3, 19, 1)
2144       GO TO 750
2145 C
2146 750   IDID=-33
2147       IF(INFO(1).EQ.-1) THEN
2148          CALL XERMSG ('SLATEC', 'DDASSL',
2149      *      'REPEATED OCCURRENCES OF ILLEGAL INPUT$$' //
2150      *      'RUN TERMINATED. APPARENT INFINITE LOOP', -999, 2)
2151       ENDIF
2152 C
2153       INFO(1)=-1
2154       RETURN
2155 C-----------END OF SUBROUTINE DDASSL------------------------------------
2156       END
2157
2158       SUBROUTINE DDASTP (X, Y, YPRIME, NEQ, RES, JAC, H, WT, JSTART,
2159      +   IDID, RPAR, IPAR, PHI, DELTA, E, WM, IWM, ALPHA, BETA, GAMMA,
2160      +   PSI, SIGMA, CJ, CJOLD, HOLD, S, HMIN, UROUND, IPHASE, JCALC,
2161      +   K, KOLD, NS, NONNEG, NTEMP)
2162       
2163       include 'stack.h'
2164
2165 C***BEGIN PROLOGUE  DDASTP
2166 C***SUBSIDIARY
2167 C***PURPOSE  Perform one step of the DDASSL integration.
2168 C***LIBRARY   SLATEC (DASSL)
2169 C***TYPE      DOUBLE PRECISION (SDASTP-S, DDASTP-D)
2170 C***AUTHOR  PETZOLD, LINDA R., (LLNL)
2171 C***DESCRIPTION
2172 C-----------------------------------------------------------------------
2173 C     DDASTP SOLVES A SYSTEM OF DIFFERENTIAL/
2174 C     ALGEBRAIC EQUATIONS OF THE FORM
2175 C     G(X,Y,YPRIME) = 0,  FOR ONE STEP (NORMALLY
2176 C     FROM X TO X+H).
2177 C
2178 C     THE METHODS USED ARE MODIFIED DIVIDED
2179 C     DIFFERENCE,FIXED LEADING COEFFICIENT
2180 C     FORMS OF BACKWARD DIFFERENTIATION
2181 C     FORMULAS. THE CODE ADJUSTS THE STEPSIZE
2182 C     AND ORDER TO CONTROL THE LOCAL ERROR PER
2183 C     STEP.
2184 C
2185 C
2186 C     THE PARAMETERS REPRESENT
2187 C     X  --        INDEPENDENT VARIABLE
2188 C     Y  --        SOLUTION VECTOR AT X
2189 C     YPRIME --    DERIVATIVE OF SOLUTION VECTOR
2190 C                  AFTER SUCCESSFUL STEP
2191 C     NEQ --       NUMBER OF EQUATIONS TO BE INTEGRATED
2192 C     RES --       EXTERNAL USER-SUPPLIED SUBROUTINE
2193 C                  TO EVALUATE THE RESIDUAL.  THE CALL IS
2194 C                  CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
2195 C                  X,Y,YPRIME ARE INPUT.  DELTA IS OUTPUT.
2196 C                  ON INPUT, IRES=0.  RES SHOULD ALTER IRES ONLY
2197 C                  IF IT ENCOUNTERS AN ILLEGAL VALUE OF Y OR A
2198 C                  STOP CONDITION.  SET IRES=-1 IF AN INPUT VALUE
2199 C                  OF Y IS ILLEGAL, AND DDASTP WILL TRY TO SOLVE
2200 C                  THE PROBLEM WITHOUT GETTING IRES = -1.  IF
2201 C                  IRES=-2, DDASTP RETURNS CONTROL TO THE CALLING
2202 C                  PROGRAM WITH IDID = -11.
2203 C     JAC --       EXTERNAL USER-SUPPLIED ROUTINE TO EVALUATE
2204 C                  THE ITERATION MATRIX (THIS IS OPTIONAL)
2205 C                  THE CALL IS OF THE FORM
2206 C                  CALL JAC(X,Y,YPRIME,PD,CJ,RPAR,IPAR)
2207 C                  PD IS THE MATRIX OF PARTIAL DERIVATIVES,
2208 C                  PD=DG/DY+CJ*DG/DYPRIME
2209 C     H --         APPROPRIATE STEP SIZE FOR NEXT STEP.
2210 C                  NORMALLY DETERMINED BY THE CODE
2211 C     WT --        VECTOR OF WEIGHTS FOR ERROR CRITERION.
2212 C     JSTART --    INTEGER VARIABLE SET 0 FOR
2213 C                  FIRST STEP, 1 OTHERWISE.
2214 C     IDID --      COMPLETION CODE WITH THE FOLLOWING MEANINGS:
2215 C                  IDID= 1 -- THE STEP WAS COMPLETED SUCCESSFULLY
2216 C                  IDID=-6 -- THE ERROR TEST FAILED REPEATEDLY
2217 C                  IDID=-7 -- THE CORRECTOR COULD NOT CONVERGE
2218 C                  IDID=-8 -- THE ITERATION MATRIX IS SINGULAR
2219 C                  IDID=-9 -- THE CORRECTOR COULD NOT CONVERGE.
2220 C                             THERE WERE REPEATED ERROR TEST
2221 C                             FAILURES ON THIS STEP.
2222 C                  IDID=-10-- THE CORRECTOR COULD NOT CONVERGE
2223 C                             BECAUSE IRES WAS EQUAL TO MINUS ONE
2224 C                  IDID=-11-- IRES EQUAL TO -2 WAS ENCOUNTERED,
2225 C                             AND CONTROL IS BEING RETURNED TO
2226 C                             THE CALLING PROGRAM
2227 C     RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS THAT
2228 C                  ARE USED FOR COMMUNICATION BETWEEN THE
2229 C                  CALLING PROGRAM AND EXTERNAL USER ROUTINES
2230 C                  THEY ARE NOT ALTERED BY DDASTP
2231 C     PHI --       ARRAY OF DIVIDED DIFFERENCES USED BY
2232 C                  DDASTP. THE LENGTH IS NEQ*(K+1),WHERE
2233 C                  K IS THE MAXIMUM ORDER
2234 C     DELTA,E --   WORK VECTORS FOR DDASTP OF LENGTH NEQ
2235 C     WM,IWM --    REAL AND INTEGER ARRAYS STORING
2236 C                  MATRIX INFORMATION SUCH AS THE MATRIX
2237 C                  OF PARTIAL DERIVATIVES,PERMUTATION
2238 C                  VECTOR,AND VARIOUS OTHER INFORMATION.
2239 C
2240 C     THE OTHER PARAMETERS ARE INFORMATION
2241 C     WHICH IS NEEDED INTERNALLY BY DDASTP TO
2242 C     CONTINUE FROM STEP TO STEP.
2243 C
2244 C-----------------------------------------------------------------------
2245 C***ROUTINES CALLED  DDAJAC, DDANRM, DDASLV, DDATRP
2246 C***REVISION HISTORY  (YYMMDD)
2247 C   830315  DATE WRITTEN
2248 C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
2249 C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
2250 C   901026  Added explicit declarations for all variables and minor
2251 C           cosmetic changes to prologue.  (FNF)
2252 C***END PROLOGUE  DDASTP
2253 C
2254       INTEGER  NEQ, JSTART, IDID, IPAR(*), IWM(*), IPHASE, JCALC, K,
2255      *   KOLD, NS, NONNEG, NTEMP
2256       DOUBLE PRECISION
2257      *   X, Y(*), YPRIME(*), H, WT(*), RPAR(*), PHI(NEQ,*), DELTA(*),
2258      *   E(*), WM(*), ALPHA(*), BETA(*), GAMMA(*), PSI(*), SIGMA(*), CJ,
2259      *   CJOLD, HOLD, S, HMIN, UROUND
2260       EXTERNAL  RES, JAC
2261 C
2262       EXTERNAL  DDAJAC, DDANRM, DDASLV, DDATRP
2263       DOUBLE PRECISION  DDANRM
2264 C
2265       INTEGER  I, IER, IRES, J, J1, KDIFF, KM1, KNEW, KP1, KP2, LCTF,
2266      *   LETF, LMXORD, LNJE, LNRE, LNST, M, MAXIT, NCF, NEF, NSF, NSP1
2267       DOUBLE PRECISION
2268      *   ALPHA0, ALPHAS, CJLAST, CK, DELNRM, ENORM, ERK, ERKM1,
2269      *   ERKM2, ERKP1, IERR, EST, HNEW, OLDNRM, PNORM, R, RATE, TEMP1,
2270      *   TEMP2, TERK, TERKM1, TERKM2, TERKP1, XOLD, XRATE
2271       LOGICAL  CONVGD
2272 C
2273       PARAMETER (LMXORD=3)
2274       PARAMETER (LNST=11)
2275       PARAMETER (LNRE=12)
2276       PARAMETER (LNJE=13)
2277       PARAMETER (LETF=14)
2278       PARAMETER (LCTF=15)
2279 C
2280       DATA MAXIT/4/
2281       DATA XRATE/0.25D0/
2282 C
2283 C
2284 C
2285 C
2286 C
2287 C-----------------------------------------------------------------------
2288 C     BLOCK 1.
2289 C     INITIALIZE. ON THE FIRST CALL,SET
2290 C     THE ORDER TO 1 AND INITIALIZE
2291 C     OTHER VARIABLES.
2292 C-----------------------------------------------------------------------
2293 C
2294 C     INITIALIZATIONS FOR ALL CALLS
2295 C***FIRST EXECUTABLE STATEMENT  DDASTP
2296       IDID=1
2297       XOLD=X
2298       NCF=0
2299       NSF=0
2300       NEF=0
2301       IF(JSTART .NE. 0) GO TO 120
2302 C
2303 C     IF THIS IS THE FIRST STEP,PERFORM
2304 C     OTHER INITIALIZATIONS
2305       IWM(LETF) = 0
2306       IWM(LCTF) = 0
2307       K=1
2308       KOLD=0
2309       HOLD=0.0D0
2310       JSTART=1
2311       PSI(1)=H
2312       CJOLD = 1.0D0/H
2313       CJ = CJOLD
2314       S = 100.D0
2315       JCALC = -1
2316       DELNRM=1.0D0
2317       IPHASE = 0
2318       NS=0
2319 120   CONTINUE
2320 C
2321 C
2322 C
2323 C
2324 C
2325 C-----------------------------------------------------------------------
2326 C     BLOCK 2
2327 C     COMPUTE COEFFICIENTS OF FORMULAS FOR
2328 C     THIS STEP.
2329 C-----------------------------------------------------------------------
2330 200   CONTINUE
2331       KP1=K+1
2332       KP2=K+2
2333       KM1=K-1
2334       XOLD=X
2335       IF(H.NE.HOLD.OR.K .NE. KOLD) NS = 0
2336       NS=MIN(NS+1,KOLD+2)
2337       NSP1=NS+1
2338       IF(KP1 .LT. NS)GO TO 230
2339 C
2340       BETA(1)=1.0D0
2341       ALPHA(1)=1.0D0
2342       TEMP1=H
2343       GAMMA(1)=0.0D0
2344       SIGMA(1)=1.0D0
2345       DO 210 I=2,KP1
2346          TEMP2=PSI(I-1)
2347          PSI(I-1)=TEMP1
2348          BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2
2349          TEMP1=TEMP2+H
2350          ALPHA(I)=H/TEMP1
2351          SIGMA(I)=(I-1)*SIGMA(I-1)*ALPHA(I)
2352          GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H
2353 210      CONTINUE
2354       PSI(KP1)=TEMP1
2355 230   CONTINUE
2356 C
2357 C     COMPUTE ALPHAS, ALPHA0
2358       ALPHAS = 0.0D0
2359       ALPHA0 = 0.0D0
2360       DO 240 I = 1,K
2361         ALPHAS = ALPHAS - 1.0D0/I
2362         ALPHA0 = ALPHA0 - ALPHA(I)
2363 240     CONTINUE
2364 C
2365 C     COMPUTE LEADING COEFFICIENT CJ
2366       CJLAST = CJ
2367       CJ = -ALPHAS/H
2368 C
2369 C     COMPUTE VARIABLE STEPSIZE ERROR COEFFICIENT CK
2370       CK = ABS(ALPHA(KP1) + ALPHAS - ALPHA0)
2371       CK = MAX(CK,ALPHA(KP1))
2372 C
2373 C     DECIDE WHETHER NEW JACOBIAN IS NEEDED
2374       TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE)
2375       TEMP2 = 1.0D0/TEMP1
2376       IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1
2377       IF (CJ .NE. CJLAST) S = 100.D0
2378 C
2379 C     CHANGE PHI TO PHI STAR
2380       IF(KP1 .LT. NSP1) GO TO 280
2381       DO 270 J=NSP1,KP1
2382          DO 260 I=1,NEQ
2383 260         PHI(I,J)=BETA(J)*PHI(I,J)
2384 270      CONTINUE
2385 280   CONTINUE
2386 C
2387 C     UPDATE TIME
2388       X=X+H
2389 C
2390 C
2391 C
2392 C
2393 C
2394 C-----------------------------------------------------------------------
2395 C     BLOCK 3
2396 C     PREDICT THE SOLUTION AND DERIVATIVE,
2397 C     AND SOLVE THE CORRECTOR EQUATION
2398 C-----------------------------------------------------------------------
2399 C
2400 C     FIRST,PREDICT THE SOLUTION AND DERIVATIVE
2401 300   CONTINUE
2402       DO 310 I=1,NEQ
2403          Y(I)=PHI(I,1)
2404 310      YPRIME(I)=0.0D0
2405       DO 330 J=2,KP1
2406          DO 320 I=1,NEQ
2407             Y(I)=Y(I)+PHI(I,J)
2408 320         YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J)
2409 330   CONTINUE
2410       PNORM = DDANRM (NEQ,Y,WT,RPAR,IPAR)
2411 C
2412 C
2413 C
2414 C     SOLVE THE CORRECTOR EQUATION USING A
2415 C     MODIFIED NEWTON SCHEME.
2416       CONVGD= .TRUE.
2417       M=0
2418       IWM(LNRE)=IWM(LNRE)+1
2419       IRES = 0
2420       CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
2421       if(iero.ne.0) return
2422       IF (IRES .LT. 0) GO TO 380
2423 C
2424 C
2425 C     IF INDICATED,REEVALUATE THE
2426 C     ITERATION MATRIX PD = DG/DY + CJ*DG/DYPRIME
2427 C     (WHERE G(X,Y,YPRIME)=0). SET
2428 C     JCALC TO 0 AS AN INDICATOR THAT
2429 C     THIS HAS BEEN DONE.
2430       IF(JCALC .NE. -1)GO TO 340
2431       IWM(LNJE)=IWM(LNJE)+1
2432       JCALC=0
2433       CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H,
2434      * IER,WT,E,WM,IWM,RES,IRES,UROUND,JAC,RPAR,
2435      * IPAR,NTEMP)
2436       if(iero.ne.0) return
2437       CJOLD=CJ
2438       S = 100.D0
2439       IF (IRES .LT. 0) GO TO 380
2440       IF(IER .NE. 0)GO TO 380
2441       NSF=0
2442 C
2443 C
2444 C     INITIALIZE THE ERROR ACCUMULATION VECTOR E.
2445 340   CONTINUE
2446       DO 345 I=1,NEQ
2447 345      E(I)=0.0D0
2448 C
2449 C
2450 C     CORRECTOR LOOP.
2451 350   CONTINUE
2452 C
2453 C     MULTIPLY RESIDUAL BY TEMP1 TO ACCELERATE CONVERGENCE
2454       TEMP1 = 2.0D0/(1.0D0 + CJ/CJOLD)
2455       DO 355 I = 1,NEQ
2456 355     DELTA(I) = DELTA(I) * TEMP1
2457 C
2458 C     COMPUTE A NEW ITERATE (BACK-SUBSTITUTION).
2459 C     STORE THE CORRECTION IN DELTA.
2460       CALL DDASLV(NEQ,DELTA,WM,IWM)
2461 C
2462 C     UPDATE Y,E,AND YPRIME
2463       DO 360 I=1,NEQ
2464          Y(I)=Y(I)-DELTA(I)
2465          E(I)=E(I)-DELTA(I)
2466 360      YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
2467 C
2468 C     TEST FOR CONVERGENCE OF THE ITERATION
2469       DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
2470       IF (DELNRM .LE. 100.D0*UROUND*PNORM) GO TO 375
2471       IF (M .GT. 0) GO TO 365
2472          OLDNRM = DELNRM
2473          GO TO 367
2474 365   RATE = (DELNRM/OLDNRM)**(1.0D0/M)
2475       IF (RATE .GT. 0.90D0) GO TO 370
2476       S = RATE/(1.0D0 - RATE)
2477 367   IF (S*DELNRM .LE. 0.33D0) GO TO 375
2478 C
2479 C     THE CORRECTOR HAS NOT YET CONVERGED.
2480 C     UPDATE M AND TEST WHETHER THE
2481 C     MAXIMUM NUMBER OF ITERATIONS HAVE
2482 C     BEEN TRIED.
2483       M=M+1
2484       IF(M.GE.MAXIT)GO TO 370
2485 C
2486 C     EVALUATE THE RESIDUAL
2487 C     AND GO BACK TO DO ANOTHER ITERATION
2488       IWM(LNRE)=IWM(LNRE)+1
2489       IRES = 0
2490       CALL RES(X,Y,YPRIME,DELTA,IRES,
2491      *  RPAR,IPAR)
2492       if(iero.ne.0) return
2493       IF (IRES .LT. 0) GO TO 380
2494       GO TO 350
2495 C
2496 C
2497 C     THE CORRECTOR FAILED TO CONVERGE IN MAXIT
2498 C     ITERATIONS. IF THE ITERATION MATRIX
2499 C     IS NOT CURRENT,RE-DO THE STEP WITH
2500 C     A NEW ITERATION MATRIX.
2501 370   CONTINUE
2502       IF(JCALC.EQ.0)GO TO 380
2503       JCALC=-1
2504       GO TO 300
2505 C
2506 C
2507 C     THE ITERATION HAS CONVERGED.  IF NONNEGATIVITY OF SOLUTION IS
2508 C     REQUIRED, SET THE SOLUTION NONNEGATIVE, IF THE PERTURBATION
2509 C     TO DO IT IS SMALL ENOUGH.  IF THE CHANGE IS TOO LARGE, THEN
2510 C     CONSIDER THE CORRECTOR ITERATION TO HAVE FAILED.
2511 375   IF(NONNEG .EQ. 0) GO TO 390
2512       DO 377 I = 1,NEQ
2513 377      DELTA(I) = MIN(Y(I),0.0D0)
2514       DELNRM = DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
2515       IF(DELNRM .GT. 0.33D0) GO TO 380
2516       DO 378 I = 1,NEQ
2517 378      E(I) = E(I) - DELTA(I)
2518       GO TO 390
2519 C
2520 C
2521 C     EXITS FROM BLOCK 3
2522 C     NO CONVERGENCE WITH CURRENT ITERATION
2523 C     MATRIX,OR SINGULAR ITERATION MATRIX
2524 380   CONVGD= .FALSE.
2525 390   JCALC = 1
2526       IF(.NOT.CONVGD)GO TO 600
2527 C
2528 C
2529 C
2530 C
2531 C
2532 C-----------------------------------------------------------------------
2533 C     BLOCK 4
2534 C     ESTIMATE THE ERRORS AT ORDERS K,K-1,K-2
2535 C     AS IF CONSTANT STEPSIZE WAS USED. ESTIMATE
2536 C     THE LOCAL ERROR AT ORDER K AND TEST
2537 C     WHETHER THE CURRENT STEP IS SUCCESSFUL.
2538 C-----------------------------------------------------------------------
2539 C
2540 C     ESTIMATE ERRORS AT ORDERS K,K-1,K-2
2541       ENORM = DDANRM(NEQ,E,WT,RPAR,IPAR)
2542       ERK = SIGMA(K+1)*ENORM
2543       TERK = (K+1)*ERK
2544       EST = ERK
2545       KNEW=K
2546       IF(K .EQ. 1)GO TO 430
2547       DO 405 I = 1,NEQ
2548 405     DELTA(I) = PHI(I,KP1) + E(I)
2549       ERKM1=SIGMA(K)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
2550       TERKM1 = K*ERKM1
2551       IF(K .GT. 2)GO TO 410
2552       IF(TERKM1 .LE. 0.5D0*TERK)GO TO 420
2553       GO TO 430
2554 410   CONTINUE
2555       DO 415 I = 1,NEQ
2556 415     DELTA(I) = PHI(I,K) + DELTA(I)
2557       ERKM2=SIGMA(K-1)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
2558       TERKM2 = (K-1)*ERKM2
2559       IF(MAX(TERKM1,TERKM2).GT.TERK)GO TO 430
2560 C     LOWER THE ORDER
2561 420   CONTINUE
2562       KNEW=K-1
2563       EST = ERKM1
2564 C
2565 C
2566 C     CALCULATE THE LOCAL ERROR FOR THE CURRENT STEP
2567 C     TO SEE IF THE STEP WAS SUCCESSFUL
2568 430   CONTINUE
2569       IERR = CK * ENORM
2570       IF(IERR .GT. 1.0D0)GO TO 600
2571 C
2572 C
2573 C
2574 C
2575 C
2576 C-----------------------------------------------------------------------
2577 C     BLOCK 5
2578 C     THE STEP IS SUCCESSFUL. DETERMINE
2579 C     THE BEST ORDER AND STEPSIZE FOR
2580 C     THE NEXT STEP. UPDATE THE DIFFERENCES
2581 C     FOR THE NEXT STEP.
2582 C-----------------------------------------------------------------------
2583       IDID=1
2584       IWM(LNST)=IWM(LNST)+1
2585       KDIFF=K-KOLD
2586       KOLD=K
2587       HOLD=H
2588 C
2589 C
2590 C     ESTIMATE THE ERROR AT ORDER K+1 UNLESS:
2591 C        ALREADY DECIDED TO LOWER ORDER, OR
2592 C        ALREADY USING MAXIMUM ORDER, OR
2593 C        STEPSIZE NOT CONSTANT, OR
2594 C        ORDER RAISED IN PREVIOUS STEP
2595       IF(KNEW.EQ.KM1.OR.K.EQ.IWM(LMXORD))IPHASE=1
2596       IF(IPHASE .EQ. 0)GO TO 545
2597       IF(KNEW.EQ.KM1)GO TO 540
2598       IF(K.EQ.IWM(LMXORD)) GO TO 550
2599       IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550
2600       DO 510 I=1,NEQ
2601 510      DELTA(I)=E(I)-PHI(I,KP2)
2602       ERKP1 = (1.0D0/(K+2))*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
2603       TERKP1 = (K+2)*ERKP1
2604       IF(K.GT.1)GO TO 520
2605       IF(TERKP1.GE.0.5D0*TERK)GO TO 550
2606       GO TO 530
2607 520   IF(TERKM1.LE.MIN(TERK,TERKP1))GO TO 540
2608       IF(TERKP1.GE.TERK.OR.K.EQ.IWM(LMXORD))GO TO 550
2609 C
2610 C     RAISE ORDER
2611 530   K=KP1
2612       EST = ERKP1
2613       GO TO 550
2614 C
2615 C     LOWER ORDER
2616 540   K=KM1
2617       EST = ERKM1
2618       GO TO 550
2619 C
2620 C     IF IPHASE = 0, INCREASE ORDER BY ONE AND MULTIPLY STEPSIZE BY
2621 C     FACTOR TWO
2622 545   K = KP1
2623       HNEW = H*2.0D0
2624       H = HNEW
2625       GO TO 575
2626 C
2627 C
2628 C     DETERMINE THE APPROPRIATE STEPSIZE FOR
2629 C     THE NEXT STEP.
2630 550   HNEW=H
2631       TEMP2=K+1
2632       R=(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
2633       IF(R .LT. 2.0D0) GO TO 555
2634       HNEW = 2.0D0*H
2635       GO TO 560
2636 555   IF(R .GT. 1.0D0) GO TO 560
2637       R = MAX(0.5D0,MIN(0.9D0,R))
2638       HNEW = H*R
2639 560   H=HNEW
2640 C
2641 C
2642 C     UPDATE DIFFERENCES FOR NEXT STEP
2643 575   CONTINUE
2644       IF(KOLD.EQ.IWM(LMXORD))GO TO 585
2645       DO 580 I=1,NEQ
2646 580      PHI(I,KP2)=E(I)
2647 585   CONTINUE
2648       DO 590 I=1,NEQ
2649 590      PHI(I,KP1)=PHI(I,KP1)+E(I)
2650       DO 595 J1=2,KP1
2651          J=KP1-J1+1
2652          DO 595 I=1,NEQ
2653 595      PHI(I,J)=PHI(I,J)+PHI(I,J+1)
2654       RETURN
2655 C
2656 C
2657 C
2658 C
2659 C
2660 C-----------------------------------------------------------------------
2661 C     BLOCK 6
2662 C     THE STEP IS UNSUCCESSFUL. RESTORE X,PSI,PHI
2663 C     DETERMINE APPROPRIATE STEPSIZE FOR
2664 C     CONTINUING THE INTEGRATION, OR EXIT WITH
2665 C     AN ERROR FLAG IF THERE HAVE BEEN MANY
2666 C     FAILURES.
2667 C-----------------------------------------------------------------------
2668 600   IPHASE = 1
2669 C
2670 C     RESTORE X,PHI,PSI
2671       X=XOLD
2672       IF(KP1.LT.NSP1)GO TO 630
2673       DO 620 J=NSP1,KP1
2674          TEMP1=1.0D0/BETA(J)
2675          DO 610 I=1,NEQ
2676 610         PHI(I,J)=TEMP1*PHI(I,J)
2677 620      CONTINUE
2678 630   CONTINUE
2679       DO 640 I=2,KP1
2680 640      PSI(I-1)=PSI(I)-H
2681 C
2682 C
2683 C     TEST WHETHER FAILURE IS DUE TO CORRECTOR ITERATION
2684 C     OR ERROR TEST
2685       IF(CONVGD)GO TO 660
2686       IWM(LCTF)=IWM(LCTF)+1
2687 C
2688 C
2689 C     THE NEWTON ITERATION FAILED TO CONVERGE WITH
2690 C     A CURRENT ITERATION MATRIX.  DETERMINE THE CAUSE
2691 C     OF THE FAILURE AND TAKE APPROPRIATE ACTION.
2692       IF(IER.EQ.0)GO TO 650
2693 C
2694 C     THE ITERATION MATRIX IS SINGULAR. REDUCE
2695 C     THE STEPSIZE BY A FACTOR OF 4. IF
2696 C     THIS HAPPENS THREE TIMES IN A ROW ON
2697 C     THE SAME STEP, RETURN WITH AN ERROR FLAG
2698       NSF=NSF+1
2699       R = 0.25D0
2700       H=H*R
2701       IF (NSF .LT. 3 .AND. ABS(H) .GE. HMIN) GO TO 690
2702       IDID=-8
2703       GO TO 675
2704 C
2705 C
2706 C     THE NEWTON ITERATION FAILED TO CONVERGE FOR A REASON
2707 C     OTHER THAN A SINGULAR ITERATION MATRIX.  IF IRES = -2, THEN
2708 C     RETURN.  OTHERWISE, REDUCE THE STEPSIZE AND TRY AGAIN, UNLESS
2709 C     TOO MANY FAILURES HAVE OCCURRED.
2710 650   CONTINUE
2711       IF (IRES .GT. -2) GO TO 655
2712       IDID = -11
2713       GO TO 675
2714 655   NCF = NCF + 1
2715       R = 0.25D0
2716       H = H*R
2717       IF (NCF .LT. 10 .AND. ABS(H) .GE. HMIN) GO TO 690
2718       IDID = -7
2719       IF (IRES .LT. 0) IDID = -10
2720       IF (NEF .GE. 3) IDID = -9
2721       GO TO 675
2722 C
2723 C
2724 C     THE NEWTON SCHEME CONVERGED,AND THE CAUSE
2725 C     OF THE FAILURE WAS THE ERROR ESTIMATE
2726 C     EXCEEDING THE TOLERANCE.
2727 660   NEF=NEF+1
2728       IWM(LETF)=IWM(LETF)+1
2729       IF (NEF .GT. 1) GO TO 665
2730 C
2731 C     ON FIRST ERROR TEST FAILURE, KEEP CURRENT ORDER OR LOWER
2732 C     ORDER BY ONE.  COMPUTE NEW STEPSIZE BASED ON DIFFERENCES
2733 C     OF THE SOLUTION.
2734       K = KNEW
2735       TEMP2 = K + 1
2736       R = 0.90D0*(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
2737       R = MAX(0.25D0,MIN(0.9D0,R))
2738       H = H*R
2739       IF (ABS(H) .GE. HMIN) GO TO 690
2740       IDID = -6
2741       GO TO 675
2742 C
2743 C     ON SECOND ERROR TEST FAILURE, USE THE CURRENT ORDER OR
2744 C     DECREASE ORDER BY ONE.  REDUCE THE STEPSIZE BY A FACTOR OF
2745 C     FOUR.
2746 665   IF (NEF .GT. 2) GO TO 670
2747       K = KNEW
2748       H = 0.25D0*H
2749       IF (ABS(H) .GE. HMIN) GO TO 690
2750       IDID = -6
2751       GO TO 675
2752 C
2753 C     ON THIRD AND SUBSEQUENT ERROR TEST FAILURES, SET THE ORDER TO
2754 C     ONE AND REDUCE THE STEPSIZE BY A FACTOR OF FOUR.
2755 670   K = 1
2756       H = 0.25D0*H
2757       IF (ABS(H) .GE. HMIN) GO TO 690
2758       IDID = -6
2759       GO TO 675
2760 C
2761 C
2762 C
2763 C
2764 C     FOR ALL CRASHES, RESTORE Y TO ITS LAST VALUE,
2765 C     INTERPOLATE TO FIND YPRIME AT LAST X, AND RETURN
2766 675   CONTINUE
2767       CALL DDATRP(X,X,Y,YPRIME,NEQ,K,PHI,PSI)
2768       RETURN
2769 C
2770 C
2771 C     GO BACK AND TRY THIS STEP AGAIN
2772 690   GO TO 200
2773 C
2774 C------END OF SUBROUTINE DDASTP------
2775       END
2776       SUBROUTINE DDATRP (X, XOUT, YOUT, YPOUT, NEQ, KOLD, PHI, PSI)
2777 C***BEGIN PROLOGUE  DDATRP
2778 C***SUBSIDIARY
2779 C***PURPOSE  Interpolation routine for DDASSL.
2780 C***LIBRARY   SLATEC (DASSL)
2781 C***TYPE      DOUBLE PRECISION (SDATRP-S, DDATRP-D)
2782 C***AUTHOR  PETZOLD, LINDA R., (LLNL)
2783 C***DESCRIPTION
2784 C-----------------------------------------------------------------------
2785 C     THE METHODS IN SUBROUTINE DDASTP USE POLYNOMIALS
2786 C     TO APPROXIMATE THE SOLUTION. DDATRP APPROXIMATES THE
2787 C     SOLUTION AND ITS DERIVATIVE AT TIME XOUT BY EVALUATING
2788 C     ONE OF THESE POLYNOMIALS,AND ITS DERIVATIVE,THERE.
2789 C     INFORMATION DEFINING THIS POLYNOMIAL IS PASSED FROM
2790 C     DDASTP, SO DDATRP CANNOT BE USED ALONE.
2791 C
2792 C     THE PARAMETERS ARE:
2793 C     X     THE CURRENT TIME IN THE INTEGRATION.
2794 C     XOUT  THE TIME AT WHICH THE SOLUTION IS DESIRED
2795 C     YOUT  THE INTERPOLATED APPROXIMATION TO Y AT XOUT
2796 C           (THIS IS OUTPUT)
2797 C     YPOUT THE INTERPOLATED APPROXIMATION TO YPRIME AT XOUT
2798 C           (THIS IS OUTPUT)
2799 C     NEQ   NUMBER OF EQUATIONS
2800 C     KOLD  ORDER USED ON LAST SUCCESSFUL STEP
2801 C     PHI   ARRAY OF SCALED DIVIDED DIFFERENCES OF Y
2802 C     PSI   ARRAY OF PAST STEPSIZE HISTORY
2803 C-----------------------------------------------------------------------
2804 C***ROUTINES CALLED  (NONE)
2805 C***REVISION HISTORY  (YYMMDD)
2806 C   830315  DATE WRITTEN
2807 C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
2808 C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
2809 C   901026  Added explicit declarations for all variables and minor
2810 C           cosmetic changes to prologue.  (FNF)
2811 C***END PROLOGUE  DDATRP
2812 C
2813       INTEGER  NEQ, KOLD
2814       DOUBLE PRECISION  X, XOUT, YOUT(*), YPOUT(*), PHI(NEQ,*), PSI(*)
2815 C
2816       INTEGER  I, J, KOLDP1
2817       DOUBLE PRECISION  C, D, GAMMA, TEMP1
2818 C
2819 C***FIRST EXECUTABLE STATEMENT  DDATRP
2820       KOLDP1=KOLD+1
2821       TEMP1=XOUT-X
2822       DO 10 I=1,NEQ
2823          YOUT(I)=PHI(I,1)
2824 10       YPOUT(I)=0.0D0
2825       C=1.0D0
2826       D=0.0D0
2827       GAMMA=TEMP1/PSI(1)
2828       DO 30 J=2,KOLDP1
2829          D=D*GAMMA+C/PSI(J-1)
2830          C=C*GAMMA
2831          GAMMA=(TEMP1+PSI(J-1))/PSI(J)
2832          DO 20 I=1,NEQ
2833             YOUT(I)=YOUT(I)+C*PHI(I,J)
2834 20          YPOUT(I)=YPOUT(I)+D*PHI(I,J)
2835 30       CONTINUE
2836       RETURN
2837 C
2838 C------END OF SUBROUTINE DDATRP------
2839       END
2840       SUBROUTINE DDAWTS (NEQ, IWT, RTOL, ATOL, Y, WT, RPAR, IPAR)
2841 C***BEGIN PROLOGUE  DDAWTS
2842 C***SUBSIDIARY
2843 C***PURPOSE  Set error weight vector for DDASSL.
2844 C***LIBRARY   SLATEC (DASSL)
2845 C***TYPE      DOUBLE PRECISION (SDAWTS-S, DDAWTS-D)
2846 C***AUTHOR  PETZOLD, LINDA R., (LLNL)
2847 C***DESCRIPTION
2848 C-----------------------------------------------------------------------
2849 C     THIS SUBROUTINE SETS THE ERROR WEIGHT VECTOR
2850 C     WT ACCORDING TO WT(I)=RTOL(I)*ABS(Y(I))+ATOL(I),
2851 C     I=1,-,N.
2852 C     RTOL AND ATOL ARE SCALARS IF IWT = 0,
2853 C     AND VECTORS IF IWT = 1.
2854 C-----------------------------------------------------------------------
2855 C***ROUTINES CALLED  (NONE)
2856 C***REVISION HISTORY  (YYMMDD)
2857 C   830315  DATE WRITTEN
2858 C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
2859 C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
2860 C   901026  Added explicit declarations for all variables and minor
2861 C           cosmetic changes to prologue.  (FNF)
2862 C***END PROLOGUE  DDAWTS
2863 C
2864       INTEGER  NEQ, IWT, IPAR(*)
2865       DOUBLE PRECISION  RTOL(*), ATOL(*), Y(*), WT(*), RPAR(*)
2866 C
2867       INTEGER  I
2868       DOUBLE PRECISION  ATOLI, RTOLI
2869 C
2870 C***FIRST EXECUTABLE STATEMENT  DDAWTS
2871       RTOLI=RTOL(1)
2872       ATOLI=ATOL(1)
2873       DO 20 I=1,NEQ
2874          IF (IWT .EQ.0) GO TO 10
2875            RTOLI=RTOL(I)
2876            ATOLI=ATOL(I)
2877 10         WT(I)=RTOLI*ABS(Y(I))+ATOLI
2878 20         CONTINUE
2879       RETURN
2880 C-----------END OF SUBROUTINE DDAWTS------------------------------------
2881       END