V10/cmd/lcc/port.bak/odes2.f
SUBROUTINE ODES2(F,X,NX,TSTART,TSTOP,DT,ERROR,ERRPAR,HANDLE,
1 GLBMAX,ERPUTS,KMAX,MMAX)
C
C TO SOLVE THE INITAL VALUE PROBLEM FOR
C
C DX(T)/DT = F(T,X(T)).
C
C METHOD - RATIONAL EXTRAPOLATION OF GRAGGS MODIFIED MID-POINT RULE.
C
C THE 2 EXTRA ARGUMENTS IN THIS SUBROUTINE PROVIDE USER CONTROL
C OVER THE MAXIMUM ORDER AND MAXIMUM LEVEL OF EXTRAPOLATION TO BE USED
C BY THE PROCEDURE.
C
C INPUT
C
C F - CALL F(T,X,NX,FTX) SHOULD RETURN FTX(I)=F(T,X)(I),
C FOR I=1,...,NX. IF IT CANNOT, IT SHOULD RETURN
C OKAY=.FALSE. IN COMMON /ODESF/OKAY .
C F SHOULD BE DECLARED EXTERNAL IN THE SUBPROGRAM
C CALLING ODES.
C X - THE INITIAL VALUES FOR THE SOLUTION.
C NX - THE LENGTH OF THE SOLUTION VECTOR X.
C TSTART - THE INITIAL TIME.
C TSTOP - THE FINAL TIME.
C DT - THE INITIAL TIME-STEP TO BE USED.
C THE PERFORMANCE OF ODES IS SUBSTANTIALLY
C INDEPENDENT OF THE VALUE OF DT CHOSEN BY THE USER.
C IT IS SUFFICIENT THAT THE USERS CHOICE FOR DT MERELY BE
C WITHIN SEVERAL ORDERS OF MAGNITUDE OF BEING CORRECT.
C THE VALUE OF DT WILL BE AUTOMATICALLY CHANGED BY ODES
C DURING THE INTEGRATION PROCESS IN SUCH A WAY AS TO GET
C THE SOLUTION, TO THE DESIRED ACCURACY, AT THE LEAST
C POSSIBLE COST.
C ERROR - ERROR TOLERANCE ROUTINE.
C
C LOGICAL FUNCTION ERROR(X,NX,T,DT,ERRPAR,ERPUTS,E)
C
C HAS AS INPUT
C
C X,T - X=X(T), THE APPROXIMATE SOLUTION FOR WHICH
C AN ERROR CRITERION IS DESIRED.
C NX - THE LENGTH OF X.
C DT - THE TIME-STEP USED TO OBTAIN X(T).
C ERRPAR - TWO PARAMETERS, AS GIVEN TO ODES2.
C ERPUTS - THIS VARIABLE HAS THE SAME VALUE AS ERPUTS
C IN THE CALL TO ODES2.
C E - THE REAL ABSOLUTE ERROR IN X(I) IS E(I),
C I=1,...,NX, FOR THE SINGLE CURRENT TIME-STEP.
C
C THE OUTPUT IS
C
C ERRPAR - MAY BE ALTERED IF DESIRED.
C E - E(I) IS THE REAL TOLERABLE ABSOLUTE ERROR IN
C X(I), I=1,...,NX. ALL E(I) MUST BE POSITIVE.
C
C FUNCTION VALUE
C
C ERROR - ERROR=.TRUE. IF CONVERGED.
C ERROR=.FALSE. IF NOT.
C ERRPAR - TWO PARAMETERS TO BE PASSED TO ERROR.
C HANDLE - OUTPUT ROUTINE WITH A CALLING SEQUENCE OF THE FORM
C
C HANDLE(T0,X0,T1,X1,NX,DT,TSTOP,E)
C
C HANDLE WILL BE CALLED AT THE END OF EACH TIME-STEP.
C
C THE INPUT TO HANDLE IS AS FOLLOWS
C
C X0,X1 - X0=X(T0) AND X1=X(T1).
C T0,T1 - T0=T1 INDICATES A RESTART AND X1 IS FULL OF
C GARBAGE.
C NX - THE LENGTH OF THE SOLUTION VECTOR X.
C DT - THE PROPOSED TIME-STEP FOR THE NEXT STEP.
C TSTOP - THE CURRENT FINAL TIME.
C E - E(I) GIVES THE REAL ABSOLUTE ERROR IN X1(I),
C I=1,...,NX, FOR THE SINGLE CURRENT TIME-STEP.
C
C THE OUTPUT FROM HANDLE MAY BE ANY OF
C
C X1 - MAY BE ALTERED IF DESIRED.
C DT - THE PROPOSED TIME-STEP FOR THE NEXT STEP.
C TSTOP - THE FINAL TIME VALUE.
C
C HANDLE SHOULD BE DECLARED EXTERNAL IN THE
C SUBPROGRAM CALLING ODES.
C GLBMAX - IF (GLBMAX) THEN THE GLOBAL MAXIMUM ABSOLUTE VALUE OF
C THE SOLUTION IS TO BE RECORDED,
C SEE COMMON /ODESM/ BELOW.
C
C IF THE ERROR SUBPROGRAM SUPPLIED BY THE USER IS ODESE,
C THEN GLBMAX DETERMINES WHETHER OR NOT THE GLOBAL
C MAXIMUM ABSOLUTE VALUE OF THE SOLUTION WILL BE USED IN
C THAT SUBPROGRAM.
C ERPUTS - IF (ERPUTS) THEN THE ERROR PER UNIT TIME-STEP CRITERION
C IS TO BE USED IN THE ERROR ROUTINE.
C OTHERWISE, THE ERROR PER TIME-STEP CRITERION IS TO BE
C USED IN THE ERROR ROUTINE.
C
C IF THE ERROR SUBPROGRAM SUPPLIED BY THE USER IS ODESE,
C THEN ERPUTS DETERMINES WHETHER OR NOT THE ERROR
C PER UNIT-TIME-STEP OR ERROR PER TIME-STEP ERROR
C CRITERION WILL BE USED BY THAT SUBPROGRAM.
C KMAX - THE MAXIMUM NUMBER OF COLUMNS ALLOWED IN THE
C EXTRAPOLATION PROCESS.
C MMAX - THE MAXIMUM NUMBER OF LEVELS OF EXTRAPOLATION PERMITTED.
C MMAX.GE.KMAX+2 IS REQUIRED.
C
C OUTPUT
C
C X - X=X(TSTOP), THE FINAL VALUE FOR THE SOLUTION.
C TSTOP - MAY BE ALTERED BY USER SUPPLIED ROUTINE HANDLE.
C DT - PROPOSED TIME-STEP FOR THE NEXT STEP, IF ANY.
C ERRPAR - MAY BE ALTERED BY USER SUPPLIED ROUTINE ERROR.
C
C SCRATCH SPACE OF LENGTH
C
C S(ODES2) .LE.
C
C 2*MMAX + NX*(KMAX+2)
C
C REAL WORDS +
C
C 5*KMAX + 3*MMAX + 3 +
C
C ( IF (GLBMAX) THEN 2*NX , OTHERWISE 0 ) +
C
C MAX ( 2*NX REAL + S(F) ,
C
C NX*(KMAX+1) +
C
C MAX( KMAX REAL + KMAX , S(ERROR) ) ,
C
C NX + S(HANDLE) )
C
C INTEGER WORDS IS ALLOCATED.
C
C ERROR STATES
C
C 1 - NX.LT.1.
C 2 - DT HAS THE WRONG SIGN ON INPUT.
C 3 - DT=0 ON INPUT.
C 4 - DT RETURNED BY HANDLE HAS THE WRONG SIGN.
C 5 - DT=0 WAS RETURNED BY HANDLE. (RECOVERABLE)
C 6 - THE ERROR DESIRED IN X(I) IS NOT POSITIVE. (RECOVERABLE)
C 7 - DT=0. (RECOVERABLE)
C 8 - CANNOT RAISE DT IN HANDLE WHEN .NOT.OKAY.
C 9 - KMAX.LT.1.
C 10 - KMAX.GT.MMAX-2.
C
C INTERNAL NAMED COMMON USAGE -
C
C ODESM HOLDS THE POINTER, IGMAX, TO THE REAL
C VECTOR OF CURRENT MAXIMUM ABSOLUTE VALUES ATTAINED BY EACH
C COMPONENT OF THE SOLUTION, AND THE POINTER, IGMAXO, TO THE
C REAL VECTOR OF MAXIMUM ABSOLUTE VALUES ATTAINED BY EACH
C COMPONENT OF THE SOLUTION AS OF THE LAST TIME STEP.
C IGMAX=0 MEANS THESE VECTORS ARE NOT USED AND HAVE NOT BEEN
C ALLOCATED.
C A0DESQ HOLDS THE REAL END-POINT VALUE AND THE REAL
C REMAINING ERROR FOR USE IN ODESQ.
C A0DESP HOLDS THE POINTER TO THE REAL VECTOR
C F(T0,X(T0)).
C
COMMON /ODESF/OKAY
COMMON /ODESM/IGMAX,IGMAXO
COMMON /A0DESQ/TEND,RERROR
COMMON /A0DESP/IFTX0
C
REAL X(NX),TSTART,TSTOP,DT,TEND
REAL ERRPAR(2),RERROR
LOGICAL GLBMAX,ERPUTS,OKAY
EXTERNAL F,ERROR,HANDLE
C
REAL DELTA
EXTERNAL A0DESG,A0DES0,A0DESE
C
COMMON /CSTAK/DS
DOUBLE PRECISION DS(500)
REAL RS(1000)
INTEGER IS(1000)
EQUIVALENCE (DS(1),RS(1)),(DS(1),IS(1))
C
IF (TSTART.EQ.TSTOP) GO TO 20
C
CALL ENTER(1)
C
C ... CHECK FOR ERRORS IN THE INPUT.
C
C/6S
C IF (NX.LT.1) CALL SETERR(16H ODES2 - NX.LT.1,16,1,2)
C IF ((DT/ABS(DT))*(TSTOP-TSTART).LT.0.0E0)
C 1 CALL SETERR(39H ODES2 - DT HAS THE WRONG SIGN ON INPUT,39,2,2)
C IF (TSTART+DT.EQ.TSTART)
C 1 CALL SETERR(22H ODES2 - DT=0 ON INPUT,22,3,2)
C IF (KMAX.LT.1) CALL SETERR(18H ODES2 - KMAX.LT.1,18,9,2)
C IF (KMAX.GT.MMAX-2)
C 1 CALL SETERR(23H ODES2 - KMAX.GT.MMAX-2,23,10,2)
C/7S
IF (NX.LT.1) CALL SETERR(' ODES2 - NX.LT.1',16,1,2)
IF ((DT/ABS(DT))*(TSTOP-TSTART).LT.0.0E0)
1 CALL SETERR(' ODES2 - DT HAS THE WRONG SIGN ON INPUT',39,2,2)
IF (TSTART+DT.EQ.TSTART)
1 CALL SETERR(' ODES2 - DT=0 ON INPUT',22,3,2)
IF (KMAX.LT.1) CALL SETERR(' ODES2 - KMAX.LT.1',18,9,2)
IF (KMAX.GT.MMAX-2)
1 CALL SETERR(' ODES2 - KMAX.GT.MMAX-2',23,10,2)
C/
C
C ... ALLOCATE AND LOAD N WITH THE SEQUENCE 1,2,3,4,6,8,12,... .
C
IN=ISTKGT(MMAX,2)
C
DO 10 I=1,MMAX
J=IN+I-1
IF (I.LT.4) IS(J)=I
10 IF (I.GT.3) IS(J)=2*IS(J-2)
C
C ... DECIDE IF HAVE ERROR PER UNIT-TIME-STEP OR PER TIME-STEP.
C
DELTA=0.0E0
IF (ERPUTS) DELTA=1.0E0
C
C ... LOAD THE GLOBAL MAXIMUM ABSOLUTE VALUES OF THE SOLUTION,
C ... IF NECESSARY.
C
IGMAX=0
IF (GLBMAX) IGMAX=ISTKGT(2*NX,3)
IGMAXO=IGMAX+NX
IF (GLBMAX) CALL A0DESL(X,NX,RS(IGMAX),RS(IGMAXO))
C
TEND=TSTOP
RERROR=ERRPAR(2)
C
IFTX0=ISTKGT(NX,3)
C
CALL SXTRP(TSTART,TSTOP,A0DESG,F,1.0E0,2.0E0,DELTA,X,NX,DT,
1 IS(IN),KMAX,MMAX,.FALSE.,ERROR,A0DESE,ERRPAR,HANDLE,
2 A0DES0,0.95E0)
C
IF (NERROR(NERR).NE.0) CALL ERROFF
C/6S
C IF (NERR.EQ.13) CALL SETERR
C 1 (13H ODES2 - DT=0,13,7,1)
C IF (NERR.EQ.14) CALL SETERR
C 1 (32H ODES2 - DT=0 RETURNED BY HANDLE,32,5,1)
C IF (NERR.EQ.17) CALL SETERR
C 1 (50H ODES2 - THE ERROR DESIRED IN X(I) IS NOT POSITIVE,50,6,1)
C/7S
IF (NERR.EQ.13) CALL SETERR
1 (' ODES2 - DT=0',13,7,1)
IF (NERR.EQ.14) CALL SETERR
1 (' ODES2 - DT=0 RETURNED BY HANDLE',32,5,1)
IF (NERR.EQ.17) CALL SETERR
1 (' ODES2 - THE ERROR DESIRED IN X(I) IS NOT POSITIVE',50,6,1)
C/
C
CALL LEAVE
C
20 RETURN
C
END