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