Mercurial > octave
changeset 1644:395bb6d6c096
[project @ 1995-12-14 08:32:12 by jwe]
Initial revision
| author | jwe |
|---|---|
| date | Thu, 14 Dec 1995 08:32:21 +0000 |
| parents | 5e108d51e370 |
| children | 44ed237bdc1e |
| files | libcruft/fftpack/cfftb1.f libcruft/fftpack/cfftf1.f libcruft/fftpack/cffti1.f libcruft/odepack/prepj.f libcruft/odepack/solsy.f libcruft/odepack/stode.f |
| diffstat | 6 files changed, 913 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/fftpack/cfftb1.f Thu Dec 14 08:32:21 1995 +0000 @@ -0,0 +1,62 @@ + subroutine cfftb1 (n,c,ch,wa,ifac) + implicit double precision (a-h,o-z) + dimension ch(1) ,c(1) ,wa(1) ,ifac(1) + nf = ifac(2) + na = 0 + l1 = 1 + iw = 1 + do 116 k1=1,nf + ip = ifac(k1+2) + l2 = ip*l1 + ido = n/l2 + idot = ido+ido + idl1 = idot*l1 + if (ip .ne. 4) go to 103 + ix2 = iw+idot + ix3 = ix2+idot + if (na .ne. 0) go to 101 + call passb4 (idot,l1,c,ch,wa(iw),wa(ix2),wa(ix3)) + go to 102 + 101 call passb4 (idot,l1,ch,c,wa(iw),wa(ix2),wa(ix3)) + 102 na = 1-na + go to 115 + 103 if (ip .ne. 2) go to 106 + if (na .ne. 0) go to 104 + call passb2 (idot,l1,c,ch,wa(iw)) + go to 105 + 104 call passb2 (idot,l1,ch,c,wa(iw)) + 105 na = 1-na + go to 115 + 106 if (ip .ne. 3) go to 109 + ix2 = iw+idot + if (na .ne. 0) go to 107 + call passb3 (idot,l1,c,ch,wa(iw),wa(ix2)) + go to 108 + 107 call passb3 (idot,l1,ch,c,wa(iw),wa(ix2)) + 108 na = 1-na + go to 115 + 109 if (ip .ne. 5) go to 112 + ix2 = iw+idot + ix3 = ix2+idot + ix4 = ix3+idot + if (na .ne. 0) go to 110 + call passb5 (idot,l1,c,ch,wa(iw),wa(ix2),wa(ix3),wa(ix4)) + go to 111 + 110 call passb5 (idot,l1,ch,c,wa(iw),wa(ix2),wa(ix3),wa(ix4)) + 111 na = 1-na + go to 115 + 112 if (na .ne. 0) go to 113 + call passb (nac,idot,ip,l1,idl1,c,c,c,ch,ch,wa(iw)) + go to 114 + 113 call passb (nac,idot,ip,l1,idl1,ch,ch,ch,c,c,wa(iw)) + 114 if (nac .ne. 0) na = 1-na + 115 l1 = l2 + iw = iw+(ip-1)*idot + 116 continue + if (na .eq. 0) return + n2 = n+n + do 117 i=1,n2 + c(i) = ch(i) + 117 continue + return + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/fftpack/cfftf1.f Thu Dec 14 08:32:21 1995 +0000 @@ -0,0 +1,62 @@ + subroutine cfftf1 (n,c,ch,wa,ifac) + implicit double precision (a-h,o-z) + dimension ch(1) ,c(1) ,wa(1) ,ifac(1) + nf = ifac(2) + na = 0 + l1 = 1 + iw = 1 + do 116 k1=1,nf + ip = ifac(k1+2) + l2 = ip*l1 + ido = n/l2 + idot = ido+ido + idl1 = idot*l1 + if (ip .ne. 4) go to 103 + ix2 = iw+idot + ix3 = ix2+idot + if (na .ne. 0) go to 101 + call passf4 (idot,l1,c,ch,wa(iw),wa(ix2),wa(ix3)) + go to 102 + 101 call passf4 (idot,l1,ch,c,wa(iw),wa(ix2),wa(ix3)) + 102 na = 1-na + go to 115 + 103 if (ip .ne. 2) go to 106 + if (na .ne. 0) go to 104 + call passf2 (idot,l1,c,ch,wa(iw)) + go to 105 + 104 call passf2 (idot,l1,ch,c,wa(iw)) + 105 na = 1-na + go to 115 + 106 if (ip .ne. 3) go to 109 + ix2 = iw+idot + if (na .ne. 0) go to 107 + call passf3 (idot,l1,c,ch,wa(iw),wa(ix2)) + go to 108 + 107 call passf3 (idot,l1,ch,c,wa(iw),wa(ix2)) + 108 na = 1-na + go to 115 + 109 if (ip .ne. 5) go to 112 + ix2 = iw+idot + ix3 = ix2+idot + ix4 = ix3+idot + if (na .ne. 0) go to 110 + call passf5 (idot,l1,c,ch,wa(iw),wa(ix2),wa(ix3),wa(ix4)) + go to 111 + 110 call passf5 (idot,l1,ch,c,wa(iw),wa(ix2),wa(ix3),wa(ix4)) + 111 na = 1-na + go to 115 + 112 if (na .ne. 0) go to 113 + call passf (nac,idot,ip,l1,idl1,c,c,c,ch,ch,wa(iw)) + go to 114 + 113 call passf (nac,idot,ip,l1,idl1,ch,ch,ch,c,c,wa(iw)) + 114 if (nac .ne. 0) na = 1-na + 115 l1 = l2 + iw = iw+(ip-1)*idot + 116 continue + if (na .eq. 0) return + n2 = n+n + do 117 i=1,n2 + c(i) = ch(i) + 117 continue + return + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/fftpack/cffti1.f Thu Dec 14 08:32:21 1995 +0000 @@ -0,0 +1,61 @@ + subroutine cffti1 (n,wa,ifac) + implicit double precision (a-h,o-z) + dimension wa(1) ,ifac(1) ,ntryh(4) + data ntryh(1),ntryh(2),ntryh(3),ntryh(4)/3,4,2,5/ + nl = n + nf = 0 + j = 0 + 101 j = j+1 + if (j-4) 102,102,103 + 102 ntry = ntryh(j) + go to 104 + 103 ntry = ntry+2 + 104 nq = nl/ntry + nr = nl-ntry*nq + if (nr) 101,105,101 + 105 nf = nf+1 + ifac(nf+2) = ntry + nl = nq + if (ntry .ne. 2) go to 107 + if (nf .eq. 1) go to 107 + do 106 i=2,nf + ib = nf-i+2 + ifac(ib+2) = ifac(ib+1) + 106 continue + ifac(3) = 2 + 107 if (nl .ne. 1) go to 104 + ifac(1) = n + ifac(2) = nf + tpi = 6.28318530717959d0 + argh = tpi/dble(n) + i = 2 + l1 = 1 + do 110 k1=1,nf + ip = ifac(k1+2) + ld = 0 + l2 = l1*ip + ido = n/l2 + idot = ido+ido+2 + ipm = ip-1 + do 109 j=1,ipm + i1 = i + wa(i-1) = 1. + wa(i) = 0. + ld = ld+l1 + fi = 0. + argld = dble(ld)*argh + do 108 ii=4,idot,2 + i = i+2 + fi = fi+1. + arg = fi*argld + wa(i-1) = cos(arg) + wa(i) = sin(arg) + 108 continue + if (ip .le. 5) go to 109 + wa(i1-1) = wa(i-1) + wa(i1) = wa(i) + 109 continue + l1 = l2 + 110 continue + return + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odepack/prepj.f Thu Dec 14 08:32:21 1995 +0000 @@ -0,0 +1,177 @@ + SUBROUTINE PREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, + 1 F, JAC, IERR) +CLLL. OPTIMIZE + EXTERNAL F, JAC + INTEGER NEQ, NYH, IWM + INTEGER IOWND, IOWNS, + 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, + 2 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP, + 1 MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1 + DOUBLE PRECISION Y, YH, EWT, FTEM, SAVF, WM + DOUBLE PRECISION ROWNS, + 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND + DOUBLE PRECISION CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ, + 1 VNORM + DIMENSION NEQ(1), Y(1), YH(NYH,1), EWT(1), FTEM(1), SAVF(1), + 1 WM(1), IWM(1) + COMMON /LS0001/ ROWNS(209), + 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, + 3 IOWND(14), IOWNS(6), + 4 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, + 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU +C----------------------------------------------------------------------- +C PREPJ IS CALLED BY STODE TO COMPUTE AND PROCESS THE MATRIX +C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN. +C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF +C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5. +C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED. +C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN +C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION +C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE +C BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5. +C +C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION +C WITH PREPJ USES THE FOLLOWING.. +C Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY. +C FTEM = WORK ARRAY OF LENGTH N (ACOR IN STODE). +C SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y. +C WM = REAL WORK SPACE FOR MATRICES. ON OUTPUT IT CONTAINS THE +C INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION +C OF P IF MITER IS 1, 2 , 4, OR 5. +C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3). +C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA.. +C WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS. +C WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3. +C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT +C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND +C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5. +C EL0 = EL(1) (INPUT). +C IERPJ = OUTPUT ERROR FLAG, = 0 IF NO TROUBLE, .GT. 0 IF +C P MATRIX FOUND TO BE SINGULAR. +C JCUR = OUTPUT FLAG = 1 TO INDICATE THAT THE JACOBIAN MATRIX +C (OR APPROXIMATION) IS NOW CURRENT. +C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND, +C MITER, N, NFE, AND NJE. +C----------------------------------------------------------------------- + NJE = NJE + 1 + IERPJ = 0 + JCUR = 1 + HL0 = H*EL0 + GO TO (100, 200, 300, 400, 500), MITER +C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. ----------------------- + 100 LENP = N*N + DO 110 I = 1,LENP + 110 WM(I+2) = 0.0D0 + CALL JAC (NEQ, TN, Y, 0, 0, WM(3), N) + CON = -HL0 + DO 120 I = 1,LENP + 120 WM(I+2) = WM(I+2)*CON + GO TO 240 +C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. -------------------- + 200 FAC = VNORM (N, SAVF, EWT) + R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC + IF (R0 .EQ. 0.0D0) R0 = 1.0D0 + SRUR = WM(1) + J1 = 2 + DO 230 J = 1,N + YJ = Y(J) + R = DMAX1(SRUR*DABS(YJ),R0/EWT(J)) + Y(J) = Y(J) + R + FAC = -HL0/R + IERR = 0 + CALL F (NEQ, TN, Y, FTEM, IERR) + IF (IERR .LT. 0) RETURN + DO 220 I = 1,N + 220 WM(I+J1) = (FTEM(I) - SAVF(I))*FAC + Y(J) = YJ + J1 = J1 + N + 230 CONTINUE + NFE = NFE + N +C ADD IDENTITY MATRIX. ------------------------------------------------- + 240 J = 3 + NP1 = N + 1 + DO 250 I = 1,N + WM(J) = WM(J) + 1.0D0 + 250 J = J + NP1 +C DO LU DECOMPOSITION ON P. -------------------------------------------- + CALL DGEFA (WM(3), N, N, IWM(21), IER) + IF (IER .NE. 0) IERPJ = 1 + RETURN +C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. --------- + 300 WM(2) = HL0 + R = EL0*0.1D0 + DO 310 I = 1,N + 310 Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2)) + IERR = 0 + CALL F (NEQ, TN, Y, WM(3), IERR) + IF (IERR .LT. 0) RETURN + NFE = NFE + 1 + DO 320 I = 1,N + R0 = H*SAVF(I) - YH(I,2) + DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I)) + WM(I+2) = 1.0D0 + IF (DABS(R0) .LT. UROUND/EWT(I)) GO TO 320 + IF (DABS(DI) .EQ. 0.0D0) GO TO 330 + WM(I+2) = 0.1D0*R0/DI + 320 CONTINUE + RETURN + 330 IERPJ = 1 + RETURN +C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. ----------------------- + 400 ML = IWM(1) + MU = IWM(2) + ML3 = ML + 3 + MBAND = ML + MU + 1 + MEBAND = MBAND + ML + LENP = MEBAND*N + DO 410 I = 1,LENP + 410 WM(I+2) = 0.0D0 + CALL JAC (NEQ, TN, Y, ML, MU, WM(ML3), MEBAND) + CON = -HL0 + DO 420 I = 1,LENP + 420 WM(I+2) = WM(I+2)*CON + GO TO 570 +C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ---------------- + 500 ML = IWM(1) + MU = IWM(2) + MBAND = ML + MU + 1 + MBA = MIN0(MBAND,N) + MEBAND = MBAND + ML + MEB1 = MEBAND - 1 + SRUR = WM(1) + FAC = VNORM (N, SAVF, EWT) + R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC + IF (R0 .EQ. 0.0D0) R0 = 1.0D0 + DO 560 J = 1,MBA + DO 530 I = J,N,MBAND + YI = Y(I) + R = DMAX1(SRUR*DABS(YI),R0/EWT(I)) + 530 Y(I) = Y(I) + R + IERR = 0 + CALL F (NEQ, TN, Y, FTEM, IERR) + IF (IERR .LT. 0) RETURN + DO 550 JJ = J,N,MBAND + Y(JJ) = YH(JJ,1) + YJJ = Y(JJ) + R = DMAX1(SRUR*DABS(YJJ),R0/EWT(JJ)) + FAC = -HL0/R + I1 = MAX0(JJ-MU,1) + I2 = MIN0(JJ+ML,N) + II = JJ*MEB1 - ML + 2 + DO 540 I = I1,I2 + 540 WM(II+I) = (FTEM(I) - SAVF(I))*FAC + 550 CONTINUE + 560 CONTINUE + NFE = NFE + MBA +C ADD IDENTITY MATRIX. ------------------------------------------------- + 570 II = MBAND + 2 + DO 580 I = 1,N + WM(II) = WM(II) + 1.0D0 + 580 II = II + MEBAND +C DO LU DECOMPOSITION OF P. -------------------------------------------- + CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER) + IF (IER .NE. 0) IERPJ = 1 + RETURN +C----------------------- END OF SUBROUTINE PREPJ ----------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odepack/solsy.f Thu Dec 14 08:32:21 1995 +0000 @@ -0,0 +1,68 @@ + SUBROUTINE SOLSY (WM, IWM, X, TEM) +CLLL. OPTIMIZE + INTEGER IWM + INTEGER IOWND, IOWNS, + 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, + 2 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + INTEGER I, MEBAND, ML, MU + DOUBLE PRECISION WM, X, TEM + DOUBLE PRECISION ROWNS, + 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND + DOUBLE PRECISION DI, HL0, PHL0, R + DIMENSION WM(1), IWM(1), X(1), TEM(1) + COMMON /LS0001/ ROWNS(209), + 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, + 3 IOWND(14), IOWNS(6), + 4 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, + 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU +C----------------------------------------------------------------------- +C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING FROM +C A CHORD ITERATION. IT IS CALLED IF MITER .NE. 0. +C IF MITER IS 1 OR 2, IT CALLS DGESL TO ACCOMPLISH THIS. +C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL +C MATRIX, AND THEN COMPUTES THE SOLUTION. +C IF MITER IS 4 OR 5, IT CALLS DGBSL. +C COMMUNICATION WITH SOLSY USES THE FOLLOWING VARIABLES.. +C WM = REAL WORK SPACE CONTAINING THE INVERSE DIAGONAL MATRIX IF +C MITER = 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE. +C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3). +C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA.. +C WM(1) = SQRT(UROUND) (NOT USED HERE), +C WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER = 3. +C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT +C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND +C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5. +C X = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION VECTOR +C ON OUTPUT, OF LENGTH N. +C TEM = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION. +C IERSL = OUTPUT FLAG (IN COMMON). IERSL = 0 IF NO TROUBLE OCCURRED. +C IERSL = 1 IF A SINGULAR MATRIX AROSE WITH MITER = 3. +C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N. +C----------------------------------------------------------------------- + IERSL = 0 + GO TO (100, 100, 300, 400, 400), MITER + 100 CALL DGESL (WM(3), N, N, IWM(21), X, 0) + RETURN +C + 300 PHL0 = WM(2) + HL0 = H*EL0 + WM(2) = HL0 + IF (HL0 .EQ. PHL0) GO TO 330 + R = HL0/PHL0 + DO 320 I = 1,N + DI = 1.0D0 - R*(1.0D0 - 1.0D0/WM(I+2)) + IF (DABS(DI) .EQ. 0.0D0) GO TO 390 + 320 WM(I+2) = 1.0D0/DI + 330 DO 340 I = 1,N + 340 X(I) = WM(I+2)*X(I) + RETURN + 390 IERSL = 1 + RETURN +C + 400 ML = IWM(1) + MU = IWM(2) + MEBAND = 2*ML + MU + 1 + CALL DGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0) + RETURN +C----------------------- END OF SUBROUTINE SOLSY ----------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odepack/stode.f Thu Dec 14 08:32:21 1995 +0000 @@ -0,0 +1,483 @@ + SUBROUTINE STODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, + 1 WM, IWM, F, JAC, PJAC, SLVS, IERR) +CLLL. OPTIMIZE + EXTERNAL F, JAC, PJAC, SLVS + INTEGER NEQ, NYH, IWM + INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, + 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, + 2 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ + DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, ACOR, WM + DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO, + 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND + DOUBLE PRECISION DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP, + 1 R, RH, RHDN, RHSM, RHUP, TOLD, VNORM + DIMENSION NEQ(1), Y(1), YH(NYH,1), YH1(1), EWT(1), SAVF(1), + 1 ACOR(1), WM(1), IWM(1) + COMMON /LS0001/ CONIT, CRATE, EL(13), ELCO(13,12), + 1 HOLD, RMAX, TESCO(3,12), + 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, IOWND(14), + 3 IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, + 4 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, + 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU +C----------------------------------------------------------------------- +C STODE PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE +C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS. +C NOTE.. STODE IS INDEPENDENT OF THE VALUE OF THE ITERATION METHOD +C INDICATOR MITER, WHEN THIS IS .NE. 0, AND HENCE IS INDEPENDENT +C OF THE TYPE OF CHORD METHOD USED, OR THE JACOBIAN STRUCTURE. +C COMMUNICATION WITH STODE IS DONE WITH THE FOLLOWING VARIABLES.. +C +C NEQ = INTEGER ARRAY CONTAINING PROBLEM SIZE IN NEQ(1), AND +C PASSED AS THE NEQ ARGUMENT IN ALL CALLS TO F AND JAC. +C Y = AN ARRAY OF LENGTH .GE. N USED AS THE Y ARGUMENT IN +C ALL CALLS TO F AND JAC. +C YH = AN NYH BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES +C AND THEIR APPROXIMATE SCALED DERIVATIVES, WHERE +C LMAX = MAXORD + 1. YH(I,J+1) CONTAINS THE APPROXIMATE +C J-TH DERIVATIVE OF Y(I), SCALED BY H**J/FACTORIAL(J) +C (J = 0,1,...,NQ). ON ENTRY FOR THE FIRST STEP, THE FIRST +C TWO COLUMNS OF YH MUST BE SET FROM THE INITIAL VALUES. +C NYH = A CONSTANT INTEGER .GE. N, THE FIRST DIMENSION OF YH. +C YH1 = A ONE-DIMENSIONAL ARRAY OCCUPYING THE SAME SPACE AS YH. +C EWT = AN ARRAY OF LENGTH N CONTAINING MULTIPLICATIVE WEIGHTS +C FOR LOCAL ERROR MEASUREMENTS. LOCAL ERRORS IN Y(I) ARE +C COMPARED TO 1.0/EWT(I) IN VARIOUS ERROR TESTS. +C SAVF = AN ARRAY OF WORKING STORAGE, OF LENGTH N. +C ALSO USED FOR INPUT OF YH(*,MAXORD+2) WHEN JSTART = -1 +C AND MAXORD .LT. THE CURRENT ORDER NQ. +C ACOR = A WORK ARRAY OF LENGTH N, USED FOR THE ACCUMULATED +C CORRECTIONS. ON A SUCCESSFUL RETURN, ACOR(I) CONTAINS +C THE ESTIMATED ONE-STEP LOCAL ERROR IN Y(I). +C WM,IWM = REAL AND INTEGER WORK ARRAYS ASSOCIATED WITH MATRIX +C OPERATIONS IN CHORD ITERATION (MITER .NE. 0). +C PJAC = NAME OF ROUTINE TO EVALUATE AND PREPROCESS JACOBIAN MATRIX +C AND P = I - H*EL0*JAC, IF A CHORD METHOD IS BEING USED. +C SLVS = NAME OF ROUTINE TO SOLVE LINEAR SYSTEM IN CHORD ITERATION. +C CCMAX = MAXIMUM RELATIVE CHANGE IN H*EL0 BEFORE PJAC IS CALLED. +C H = THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. +C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE +C PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS +C SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM. +C HMIN = THE MINIMUM ABSOLUTE VALUE OF THE STEP SIZE H TO BE USED. +C HMXI = INVERSE OF THE MAXIMUM ABSOLUTE VALUE OF H TO BE USED. +C HMXI = 0.0 IS ALLOWED AND CORRESPONDS TO AN INFINITE HMAX. +C HMIN AND HMXI MAY BE CHANGED AT ANY TIME, BUT WILL NOT +C TAKE EFFECT UNTIL THE NEXT CHANGE OF H IS CONSIDERED. +C TN = THE INDEPENDENT VARIABLE. TN IS UPDATED ON EACH STEP TAKEN. +C JSTART = AN INTEGER USED FOR INPUT ONLY, WITH THE FOLLOWING +C VALUES AND MEANINGS.. +C 0 PERFORM THE FIRST STEP. +C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST. +C -1 TAKE THE NEXT STEP WITH A NEW VALUE OF H, MAXORD, +C N, METH, MITER, AND/OR MATRIX PARAMETERS. +C -2 TAKE THE NEXT STEP WITH A NEW VALUE OF H, +C BUT WITH OTHER INPUTS UNCHANGED. +C ON RETURN, JSTART IS SET TO 1 TO FACILITATE CONTINUATION. +C KFLAG = A COMPLETION CODE WITH THE FOLLOWING MEANINGS.. +C 0 THE STEP WAS SUCCESFUL. +C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED. +C -2 CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED. +C -3 FATAL ERROR IN PJAC OR SLVS. +C A RETURN WITH KFLAG = -1 OR -2 MEANS EITHER +C ABS(H) = HMIN OR 10 CONSECUTIVE FAILURES OCCURRED. +C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF TN AND +C THE YH ARRAY ARE AS OF THE BEGINNING OF THE LAST +C STEP, AND H IS THE LAST STEP SIZE ATTEMPTED. +C MAXORD = THE MAXIMUM ORDER OF INTEGRATION METHOD TO BE ALLOWED. +C MAXCOR = THE MAXIMUM NUMBER OF CORRECTOR ITERATIONS ALLOWED. +C MSBP = MAXIMUM NUMBER OF STEPS BETWEEN PJAC CALLS (MITER .GT. 0). +C MXNCF = MAXIMUM NUMBER OF CONVERGENCE FAILURES ALLOWED. +C METH/MITER = THE METHOD FLAGS. SEE DESCRIPTION IN DRIVER. +C N = THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS. +C IERR = ERROR FLAG FROM USER-SUPPLIED FUNCTION +C----------------------------------------------------------------------- + KFLAG = 0 + TOLD = TN + NCF = 0 + IERPJ = 0 + IERSL = 0 + JCUR = 0 + ICF = 0 + DELP = 0.0D0 + IF (JSTART .GT. 0) GO TO 200 + IF (JSTART .EQ. -1) GO TO 100 + IF (JSTART .EQ. -2) GO TO 160 +C----------------------------------------------------------------------- +C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE +C INITIALIZED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED +C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL +C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE +C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2 +C FOR THE NEXT INCREASE. +C----------------------------------------------------------------------- + LMAX = MAXORD + 1 + NQ = 1 + L = 2 + IALTH = 2 + RMAX = 10000.0D0 + RC = 0.0D0 + EL0 = 1.0D0 + CRATE = 0.7D0 + HOLD = H + MEO = METH + NSLP = 0 + IPUP = MITER + IRET = 3 + GO TO 140 +C----------------------------------------------------------------------- +C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1. +C IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE. +C IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1), +C IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP. +C IF THE CALLER HAS CHANGED METH, CFODE IS CALLED TO RESET +C THE COEFFICIENTS OF THE METHOD. +C IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT +C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY. +C IF H IS TO BE CHANGED, YH MUST BE RESCALED. +C IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1 +C TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS. +C----------------------------------------------------------------------- + 100 IPUP = MITER + LMAX = MAXORD + 1 + IF (IALTH .EQ. 1) IALTH = 2 + IF (METH .EQ. MEO) GO TO 110 + CALL CFODE (METH, ELCO, TESCO) + MEO = METH + IF (NQ .GT. MAXORD) GO TO 120 + IALTH = L + IRET = 1 + GO TO 150 + 110 IF (NQ .LE. MAXORD) GO TO 160 + 120 NQ = MAXORD + L = LMAX + DO 125 I = 1,L + 125 EL(I) = ELCO(I,NQ) + NQNYH = NQ*NYH + RC = RC*EL(1)/EL0 + EL0 = EL(1) + CONIT = 0.5D0/DBLE(NQ+2) + DDN = VNORM (N, SAVF, EWT)/TESCO(1,L) + EXDN = 1.0D0/DBLE(L) + RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0) + RH = DMIN1(RHDN,1.0D0) + IREDO = 3 + IF (H .EQ. HOLD) GO TO 170 + RH = DMIN1(RH,DABS(H/HOLD)) + H = HOLD + GO TO 175 +C----------------------------------------------------------------------- +C CFODE IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE +C CURRENT METH. THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET +C WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM. +C----------------------------------------------------------------------- + 140 CALL CFODE (METH, ELCO, TESCO) + 150 DO 155 I = 1,L + 155 EL(I) = ELCO(I,NQ) + NQNYH = NQ*NYH + RC = RC*EL(1)/EL0 + EL0 = EL(1) + CONIT = 0.5D0/DBLE(NQ+2) + GO TO (160, 170, 200), IRET +C----------------------------------------------------------------------- +C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST +C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH IS SET TO +C L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS +C FORCED BY A CONVERGENCE OR ERROR TEST FAILURE. +C----------------------------------------------------------------------- + 160 IF (H .EQ. HOLD) GO TO 200 + RH = H/HOLD + H = HOLD + IREDO = 3 + GO TO 175 + 170 RH = DMAX1(RH,HMIN/DABS(H)) + 175 RH = DMIN1(RH,RMAX) + RH = RH/DMAX1(1.0D0,DABS(H)*HMXI*RH) + R = 1.0D0 + DO 180 J = 2,L + R = R*RH + DO 180 I = 1,N + 180 YH(I,J) = YH(I,J)*R + H = H*RH + RC = RC*RH + IALTH = L + IF (IREDO .EQ. 0) GO TO 690 +C----------------------------------------------------------------------- +C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY +C MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX. +C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1). +C WHEN RC DIFFERS FROM 1 BY MORE THAN CCMAX, IPUP IS SET TO MITER +C TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED. +C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY MSBP STEPS. +C----------------------------------------------------------------------- + 200 IF (DABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER + IF (NST .GE. NSLP+MSBP) IPUP = MITER + TN = TN + H + I1 = NQNYH + 1 + DO 215 JB = 1,NQ + I1 = I1 - NYH +CDIR$ IVDEP + DO 210 I = I1,NQNYH + 210 YH1(I) = YH1(I) + YH1(I+NYH) + 215 CONTINUE +C----------------------------------------------------------------------- +C UP TO MAXCOR CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS +C MADE ON THE R.M.S. NORM OF EACH CORRECTION, WEIGHTED BY THE ERROR +C WEIGHT VECTOR EWT. THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE +C VECTOR ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP. +C----------------------------------------------------------------------- + 220 M = 0 + DO 230 I = 1,N + 230 Y(I) = YH(I,1) + IERR = 0 + CALL F (NEQ, TN, Y, SAVF, IERR) + IF (IERR .LT. 0) RETURN + NFE = NFE + 1 + IF (IPUP .LE. 0) GO TO 250 +C----------------------------------------------------------------------- +C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED AND +C PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION. IPUP IS SET +C TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE. +C----------------------------------------------------------------------- + IERR = 0 + CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC, + 1 IERR) + IF (IERR .LT. 0) RETURN + IPUP = 0 + RC = 1.0D0 + NSLP = NST + CRATE = 0.7D0 + IF (IERPJ .NE. 0) GO TO 430 + 250 DO 260 I = 1,N + 260 ACOR(I) = 0.0D0 + 270 IF (MITER .NE. 0) GO TO 350 +C----------------------------------------------------------------------- +C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM +C THE RESULT OF THE LAST FUNCTION EVALUATION. +C----------------------------------------------------------------------- + DO 290 I = 1,N + SAVF(I) = H*SAVF(I) - YH(I,2) + 290 Y(I) = SAVF(I) - ACOR(I) + DEL = VNORM (N, Y, EWT) + DO 300 I = 1,N + Y(I) = YH(I,1) + EL(1)*SAVF(I) + 300 ACOR(I) = SAVF(I) + GO TO 400 +C----------------------------------------------------------------------- +C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR, +C AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND +C P AS COEFFICIENT MATRIX. +C----------------------------------------------------------------------- + 350 DO 360 I = 1,N + 360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I)) + CALL SLVS (WM, IWM, Y, SAVF) + IF (IERSL .LT. 0) GO TO 430 + IF (IERSL .GT. 0) GO TO 410 + DEL = VNORM (N, Y, EWT) + DO 380 I = 1,N + ACOR(I) = ACOR(I) + Y(I) + 380 Y(I) = YH(I,1) + EL(1)*ACOR(I) +C----------------------------------------------------------------------- +C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE +C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. +C----------------------------------------------------------------------- + 400 IF (M .NE. 0) CRATE = DMAX1(0.2D0*CRATE,DEL/DELP) + DCON = DEL*DMIN1(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT) + IF (DCON .LE. 1.0D0) GO TO 450 + M = M + 1 + IF (M .EQ. MAXCOR) GO TO 410 + IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410 + DELP = DEL + IERR = 0 + CALL F (NEQ, TN, Y, SAVF, IERR) + IF (IERR .LT. 0) RETURN + NFE = NFE + 1 + GO TO 270 +C----------------------------------------------------------------------- +C THE CORRECTOR ITERATION FAILED TO CONVERGE. +C IF MITER .NE. 0 AND THE JACOBIAN IS OUT OF DATE, PJAC IS CALLED FOR +C THE NEXT TRY. OTHERWISE THE YH ARRAY IS RETRACTED TO ITS VALUES +C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF H CANNOT BE +C REDUCED OR MXNCF FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2. +C----------------------------------------------------------------------- + 410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430 + ICF = 1 + IPUP = MITER + GO TO 220 + 430 ICF = 2 + NCF = NCF + 1 + RMAX = 2.0D0 + TN = TOLD + I1 = NQNYH + 1 + DO 445 JB = 1,NQ + I1 = I1 - NYH +CDIR$ IVDEP + DO 440 I = I1,NQNYH + 440 YH1(I) = YH1(I) - YH1(I+NYH) + 445 CONTINUE + IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680 + IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 670 + IF (NCF .EQ. MXNCF) GO TO 670 + RH = 0.25D0 + IPUP = MITER + IREDO = 1 + GO TO 170 +C----------------------------------------------------------------------- +C THE CORRECTOR HAS CONVERGED. JCUR IS SET TO 0 +C TO SIGNAL THAT THE JACOBIAN INVOLVED MAY NEED UPDATING LATER. +C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500 +C IF IT FAILS. +C----------------------------------------------------------------------- + 450 JCUR = 0 + IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ) + IF (M .GT. 0) DSM = VNORM (N, ACOR, EWT)/TESCO(2,NQ) + IF (DSM .GT. 1.0D0) GO TO 500 +C----------------------------------------------------------------------- +C AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY. +C CONSIDER CHANGING H IF IALTH = 1. OTHERWISE DECREASE IALTH BY 1. +C IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR +C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP. +C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER +C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A +C FACTOR OF AT LEAST 1.1. IF NOT, IALTH IS SET TO 3 TO PREVENT +C TESTING FOR THAT MANY STEPS. +C----------------------------------------------------------------------- + KFLAG = 0 + IREDO = 0 + NST = NST + 1 + HU = H + NQU = NQ + DO 470 J = 1,L + DO 470 I = 1,N + 470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I) + IALTH = IALTH - 1 + IF (IALTH .EQ. 0) GO TO 520 + IF (IALTH .GT. 1) GO TO 700 + IF (L .EQ. LMAX) GO TO 700 + DO 490 I = 1,N + 490 YH(I,LMAX) = ACOR(I) + GO TO 700 +C----------------------------------------------------------------------- +C THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES. +C RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE +C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR +C ONE LOWER ORDER. AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE +C BY A FACTOR OF 0.2 OR LESS. +C----------------------------------------------------------------------- + 500 KFLAG = KFLAG - 1 + TN = TOLD + I1 = NQNYH + 1 + DO 515 JB = 1,NQ + I1 = I1 - NYH +CDIR$ IVDEP + DO 510 I = I1,NQNYH + 510 YH1(I) = YH1(I) - YH1(I+NYH) + 515 CONTINUE + RMAX = 2.0D0 + IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 660 + IF (KFLAG .LE. -3) GO TO 640 + IREDO = 2 + RHUP = 0.0D0 + GO TO 540 +C----------------------------------------------------------------------- +C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS +C RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED +C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY. +C IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE. +C THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN +C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE +C ADDITIONAL SCALED DERIVATIVE. +C----------------------------------------------------------------------- + 520 RHUP = 0.0D0 + IF (L .EQ. LMAX) GO TO 540 + DO 530 I = 1,N + 530 SAVF(I) = ACOR(I) - YH(I,LMAX) + DUP = VNORM (N, SAVF, EWT)/TESCO(3,NQ) + EXUP = 1.0D0/DBLE(L+1) + RHUP = 1.0D0/(1.4D0*DUP**EXUP + 0.0000014D0) + 540 EXSM = 1.0D0/DBLE(L) + RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0) + RHDN = 0.0D0 + IF (NQ .EQ. 1) GO TO 560 + DDN = VNORM (N, YH(1,L), EWT)/TESCO(1,NQ) + EXDN = 1.0D0/DBLE(NQ) + RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0) + 560 IF (RHSM .GE. RHUP) GO TO 570 + IF (RHUP .GT. RHDN) GO TO 590 + GO TO 580 + 570 IF (RHSM .LT. RHDN) GO TO 580 + NEWQ = NQ + RH = RHSM + GO TO 620 + 580 NEWQ = NQ - 1 + RH = RHDN + IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0 + GO TO 620 + 590 NEWQ = L + RH = RHUP + IF (RH .LT. 1.1D0) GO TO 610 + R = EL(L)/DBLE(L) + DO 600 I = 1,N + 600 YH(I,NEWQ+1) = ACOR(I)*R + GO TO 630 + 610 IALTH = 3 + GO TO 700 + 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610 + IF (KFLAG .LE. -2) RH = DMIN1(RH,0.2D0) +C----------------------------------------------------------------------- +C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS. +C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED. +C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE. +C----------------------------------------------------------------------- + IF (NEWQ .EQ. NQ) GO TO 170 + 630 NQ = NEWQ + L = NQ + 1 + IRET = 2 + GO TO 150 +C----------------------------------------------------------------------- +C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED. +C IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1. +C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE +C YH ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST +C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN +C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED, +C UNTIL IT SUCCEEDS OR H REACHES HMIN. +C----------------------------------------------------------------------- + 640 IF (KFLAG .EQ. -10) GO TO 660 + RH = 0.1D0 + RH = DMAX1(HMIN/DABS(H),RH) + H = H*RH + DO 645 I = 1,N + 645 Y(I) = YH(I,1) + IERR = 0 + CALL F (NEQ, TN, Y, SAVF, IERR) + IF (IERR .LT. 0) RETURN + NFE = NFE + 1 + DO 650 I = 1,N + 650 YH(I,2) = H*SAVF(I) + IPUP = MITER + IALTH = 5 + IF (NQ .EQ. 1) GO TO 200 + NQ = 1 + L = 2 + IRET = 3 + GO TO 150 +C----------------------------------------------------------------------- +C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD +C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP. +C----------------------------------------------------------------------- + 660 KFLAG = -1 + GO TO 720 + 670 KFLAG = -2 + GO TO 720 + 680 KFLAG = -3 + GO TO 720 + 690 RMAX = 10.0D0 + 700 R = 1.0D0/TESCO(2,NQU) + DO 710 I = 1,N + 710 ACOR(I) = ACOR(I)*R + 720 HOLD = H + JSTART = 1 + RETURN +C----------------------- END OF SUBROUTINE STODE ----------------------- + END