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