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