c This is a modified version of DVERK that uses a defect based error and c stepsize control. It produces a continous approximation to the solution c on each step. The continuous solution can be accessed on each step c by using option C(9) to return control before accepting a step. On c return from DVERK the driver program can check to ensure that the c current step will be accepted, IND=5, and if so the continuous c solution can be evaluated at any point in the current step c (xt between X and XTRIAL) using the routine INTRP. A derivative c evaluation routine, DERIV, is also available for use at this point. c c This routine uses the same calling sequence and shares the same c overview as the original code DVERK. For this reason the internal c documentation for DVERK is relevant to this modified version as c well. In particular it describes the options available and how c they can be invoked. c c Note that W must be declared to be W(n,14) and the value of YTRIAL c (on return after an attempted step) is stored in W(*,13). c C YTRIAL IS IN W(K,13) SUBROUTINE DERIV( N, X, Y, XOUT, YPOUT, HTRIAL, NW, W) C INTEGER N, NW, I C EXTERNAL FCN DOUBLE PRECISION X,Y(N),XOUT,YPOUT(N),W(NW,14),TAU,HTRIAL DOUBLE PRECISION D0,D1,D2,D3,D4,D5,D6 TAU = ( XOUT - X ) / HTRIAL D0 = -40.D0/(17.D0*HTRIAL) * (TAU - 1.D0) + * ((((20.D0*TAU - 40.4D0)*TAU + 23.1D0)*TAU + - 1.7D0)*TAU - .85D0) + - (TAU - 1.D0)**2/(17.D0*HTRIAL) + *(((1.6D3*TAU-2424.D0)*TAU+924.D0)*TAU-34.D0) D1 = -2.D3/8721.D0 * ((3.D0*TAU - 4.D0)*TAU + 1.D0) + * ((( 15.1D0*TAU - 32.542D0)*TAU + + 21.903D0)*TAU - 4.3605D0) + - TAU*(TAU - 1.D0)**2/8721.D0 + *((90600.D0*TAU - 130168.D0)*TAU + 43806.D0) D2 = 40.D0/(17.D0*HTRIAL) * TAU + * (((( 20.D0*TAU - 80.4D0)*TAU + 123.9D0)*TAU + - 88.3D0)*TAU + 25.65D0) + + TAU**2/(17.D0*HTRIAL) + *(((1.6D3*TAU-4824.D0)*TAU+4956.D0)*TAU-1766.D0) D3 = 4.D3/102.D0 * (3.D0*TAU - 2.D0)*TAU + * ((( 1.9D0*TAU - 3.698D0)*TAU + + 2.08D0)*TAU - .2565D0) + + (TAU - 1.D0)*TAU**2 / 51.D0 + * ((114.D2*TAU - 14792.D0)*TAU + 4160.D0) D4 = -8.D3/459.D0 * (TAU*(TAU - 1.D0)*(2.D0*TAU - 1.D0)) + * (( 2.3D0*TAU - 4.136D0)*TAU + 14.877D0/8.D0) + -(TAU*(TAU-1.D0))**2/918.D0 * (368.D2*TAU-33088.D0) D5 = 250.D0/918.D0 * TAU*(TAU - 1.D0)*(2.D0*TAU - 1.D0) + * ((800.D0*TAU - 800.D0)*TAU + 57.D0) + + 125.D0/918.D0*(TAU*(TAU-1.D0))**2 * (1.6D3*TAU-8.D2) D6 = -16.D3/8721.D0 * TAU*(TAU - 1.D0)*(2.D0*TAU - 1.D0) + * ((2.D2*TAU - 2.D2)*TAU + 27.D0) + - 8.D3/8721.D0 * (TAU*(TAU - 1.D0))**2 + * (4.D2*TAU - 2.D2) C DO 40 I = 1,N YPOUT(I) = D0*Y(I) + D1*W(I,1) + D2*W(I,13) + D3*W(I,9) + + D4*W(I,10) + D5*W(I,11) + D6*W(I,12) 40 CONTINUE C RETURN END C SUBROUTINE INTRP( N, X, Y, XOUT, YOUT, HTRIAL, NW, W) C INTEGER N, NW, I DOUBLE PRECISION X,Y(N),XOUT,YOUT(N),W(NW,14),TAU,HTRIAL DOUBLE PRECISION D0,D1,D2,D3,D4,D5,D6 TAU = ( XOUT - X ) / HTRIAL D0 = -20.D0/17.D0 * (TAU-1.D0)**2 + * ((((20.D0*TAU - 40.4D0)*TAU + 23.1D0)*TAU + - 1.7D0)*TAU - .85D0) D1 = -2.D3/8721.D0 * TAU * (TAU-1.D0)**2 + * ((( 15.1D0*TAU - 32.542D0)*TAU + + 21.903D0)*TAU - 4.3605D0) * HTRIAL D2 = 20.D0/17.D0 * TAU**2 + * (((( 20.D0*TAU - 80.4D0)*TAU + 123.9D0)*TAU + - 88.3D0)*TAU + 25.65D0) D3 = 4.D3/102.D0 * (TAU-1.D0) * TAU**2 + * ((( 1.9D0*TAU - 3.698D0)*TAU + + 2.08D0)*TAU - .2565D0) * HTRIAL D4 = -4.D3/459.D0 * TAU**2 * (TAU-1.D0)**2 + * (( 2.3D0*TAU - 4.136D0)*TAU + 14.877D0/8.D0) * HTRIAL D5 = 25.D3/459.D0 * TAU**2 * (TAU-1.D0)**2 + * (( 2.D0*TAU - 2.D0)*TAU + .1425D0) * HTRIAL D6 = -8.D5/8721.D0 * TAU**2 * (TAU-1.D0)**2 + * ((2.D0*TAU - 2.D0)*TAU + .27D0) * HTRIAL C DO 40 I = 1,N YOUT(I) = D0*Y(I) + D1*W(I,1) + D2*W(I,13) + D3*W(I,9) + + D4*W(I,10) + D5*W(I,11) + D6*W(I,12) 40 CONTINUE C RETURN END C c SUBROUTINE DEFECT(N,XOUT,DFCT,FCN,X,Y,C,NW,W) c DOUBLE PRECISION XOUT,DFCT(N),X,Y(N),C(30),W(NW,16),YPP(10), c + YOUT(10),YP(10) c INTEGER N,I,NW c EXTERNAL FCN C C CALL INTRP(N,X,Y,XOUT,YOUT,C(18),NW,W) C CALL FCN( XOUT,YOUT,YP) C CALL DERIV(N,X,Y,XOUT,YPP,C(18),NW,W) C DO 10 I = 1,N C DFCT(I) = YPP(I) - YP(I) c10 CONTINUE c RETURN c END C SUBROUTINE DVERK (N, FCN, X, Y, XEND, TOL, IND, C, NW, W) INTEGER N, IND, NW, K, I DOUBLE PRECISION X, Y(N), XEND, TOL, C(30), W(NW,14), TEMP DOUBLE PRECISION A1, A3, A4, A5, A6, A7, A8, A0, + BY0, BY1, BY5, BF0, BF1, BF5, + CY0, CY1, CY5, CF0, CF1, CF5, + D0, D1, D2, D3, D4, D5, D6 PARAMETER ( A1 = 49.D0/640.D0, A3 = 4375.D0/11968.D0, + A4 = 23.D0/384.D0, A5 = -33.D0/1955.D0, + A6 = 0.D0, A7 = 125.D0/8832.D0, + A8 = -43.D0/1408.D0, A0 = 1.D0/32.D0 ) PARAMETER ( BY0 = 128.D0/3125.D0, BY1 = 2592.D0/3125.D0, + BY5 = 81.D0/625.D0, BF0 = 18.D0/3125.D0, + BF1 = -162.D0/3125.D0, BF5 = 162.D0/3125.D0 ) PARAMETER ( CY0 = 5427.D0/400000.D0, CY1 = 380133.D0/400000.D0, + CY5 = 361.D0/10000.D0, CF0 = 1539.D0/800000.D0, + CF1 = -29241.D0/800000.D0, CF5 = 3249.D0/200000.D0 ) C C ADDED EXTERNAL STATEMENT FOR FCN TO AVOID A WARNING MESSAGE. C EXTERNAL FCN C*********************************************************************** C C ****************************************************************** C * BEGIN INITIALIZATION, PARAMETER CHECKING, INTERRUPT RE-ENTRIES * C ****************************************************************** C C ......ABORT IF IND OUT OF RANGE 1 TO 6 IF (IND.LT.1 .OR. IND.GT.6) GO TO 500 C C CASES - INITIAL ENTRY, NORMAL RE-ENTRY, INTERRUPT RE-ENTRIES GO TO (5, 5, 45, 1111, 2222, 2222), IND C CASE 1 - INITIAL ENTRY (IND .EQ. 1 OR 2) C .........ABORT IF N.GT.NW OR TOL.LE.0 5 IF (N.GT.NW .OR. TOL.LE.0.D0) GO TO 500 IF (IND.EQ. 2) GO TO 15 C INITIAL ENTRY WITHOUT OPTIONS (IND .EQ. 1) C SET C(1) TO C(9) EQUAL TO 0 DO 10 K = 1, 9 C(K) = 0.D0 10 CONTINUE GO TO 35 15 CONTINUE C INITIAL ENTRY WITH OPTIONS (IND .EQ. 2) C MAKE C(1) TO C(9) NON-NEGATIVE DO 20 K = 1, 9 C(K) = DABS(C(K)) 20 CONTINUE C MAKE FLOOR VALUES NON-NEGATIVE IF THEY ARE TO BE USED IF (C(1).NE.4.D0 .AND. C(1).NE.5.D0) GO TO 30 DO 25 K = 1, N C(K+30) = DABS(C(K+30)) 25 CONTINUE 30 CONTINUE 35 CONTINUE C INITIALIZE RREB, DWARF, PREV XEND, FLAG, COUNTS C(10) = 2.D0**(-56) C(11) = 1.D-35 C SET PREVIOUS XEND INITIALLY TO INITIAL VALUE OF X C(20) = X DO 40 K = 21, 24 C(K) = 0.D0 40 CONTINUE GO TO 50 C CASE 2 - NORMAL RE-ENTRY (IND .EQ. 3) C .........ABORT IF XEND REACHED, AND EITHER X CHANGED OR XEND NOT 45 IF (C(21).NE.0.D0 .AND. + (X.NE.C(20) .OR. XEND.EQ.C(20))) GO TO 500 C RE-INITIALIZE FLAG C(21) = 0.D0 GO TO 50 C CASE 3 - RE-ENTRY FOLLOWING AN INTERRUPT (IND .EQ. 4 TO 6) C TRANSFER CONTROL TO THE APPROPRIATE RE-ENTRY POINT.......... C THIS HAS ALREADY BEEN HANDLED BY THE COMPUTED GO TO . C END CASES V 50 CONTINUE C C END INITIALIZATION, ETC. C C ***************************************************************** C * LOOP THROUGH THE FOLLOWING 4 STAGES, ONCE FOR EACH TRIAL STEP * C * UNTIL THE OCCURRENCE OF ONE OF THE FOLLOWING * C * (A) THE NORMAL RETURN (WITH IND .EQ. 3) ON REACHING XEND IN * C * STAGE 4 * C * (B) AN ERROR RETURN (WITH IND .LT. 0) IN STAGE 1 OR STAGE 4 * C * (C) AN INTERRUPT RETURN (WITH IND .EQ. 4,5 OR 6), IF * C * REQUESTED, IN STAGE 1 OR STAGE 4 * C ****************************************************************** C 99999 CONTINUE C C *************************************************************** C * STAGE 1 - PREPARE - DO CALCULATIONS OF HMIN, HMAX, ETC., * C * AND SOME PARAMETER CHECKING, AND END UPWITH SUITABLE * C * VALUES OF HMAG, XTRIAL AND HTRIAL IN PREPARATION FOR TAKING * C * AN INTEGRATION STEP. * C *************************************************************** C C***********ERROR RETURN (WITH IND=-1) IF NO OF FCN EVALS TOO GREAT IF (C(7).EQ.0.D0 .OR. C(24).LT.C(7)) GO TO 100 IND = -1 RETURN 100 CONTINUE C C STORE K1 IN THE FIRST COLLUMN OF W. UNLESS THIS IS THE FIRST C STEP THIS VALUE HAS ALREADY BEEN COMPUTED. IF (IND .EQ. 6) GO TO 105 IF (IND .EQ. 5) THEN DO 102 K = 1,N W(K,1) = W(K,9) 102 CONTINUE ELSE CALL FCN(N,X,Y,W(1,1)) C(24) = C(24) + 1 ENDIF 105 CONTINUE C C CALCULATE HMIN - USE DEFAULT UNLESS VALUE PRESCRIBED C(13) = C(3) IF (C(3) .NE. 0.D0) GO TO 165 C CALCULATE DEFAULT VALUE OF HMIN C FIRST CALCULATE WEIGHTED NORM Y - C(12) - AS SPECIFIED C BY THE ERROR CONTROL INDICATOR C(1) TEMP = 0.D0 IF (C(1) .NE. 1.D0) GO TO 115 C ABSOLUTE ERROR CONTROL - WEIGHTS ARE 1 DO 110 K = 1, N TEMP = DMAX1(TEMP, DABS(Y(K))) 110 CONTINUE C(12) = TEMP GO TO 160 115 IF (C(1) .NE. 2.D0) GO TO 120 C RELATIVE ERROR CONTROL - WEIGHTS ARE 1/DABS(Y(K)) SO C WEIGHTED NORM Y IS 1 C(12) = 1.D0 GO TO 160 120 IF (C(1) .NE. 3.D0) GO TO 130 C WEIGHTS ARE 1/MAX(C(2),ABS(Y(K))) DO 125 K = 1, N TEMP = DMAX1(TEMP, DABS(Y(K))/C(2)) 125 CONTINUE C(12) = DMIN1(TEMP, 1.D0) GO TO 160 130 IF (C(1) .NE. 4.D0) GO TO 140 C WEIGHTS ARE 1/MAX(C(K+30),ABS(Y(K))) DO 135 K = 1, N TEMP = DMAX1(TEMP, DABS(Y(K))/C(K+30)) 135 CONTINUE C(12) = DMIN1(TEMP, 1.D0) GO TO 160 140 IF (C(1) .NE. 5.D0) GO TO 150 C WEIGHTS ARE 1/C(K+30) DO 145 K = 1, N TEMP = DMAX1(TEMP, DABS(Y(K))/C(K+30)) 145 CONTINUE C(12) = TEMP GO TO 160 150 CONTINUE C DEFAULT CASE - WEIGHTS ARE 1/MAX(1,ABS(Y(K))) DO 155 K = 1, N TEMP = DMAX1(TEMP, DABS(Y(K))) 155 CONTINUE C(12) = DMIN1(TEMP, 1.D0) 160 CONTINUE C(13) = 10.D0*DMAX1(C(11),C(10)*DMAX1(C(12)/TOL,DABS(X))) 165 CONTINUE C C CALCULATE SCALE - USE DEFAULT UNLESS VALUE PRESCRIBED C(15) = C(5) IF (C(5) .EQ. 0.D0) C(15) = 1.D0 C C CALCULATE HMAX - CONSIDER 4 CASES C CASE 1 BOTH HMAX AND SCALE PRESCRIBED IF (C(6).NE.0.D0 .AND. C(5).NE.0.D0) + C(16) = DMIN1(C(6), 2.D0/C(5)) C CASE 2 - HMAX PRESCRIBED, BUT SCALE NOT IF (C(6).NE.0.D0 .AND. C(5).EQ.0.D0) C(16) = C(6) C CASE 3 - HMAX NOT PRESCRIBED, BUT SCALE IS IF (C(6).EQ.0.D0 .AND. C(5).NE.0.D0) C(16) = 2.D0/C(5) C CASE 4 - NEITHER HMAX NOR SCALE IS PROVIDED IF (C(6).EQ.0.D0 .AND. C(5).EQ.0.D0) C(16) = 2.D0 C C***********ERROR RETURN (WITH IND=-2) IF HMIN .GT. HMAX IF (C(13) .LE. C(16)) GO TO 170 IND = -2 RETURN 170 CONTINUE C C CALCULATE PRELIMINARY HMAG - CONSIDER 3 CASES IF (IND .GT. 2) GO TO 175 C CASE 1 - INITIAL ENTRY - USE PRESCRIBED VALUE OF HSTART, IF C ANY, ELSE DEFAULT C(14) = C(4) IF (C(4) .EQ. 0.D0) C(14) = C(16)*TOL**(1./6.) GO TO 185 175 IF (C(23) .GT. 1.D0) GO TO 180 C CASE 2 - AFTER A SUCCESSFUL STEP, OR AT MOST ONE FAILURE, C USE MIN(2, .9*(TOL/EST)**(1/6))*HMAG, BUT AVOID POSSIBLE C OVERFLOW. THEN AVOID REDUCTION BY MORE THAN HALF. TEMP = 2.D0*C(14) IF (TOL .LT. (2.D0/.9D0)**6*C(19)) + TEMP = .9D0*(TOL/C(19))**(1./6.)*C(14) C(14) = DMAX1(TEMP, .5D0*C(14)) GO TO 185 180 CONTINUE C CASE 3 - AFTER TWO OR MORE SUCCESSIVE FAILURES C(14) = .5D0*C(14) 185 CONTINUE C C CHECK AGAINST HMAX C(14) = DMIN1(C(14), C(16)) C C CHECK AGAINST HMIN C(14) = DMAX1(C(14), C(13)) C C***********INTERRUPT NO 1 (WITH IND=4) IF REQUESTED IF (C(8) .EQ. 0.D0) GO TO 1111 IND = 4 RETURN C RESUME HERE ON RE-ENTRY WITH IND .EQ. 4 ........RE-ENTRY.. 1111 CONTINUE C C CALCULATE HMAG, XTRIAL - DEPENDING ON PRELIMINARY HMAG, XEND IF (C(14) .GE. DABS(XEND - X)) GO TO 190 C DO NOT STEP MORE THAN HALF WAY TO XEND C(14) = DMIN1(C(14), .5D0*DABS(XEND - X)) C(17) = X + DSIGN(C(14), XEND - X) GO TO 195 190 CONTINUE C HIT XEND EXACTLY C(14) = DABS(XEND - X) C(17) = XEND 195 CONTINUE C C CALCULATE HTRIAL C(18) = C(17) - X C C END STAGE 1 C C *************************************************************** C * STAGE 2 - CALCULATE YTRIAL (ADDING 10 TO NO OF FCN EVALS).* C * W(*,2), ... W(*,12) HOLD INTERMEDIATE RESULTS NEEDED IN * C * STAGE 3. W(*,13) IS TEMPORARY STORAGE UNTIL FINALLY IT HOLDS* C * YTRIAL. * C *************************************************************** C TEMP = C(18)/1398169080000.D0 C DO 200 K = 1, N W(K,13) = Y(K) + TEMP*W(K,1)*233028180000.D0 200 CONTINUE CALL FCN(N, X + C(18)/6.D0, W(1,13), W(1,2)) C DO 205 K = 1, N W(K,13) = Y(K) + TEMP*(W(K,1)*74569017600.D0 + + W(K,2)*298276070400.D0) 205 CONTINUE CALL FCN(N, X + C(18)*(4.D0/15.D0), W(1,13), W(1,3)) C DO 210 K = 1, N W(K,13) = Y(K) + TEMP*(W(K,1)*1165140900000.D0 + - W(K,2)*3728450880000.D0 + + W(K,3)*3495422700000.D0 ) 210 CONTINUE CALL FCN(N, X + C(18)*(2.D0/3.D0), W(1,13), W(1,4)) C DO 215 K = 1, N W(K,13) = Y(K) + TEMP*( - W(K,1)*3604654659375.D0 + + W(K,2)*12816549900000.D0 + - W(K,3)*9284716546875.D0 + + W(K,4)*1237962206250.D0 ) 215 CONTINUE CALL FCN(N, X + C(18)*(5.D0/6.D0), W(1,13), W(1,5)) C C DO 220 K = 1, N C W(K,13) = Y(K) + TEMP*(W(K,1)*3355605792000.D0 C + - W(K,2)*11185352640000.D0 C + + W(K,3)*9172628850000.D0 C + - W(K,4)*427218330000.D0 C + + W(K,5)*482505408000.D0) C 220 CONTINUE C CALL FCN(N, X + C(18), W(1,13), W(1,6)) C DO 225 K = 1, N W(K,13) = Y(K) + TEMP*( - W(K,1)*770204740536.D0 + + W(K,2)*2311639545600.D0 + - W(K,3)*1322092233000.D0 + - W(K,4)*453006781920.D0 + + W(K,5)*326875481856.D0) 225 CONTINUE CALL FCN(N, X + C(18)/15.D0, W(1,13), W(1,7)) C DO 230 K = 1, N W(K,13) = Y(K) + TEMP*(W(K,1)*2845924389000.D0 + - W(K,2)*9754668000000.D0 + + W(K,3)*7897110375000.D0 + - W(K,4)*192082660000.D0 + + W(K,5)*400298976000.D0 + + W(K,7)*201586000000.D0) 230 CONTINUE CALL FCN(N, X + C(18), W(1,13), W(1,8)) C C CALCULATE YTRIAL, THE EXTRAPOLATED APPROXIMATION AND STORE C IN W(*,13) DO 235 K = 1, N W(K,13) = Y(K) + TEMP*(W(K,1)*104862681000.D0 + + W(K,3)*545186250000.D0 + + W(K,4)*446637345000.D0 + + W(K,5)*188806464000.D0 + + W(K,7)*15076875000.D0 + + W(K,8)*97599465000.D0) 235 CONTINUE C C NOW CALCULATE K9,K10 AND Y5 (STORED IN W(*,2)) C CALL FCN (N, C(17), W(1,13), W(1,9) ) C DO 240 I = 1,N W(I,2) = Y(I) + C(18) * ( A1*W(I,1) + A3*W(I,3) + A4*W(I,4) + + A5*W(I,5) + A7*W(I,7) + A8*W(I,8) + + A0*W(I,9) ) 240 CONTINUE C CALL FCN (N, X + .5D0*C(18), W(1,2), W(1,10) ) C CNOW USE THESE VALUES TO CALCULATE K11 AND K12 USING C W(*,14) AS TEMPORARY STROAGE. C DO 250 I = 1,N W(I,14) = BY0*Y(I) + BY1*W(I,13) + BY5*W(I,2) + + C(18) * ( BF0*W(I,1) + BF1*W(I,9) + BF5*W(I,10)) 250 CONTINUE C CALL FCN (N, X + .90D0*C(18), W(1,14), W(1,11) ) C DO 260 I = 1,N W(I,14) = CY0*Y(I) + CY1*W(I,13) + CY5*W(I,2) + + C(18) * ( CF0*W(I,1) + CF1*W(I,9) + CF5*W(I,10)) 260 CONTINUE C CALL FCN (N, X + .95D0*C(18), W(1,14), W(1,12) ) C C ADD 11 TO THE NO OF FCN EVALS C(24) = C(24) + 10.D0 C C END STAGE 2 C C *************************************************************** C * STAGE 3 - CALCULATE THE ERROR ESTIMATE EST. FIRST CALCULATE * C * THE UNWEIGHTED ABSOLUTE ERROR ESTIMATE VECTOR AND STORE * C * IT IN W(*,2). THEN CALCULATE THE WEIGHTED MAX NORM OF W(*,2)* C * AS SPECIFIED BY THE ERROR CONTROL INDICATOR C(1). * C *************************************************************** C C CALCULATE THE DEFECT AT X + .5*H AND USE IT AS AN ESTIMA C OF THE MAXIMUM DEFECT ON THIS STEP. TEMP = .1943D0 CALL INTRP(N,X,Y,X + TEMP*C(18),W(1,2),C(18),NW,W) CALL FCN(N,X + TEMP*C(18),W(1,2),W(1,14)) CALL DERIV(N,X,Y,X + TEMP*C(18),W(1,2),C(18),NW,W) DO 300 K = 1,N W(K,2) = W(K,2) - W(K,14) 300 CONTINUE C(24) = C(24) + 1.D0 C C CALCULATE THE WEIGHTED MAX NORM OF W(*,2) AS SPECIFIED BY C THE ERROR CONTROL INDICATOR C(1) TEMP = 0.D0 IF (C(1) .NE. 1.D0) GO TO 310 C ABSOLUTE ERROR CONTROL DO 305 K = 1, N TEMP = DMAX1(TEMP,DABS(W(K,2))) 305 CONTINUE GO TO 360 310 IF (C(1) .NE. 2.D0) GO TO 320 C RELATIVE ERROR CONTROL DO 315 K = 1, N TEMP = DMAX1(TEMP, DABS(W(K,2)/Y(K))) 315 CONTINUE GO TO 360 320 IF (C(1) .NE. 3.D0) GO TO 330 C WEIGHTS ARE 1/MAX(C(2),ABS(Y(K))) DO 325 K = 1, N TEMP = DMAX1(TEMP, DABS(W(K,2)) + / DMAX1(C(2), DABS(Y(K))) ) 325 CONTINUE GO TO 360 330 IF (C(1) .NE. 4.D0) GO TO 340 C WEIGHTS ARE 1/MAX(C(K+30),ABS(Y(K))) DO 335 K = 1, N TEMP = DMAX1(TEMP, DABS(W(K,2)) + / DMAX1(C(K+30), DABS(Y(K))) ) 335 CONTINUE GO TO 360 340 IF (C(1) .NE. 5.D0) GO TO 350 C WEIGHTS ARE 1/C(K+30) DO 345 K = 1, N TEMP = DMAX1(TEMP, DABS(W(K,2)/C(K+30))) 345 CONTINUE GO TO 360 350 CONTINUE C DEFAULT CASE - WEIGHTS ARE 1/MAX(1,ABS(Y(K))) DO 355 K = 1, N TEMP = DMAX1(TEMP, DABS(W(K,2)) + / DMAX1(1.D0, DABS(Y(K))) ) 355 CONTINUE 360 CONTINUE C C CALCULATE EST - (THE WEIGHTED MAX NORM OF W(*,2))*HMAG*SCALE C(19) = TEMP C C END STAGE 3 C C *************************************************************** C * STAGE 4 - MAKE DECISIONS. * C *************************************************************** C C SET IND=5 IF STEP ACCEPTABLE, ELSE SET IND=6 IND = 5 IF (C(19) .GT. TOL) IND = 6 C C***********INTERRUPT NO 2 IF REQUESTED IF (C(9) .EQ. 0.D0) GO TO 2222 RETURN C RESUME HERE ON RE-ENTRY WITH IND .EQ. 5 OR 6 ...RE-ENTRY.. 2222 CONTINUE C IF (IND .EQ. 6) GO TO 410 C STEP ACCEPTED (IND .EQ. 5), SO UPDATE X, Y FROM XTRIAL, C YTRIAL, ADD 1 TO THE NO OF SUCCESSFUL STEPS, AND SET C THE NO OF SUCCESSIVE FAILURES TO ZERO X = C(17) DO 400 K = 1, N Y(K) = W(K,13) 400 CONTINUE C(22) = C(22) + 1.D0 C(23) = 0.D0 C**************RETURN(WITH IND=3, XEND SAVED, FLAG SET) IF X .EQ. XEND IF (X .NE. XEND) GO TO 405 IND = 3 C(20) = XEND C(21) = 1.D0 RETURN 405 CONTINUE GO TO 420 410 CONTINUE C STEP NOT ACCEPTED (IND .EQ. 6), SO ADD 1 TO THE NO OF C SUCCESSIVE FAILURES C(23) = C(23) + 1.D0 C**************ERROR RETURN (WITH IND=-3) IF HMAG .LE. HMIN IF (C(14) .GT. C(13)) GO TO 415 IND = -3 RETURN 415 CONTINUE 420 CONTINUE C C END STAGE 4 C GO TO 99999 C END LOOP C C BEGIN ABORT ACTION 500 CONTINUE C WRITE(6,505) IND, TOL, X, N, C(13), XEND, NW, C(16), C(20), + C(22), C(23), C(24), (Y(K), K = 1, N) 505 FORMAT( /// 1H0, 58HCOMPUTATION STOPPED IN DVERK WITH THE FOLLOWIN +G VALUES - + / 1H0, 5HIND =, I4, 5X, 6HTOL =, 1PD13.6, 5X, 11HX =, +1PD22.15 + / 1H , 5HN =, I4, 5X, 6HHMIN =, 1PD13.6, 5X, 11HXEND =, +1PD22.15 + / 1H , 5HNW =, I4, 5X, 6HHMAX =, 1PD13.6, 5X, 11HPREV XEND =, +1PD22.15 + / 1H0, 14X, 27HNO OF SUCCESSFUL STEPS = , 0PF8.0 + / 1H , 14X, 27HNO OF SUCCESSIVE FAILURES =, 0PF8.0 + / 1H , 14X, 27HNO OF FUNCTION EVALS = , 0PF8.0 + / 1H0, 23HTHE COMPONENTS OF Y ARE + // (1H , 1P5D24.15) ) C STOP C C END ABORT ACTION C END C