C----------------------------------------------------------------------- C C THE FOUR USER CALLABLE ROUTINES DOCUMENTED BELOW ARE : C C 1. RKNSET - SETUP ROUTINE - THE USER MUST CALL THIS ROUTINE PRIOR TO C THE FIRST CALL OF THE CORE INTEGRATOR TO SET INTEGRATION C TOLERANCES AND OPTIONAL INPUTS. C C 2. RKNINT - CORE INTEGRATION ROUTINE C C 3. RKNDIA - DIAGNOSTIC ROUTINE - THE USER MAY CALL THIS ROUTINE TO C EXTRACT OPTIONAL OUTPUT ABOUT THE INTEGRATION. C C 4. EVAL - INTERPOLATION ROUTINE - THE USER MAY CALL THIS ROUTINE C ONLY WHEN THE LOWER ORDER FORMULAS ARE EMPLOYED. C C C IT IS LIKELY THAT THE CODE WILL BE INCLUDED IN THE D02L SUBCHAPTER C AT MARK 13 OF THE NAG FORTRAN LIBRARY - NUMERICAL ALGORITHMS GROUP, C 256 BANBURY ROAD, OXFORD, ENGLAND, OX2 7DE. C C THE CODE IS AVAILABLE FROM THE AUTHORS ON IBM PC XT FORMAT FLOPPY C DISKETTE. C C C----------------------------------------------------------------------- C C SUBROUTINE RKNSET(NEQ,H,TOL,THRES,THRESP,MAXSTP,START, C + ONESTP,HIGH,LWORK,RWORK,IFAIL) C C*********************************************************************** C ABSTRACT C*********************************************************************** C C THIS IS THE SETUP ROUTINE FOR THE NYSTROM CODE RKNINT. IT MUST BE C CALLED BEFORE THE FIRST CALL TO RKNINT FOR ANY PROBLEM AND CAN BE C USED TO CHANGE OPTIONAL INPUTS BETWEEN CONTINUATION CALLS. C C THE MACHINE-DEPENDENT CONSTANT UROUND IS SET IN RKNSET AND PASSED C IN THE WORKING STORAGE ARRAY RWORK TO OTHER ROUTINES. UROUND IS C DEFINED AS THE LARGEST POSITIVE NUMBER SUCH THAT C UROUND = 1 + UROUND IN THE ARITHMETIC BEING USED. C C ASSOCIATED ROUTINES ARE RKNINT, THE NYSTROM INTEGRATOR, THE C DIAGNOSTIC ROUTINE RKNDIA AND THE DENSE OUTPUT ROUTINE, EVAL. C C*********************************************************************** C OVERVIEW OF ARGUMENTS C*********************************************************************** C C NEQ - NUMBER OF DIFFERENTIAL EQUATIONS TO BE SOLVED. C C H - OPTIONAL STEPSIZE. C C TOL - ERROR CONTROL PARAMETER. C C THRES, THRESP - OPTIONAL ERROR CONTROL PARAMETERS. C C MAXSTP - OPTIONAL MAXIMUM NUMBER OF ATTEMPTED STEPS. C C START - LOGICAL VARIABLE INDICATING INITIALIZATION. C C ONESTP - LOGICAL VARIABLE INDICATING MODE IN WHICH RKNINT IS TO RUN. C C HIGH - LOGICAL VARIABLE INDICATING WHICH FORMULA PAIR IS TO C BE USED. C C LWORK - DIMENSION OF RWORK IN THE CALLING PROGRAM. C C RWORK - ARRAY USED TO PASS OPTIONAL PARAMETERS TO RKNINT. C C IFAIL - ON EXIT, INDICATES STATUS OF RKNSET. C C*********************************************************************** C INITIALIZATION CALL TO RKNSET C*********************************************************************** C C TO INITIALIZE RKNINT FOR THE START OF A NEW PROBLEM, CALL RKNSET WITH C START = .TRUE. C C NEQ - INTEGER. NUMBER OF DIFFERENTIAL EQUATIONS. MUST BE SPECIFIED. C C H - REAL. IF H .EQ. 0 THEN RKNINT WILL SELECT A SUITABLE C STARTING STEPSIZE C ELSE ABS(H) IN THE APPROPRIATE DIRECTION C WILL BE USED. C C THE FIRST STEP TAKEN BY RKNINT IS A CRITICAL ONE BECAUSE IT MUST C REFLECT HOW FAST Y AND Y' CHANGE NEAR THE INITIAL POINT. THE USER C CAN PROVIDE RKNINT WITH SOME IDEA OF THE SCALE OF THE PROBLEM VIA C THIS PARAMETER - ANYTHING REASONABLE WILL DO. A MISTAKE IN THE C SIGN OF THE GUESSED H WILL BE FORGIVEN AS ONLY ITS SIZE IS C IMPORTANT. C C**** IF YOU DO NOT HAVE AN IDEA AS TO THE SCALE OF THE PROBLEM C**** SET H = 0 ON THE INITIAL CALL TO RKNSET AND RKNINT WILL C**** GENERATE A SUITABLE VALUE OF H AT SOME ADDITIONAL EXPENSE. C C TOL - REAL. C THE USER MUST SET THIS ERROR PARAMETER TO SPECIFY HOW C ACCURATELY HE WANTS THE CODE TO COMPUTE THE SOLUTION. C THE RELATIVE ERROR TOLERANCE TOL MUST SATISFY C C 1 .GE. TOL .GE. 1O * UROUND. C C TOL AND THE VALUES IN THE ARRAYS THRES(*) AND THRESP(*) C ARE USED IN THE LOCAL ERROR TEST. AT EACH STEP THIS C TEST REQUIRES ROUGHLY THAT C C ABS(LOCAL ERROR IN Y) .LE. TOL * MAX(ABS(Y),THRES) C C FOR EACH COMPONENT OF Y AND THRES, AND C C ABS(LOCAL ERROR IN YP) .LE. TOL * MAX(ABS(YP),THRESP) C C FOR EACH COMPONENT OF YP AND THRESP. THE DEFAULT VALUE FOR C EACH COMPONENT OF THRES AND THRESP IS THRDEF, CURRENTLY C DEFINED TO BE 50*UROUND. THIS CHOICE YIELDS ESSENTIALLY C A PURE RELATIVE ERROR TEST BUT THE USE OF THE SMALL POSITIVE C THRES(HOLD) WILL HELP TO AVOID TROUBLE SHOULD A SOLUTION C COMPONENT VANISH ON A PARTICULAR STEP. NOTE THAT THERE IS C NO DEFAULT VALUE FOR TOL. C C A NEARLY PURE RELATIVE ERROR TEST MAY BE TOO DEMANDING C FOR COMPONENTS WHICH BECOME VERY SMALL. OFTEN COMPONENTS C WHICH BECOME SMALLER IN MAGNITUDE THAN A THRES(HOLD) ARE NOT C PHYSICALLY INTERESTING. A SUITABLE THRES(HOLD) VALUE DEPENDS C THE SCALING OF THE PROBLEM AND HENCE REQUIRES THE USER'S INSI C A PRELIMINARY INTEGRATION MIGHT BE NECESSARY TO LEARN ENOUGH C ABOUT THE PROBLEM TO SELECT SUITABLE ERROR PARAMETERS FOR THE C DEFINITIVE INTEGRATION. C C THE CODE WILL NOT ATTEMPT TO COMPUTE A SOLUTION AT AN ACCURAC C UNREASONABLE FOR THE COMPUTER BEING USED. IT WILL ADVISE YOU C SHOULD THIS SITUATION ARISE. TO CONTINUE THE INTEGRATION, YO C WILL HAVE TO INCREASE TOL AND/OR THRES AND/OR THRESP. C C PRACTICALLY ALL PRESENT DAY CODES, INCLUDING THIS ONE, CONTRO C AN ESTIMATE OF THE LOCAL ERROR AT EACH STEP AND DO NOT EVEN C ATTEMPT TO CONTROL THE GLOBAL ERROR DIRECTLY. IN PARTICULAR, C CONTROLLING THE LOCAL ERROR IN A RELATIVE SENSE DOES NOT NECE C RESULT IN CONTROL OF THE GLOBAL ERROR IN A RELATIVE SENSE. C USUALLY, BUT NOT ALWAYS, THE TRUE ACCURACY OF THE APPROXIMATE C SOLUTION IS COMPARABLE TO THE ERROR TOLERANCES. ALSO, THIS C CODE WILL USUALLY, BUT NOT ALWAYS, YIELD A MORE ACCURATE C SOLUTION IF THE TOLERANCES ARE REDUCED AND THE INTEGRATION IS C REPEATED. BY COMPARING TWO SUCH SOLUTIONS IT IS POSSIBLE TO C OBTAIN A FAIRLY RELIABLE IDEA OF THE TRUE ERROR IN THE C SOLUTION AT THE LARGER TOLERANCES. C C THRES - REAL ARRAY OF SIZE .GE. NEQ, USED IN ERROR CONTROL FOR C COMPONENTS OF Y. C IF THRES(1) .LE. 0 THEN RKNSET USES THE DEFAULT (THRDEF) FOR C ALL COMPONENTS OF THRES. C ELSE IT CHECKS EACH COMPONENT. C IF THRES(K) .LE. 0 THEN IFAIL = 1 C ELSE ACCEPT THRES(K C C THRESP - REAL ARRAY OF SIZE .GE. NEQ, USED IN ERROR CONTROL FOR C COMPONENTS OF YP. THRESP IS DEFINED IN THE SAME WAY AS THRES C C MAXSTP - MAXIMUM NUMBER OF ATTEMPTED STEPS ALLOWED ON **ANY CALL** TO C RKNINT. IF MAXSTP .LE. 0 THEN A DEFAULT OF 1000 IS USED C ELSE MAXSTP IS USED. C C START - LOGICAL. SET .TRUE. FOR INITIALIZATION CALL. C C ONESTP - LOGICAL. SET .TRUE. IF YOU WISH TO RUN RKNINT IN STEP-ORIENT C MODE (CONTROL RETURNS TO THE CALLING PROGRAM AFTER EACH STEP) C SET .FALSE. FOR INTERVAL MODE (INTEGRATION TO TEND BEFORE C CONTROL RETURNS TO CALLING PROGRAM). C C HIGH - LOGICAL. SET .TRUE. IF YOU WISH TO USE THE 12(10) FORMULA C PAIR AND .FALSE. IF YOU WISH TO USE THE 6(4) PAIR. C C THE 12(10) PAIR REQUIRES MORE STORAGE THAN THE 6(4) PAIR. THE C HIGH ORDER PAIR WILL PROBABLY TAKE LESS WORK IF THE ERROR C REQUIREMENT IS STRINGENT WHILE THE LOWER ORDER PAIR WILL USUALL C BE MORE EFFICIENT AT RELAXED TOLERANCES. IF THE DENSE OUTPUT C ROUTINE EVAL IS TO BE USED, THE 6(4) PAIR MUST BE SELECTED. C C LWORK - INTEGER. DIMENSION OF RWORK IN THE CALLING PROGRAM. C IF (HIGH) THEN LWORK MUST BE .GE. 14 + 20*NEQ C ELSE LWORK MUST BE .GE. 14 + 11*NEQ. C C RWORK - REAL ARRAY OF SIZE LWORK. USED TO PASS OPTIONAL C INPUTS TO RKNINT. THE SAME ARRAY MUST BE PASSED TO RKNINT C AND MUST NOT BE CHANGED BY THE USER. C C IFAIL - INTEGER. USED FOR OUTPUT ONLY. C C*********************************************************************** C ON EXIT FROM RKNSET: C*********************************************************************** C C **** START ALWAYS HAS THE VALUE .FALSE. ON RETURN FROM RKNSET. C C IFAIL = 0 == SUCCESSFUL INITIALIZATION. C C IFAIL = 1 == THE USER INDICATED THAT HE WISHED TO SET THRES (THRESP) C BY SETTING THRES(1) (THRESP(1)) .GT. 0.0 BUT C AT LEAST ONE COMPONENT OF THRES (THRESP) WAS .LE. 0.0. C C IFAIL = 2 == LWORK NOT LARGE ENOUGH. C C IFAIL = 3 == (TOL .GT. 1) OR (TOL .LT. 10*UROUND). C C*********************************************************************** C CHANGING OPTIONAL INPUTS BEFORE A CONTINUATION CALL TO RKNINT C*********************************************************************** C C NEQ - MUST NOT BE CHANGED. C C H - IF H .EQ. 0 THEN THE PREDICTED STEPSIZE CALCULATED IN C RKNINT WILL BE ATTEMPTED ON THE NEXT STEP C ELSE RKNINT WILL ATTEMPT A STEP OF SIZE C ABS(H) IN THE APPROPRIATE DIRECTION. C C TOL, THRES, THRESP - THE VALUES OF TOL, THRES AND THRESP MAY BE CHANGE C CHECKING IS DONE AS ON AN INITIALIZATION CALL. C C MAXSTP - CAN BE CHANGED. SET AS ON AN INITIALIZATION CALL. C C START - MUST BE .FALSE. THE USER NEED NOT RESET START AS THE C VALUE .FALSE. WAS RETURNED BY THE INITIALIZATION CALL. C C ONESTP - CAN BE CHANGED. SET AS ON AN INITIALIZATION CALL. C C HIGH - CAN BE CHANGED. YOU MUST SUPPLY SUFFICIENT STORAGE FOR C THE FORMULA PAIR SELECTED. C C LWORK - CHECKED ON EACH CALL TO RKNSET. C C RWORK - MUST NOT BE CHANGED BY THE USER. C C ON EXIT, POSSIBLE VALUES OF IFAIL ARE THE SAME AS FOR THE C INITIALIZATION CALL. C C*********************************************************************** C SUBROUTINE RKNINT(F,NEQ,T,TEND,Y,YP,YDP,RWORK,IFAIL) C C ********************************************************************** C ** ABSTRACT ** C ************** C C SUBROUTINE RKNINT USES RUNGE-KUTTA-NYSTROM METHODS TO C INTEGRATE A SYSTEM OF SECOND ORDER ORDINARY DIFFERENTIAL C EQUATIONS OF THE FORM C Y'' = F(T,Y). C TWO SETS OF FORMULAS ARE PROVIDED - A 6(4) (LOW-ORDER) PAIR C THAT USES 6 STAGES AND A 12(10) (HIGH-ORDER) PAIR THAT USES C 17 STAGES. THE USER SELECTS THE PAIR TO BE USED. C C THE 6(4) FORMULAS ARE DESCRIBED IN C DORMAND,J.R., PRINCE,P.J. (1986) RUNGE-KUTTA NYSTROM TRIPLES. C REPORT TP-CS-86-05, DEPT. OF COMPUTER SCIENCE, TEESSIDE C POLYTECHNIC, MIDDLESBROUGH, CLEVELAND. C C THE 12(10) FORMULAS ARE DESCRIBED IN C DORMAND,J.R., EL-MIKKAWY,M.E.A., PRINCE,P.J. (1986) HIGH ORDER C EMBEDDED RUNGE-KUTTA-NYSTROM FORMULAE. C REPORT TP-MR-86-02, DEPT. OF MATHEMATICS, TEESSIDE C POLYTECHNIC, MIDDLESBROUGH, CLEVELAND. C C GIVEN INITIAL VALUES FOR Y AND Y' AT AN INITIAL TIME T, THE C ROUTINE INTEGRATES FROM T TO TEND. WITH THE SOLUTION Y AND C Y' OF THAT PROBLEM, THE USER CAN RESET TEND AND CONTINUE THE C INTEGRATION. IT IS EASY TO MAKE THE ROUTINE RETURN THE C INTERMEDIATE VALUES OF Y AND Y' AT EACH STEP ON THE WAY TO TEND. C C A DENSE OUTPUT FACILITY IS PROVIDED WITH THE LOW-ORDER PAIR, C THROUGH A THIRD FORMULA OF ORDER 6 WHICH USES THE SAME STAGES. C RKNINT RETURNS WITH SUFFICIENT INFORMATION TO ALLOW C THE USER TO CALL A SEPARATE EVALUATION ROUTINE WHICH CAN C COMPUTE APPROXIMATIONS TO Y AND Y' AT ANY T VALUE IN THE LAST C COMPLETED STEP. C C RKNINT CALLS THE ROUTINES COEF64 AND COEFHI. C ASSOCIATED ROUTINES ARE THE SETUP ROUTINE RKNSET, THE DIAGNOSTIC C ROUTINE RKNDIA AND THE OUTPUT ROUTINE EVAL. C C*********************************************************************** C** OVERVIEW OF THE ARGUMENTS FOR RKNINT C ********************************************************************** C C F - SUBROUTINE SUPPLIED BY THE USER WHICH DESCRIBES THE SYSTEM OF C DIFFERENTIAL EQUATIONS C C NEQ - NUMBER OF DIFFERENTIAL EQUATIONS C C T - CURRENT VALUE OF THE INDEPENDENT VARIABLE C C TEND - THE CODE WILL ATTEMPT TO INTEGRATE TO TEND C C Y, YP - ARRAYS CONTAINING COMPONENTS OF THE SOLUTION AND ITS C DERIVATIVE AT T C C YDP - AFTER INITIALIZATION, THIS ARRAY CONTAINS F(T,Y). C C RWORK - A WORKING STORAGE ARRAY. THE SETUP ROUTINE RKNSET PASSES C INFORMATION TO RKNINT USING RWORK. THE ARRAY MUST BE OF SI C A) FOR THE 6(4) PAIR, .GE. 14 + 11*NEQ C B) FOR THE 12(10) PAIR, .GE. 14 + 20*NEQ. C WHEN THE 6(4) PAIR IS USED, THIS ARRAY CONTAINS THE C INFORMATION NECESSARY TO COMPUTE THE INTERPOLANT. C C IFAIL - INTEGER REPORTING THE STATUS OF THE INTEGRATION. C C*********************************************************************** C** HOW TO INITIALIZE RKNINT C*********************************************************************** C C *** BEFORE CALLING RKNINT FOR THE FIRST TIME IN AN INTEGRATION, THE C *** SETUP ROUTINE RKNSET MUST BE CALLED. USING RKNSET, THE USER C *** ALLOCATES SUFFICIENT STORAGE FOR THE VARIABLES USED BY THE CODE, C *** SELECTS THE FORMULA PAIR, MODE OF INTEGRATION AND ACCURACY C *** REQUIREMENT AND SETS APPROPRIATE VALUES FOR INITIALIZATION OF C *** RKNINT. C *** SEE THE DOCUMENTATION FOR RKNSET. C C BEFORE CALLING THE SUBROUTINE RKNINT, THE USER MUST DEFINE THE C PROBLEM TO BE SOLVED. C C F - PROVIDE A SUBROUTINE OF THE FORM C F(NEQ,T,Y,YDP) C TO DEFINE THE SYSTEM OF DIFFERENTIAL EQUATIONS, WHERE Y AND YDP C ARE DIMENSIONED APPROPRIATELY IN F. GIVEN THE VALUES OF T AND C Y, THE ROUTINE F MUST EVALUATE THE NEQ COMPONENTS OF Y'' AND C RETURN THE VALUES IN YDP. C C F MUST NOT ALTER T OR Y. THE NAME F MUST BE DECLARED EXTERNAL C IN THE CALLING PROGRAM. C C NEQ - INTEGER. NUMBER OF DIFFERENTIAL EQUATIONS. C C T - REAL. DEFINES THE INITIAL POINT OF THE INTEGRATION. C C TEND - REAL. RKNINT WILL ATTEMPT TO INTEGRATE FROM T TO TEND. C INTEGRATION EITHER FORWARD IN T (TEND .GT. T) OR BACKWARD C IN T (TEND .LT. T) IS PERMITTED. C C IF T .EQ. TEND, NO ACTION WILL BE TAKEN. C C THE CODE ADVANCES THE SOLUTION FROM T TO TEND USING STEPSIZE C THAT ARE AUTOMATICALLY SELECTED SO AS TO ACHIEVE THE REQUEST C ACCURACY AS EFFICIENTLY AS POSSIBLE. C IF INDICATED, THE CODE WILL RETURN WITH THE SOLUTION C FOLLOWING EACH INTERMEDIATE STEP (STEP-ORIENTED MODE) BUT TE C MUST STILL BE PROVIDED. C C NOTE THAT RKNINT NEVER INTEGRATES PAST TEND - IT WILL HIT TE C EXACTLY. THIS FEATURE IS USEFUL IF THERE ARE KNOWN VALUES O C SUCH THAT INTEGRATION PAST T SHOULD NOT BE ATTEMPTED. C C Y - REAL ARRAY OF SIZE .GE. NEQ. ON THE FIRST CALL, SET TO THE VAL C OF Y AT THE INITIAL POINT. C C YP - REAL ARRAY OF SIZE .GE. NEQ. ON THE FIRST CALL, SET TO THE VA C OF Y' AT THE INITIAL POINT. C C YDP - REAL ARRAY OF SIZE .GE. NEQ. USED ONLY FOR OUTPUT. C C RWORK - REAL. WORKING STORAGE ARRAY. C C IFAIL - INTEGER. ERROR INDICATOR. MUST BE SET TO ZERO ON ENTRY. C C*********************************************************************** C** OUTPUT - AFTER ANY EXIT FROM RKNINT C*********************************************************************** C C THE PRINCIPAL AIM OF THE CODE IS TO RETURN A COMPUTED SOLUTION AT T C ALTHOUGH IT IS POSSIBLE TO OBTAIN INTERMEDIATE RESULTS ALONG THE WA C TO FIND OUT WHETHER THE CODE ACHIEVED ITS GOAL OR IF THE INTEGRATIO C WAS INTERRUPTED, CHECK T AND IFAIL. C C T - THE INTEGRATION WAS SUCCESSFULLY ADVANCED TO THE OUTPUT VALUE C OF T. C C Y, YP - CONTAIN THE APPROXIMATE SOLUTION AND DERIVATIVE AT T. C C YDP - CONTAINS THE VALUE OF F AT T, Y. C C RWORK - CONTAINS THE INTERMEDIATE STAGES FROM THE LAST STEP C AS WELL AS THE VALUES OF Y, YP AND YDP AT THE START C OF THE LAST STEP. THE USER MUST NOT ALTER THIS ARRAY C BEFORE A CALL TO THE DENSE OUTPUT ROUTINE, EVAL. IN C PARTICULAR, A CALL TO RKNSET SHOULD NOT BE MADE BETWEEN C A CALL TO RKNINT AND A CALL TO EVAL. C C IFAIL - REPORTS THE PROGRESS OF THE CODE. C C 0 = SUCCESSFUL INTEGRATION TO THE VALUE GIVEN BY T. C C > 0 = FAILURE OF SOME SORT. C 1 = ON ANY CALL, T AND TEND WERE THE SAME ON INPUT. C ON A CONTINUATION CALL, EITHER NEQ OR THE DIRECTION C OF INTEGRATION WAS CHANGED. C 2 = THE NUMBER OF ATTEMPTED STEPS WAS GREATER THAN THE C MAXIMUM ALLOWED. C 3 = THE LOCAL ERROR TEST SPECIFIED BY TOL, THRES AND C THRESP IS TOO STRINGENT. C 4 = THE CODE IS BEING USED INEFFICIENTLY. THE STEP SIZE HAS C BEEN REDUCED BY A SIGNIFICANT AMOUNT TOO OFTEN IN ORDER C TO HIT THE OUTPUT POINTS SPECIFIED BY TEND. C C*********************************************************************** C** CONTINUATION CALLS C*********************************************************************** C C THIS CODE IS ORGANIZED SO THAT SUBSEQUENT CALLS TO CONTINUE THE C INTEGRATION INVOLVE LITTLE (IF ANY) ADDITIONAL EFFORT ON THE USER'S C PART. YOU MUST MONITOR IFAIL IN ORDER TO DETERMINE WHAT TO DO NEXT C C DO NOT ALTER ANY QUANTITY NOT SPECIFICALLY PERMITTED BELOW, IN C PARTICULAR DO NOT CHANGE NEQ, RWORK OR THE DIFFERENTIAL EQUATION C IN SUBROUTINE F. ANY SUCH ALTERATION CONSTITUTES A NEW PROBLEM C AND MUST BE TREATED AS SUCH, I.E., YOU MUST START AFRESH. C C OPTIONAL INPUTS SUCH AS THE FORMULA PAIR, TOL, THRES, C THRESP AND H CAN BE CHANGED BY A CALL TO RKNSET (SEE SECTION ON C OPTIONAL INPUTS AND OUTPUTS). NOTE THAT C INCREASING A TOLERANCE MAKES THE EQUATION EASIER TO INTEGRATE. C DECREASING A TOLERANCE WILL MAKE IT HARDER AND MAY NOT BE A REASONA C CHOICE ON A CONTINUATION CALL. C C ** FOLLOWING A COMPLETED TASK (IFAIL = 0) ** C CALL THE CODE AGAIN TO CONTINUE THE INTEGRATION TO A NEW C VALUE OF TEND. C C ** FOLLOWING AN INTERRUPTED TASK ** C IF C IFAIL = 1, THERE IS AN ERROR IN A PARAMETER. CHECK INPUT VALUES. C C IFAIL = 2, THE CODE HAS ATTEMPTED MAXSTP STEPS. YOU SHOULD C CONSIDER THE POSSIBILITY THAT THE ERROR PARAMETERS C ARE NOT APPROPRIATE. IN PARTICULAR, IF A SOLUTION C COMPONENT VANISHES AND THE CORRESPONDING COMPONENT C OF THRES OR THRESP IS TOO SMALL, THE CODE MAY NEED A C OF STEPS. YOU MAY ALTER TOL, THRES OR THRESP IF YOU C WISH. IF YOU WANT TO CONTINUE YOU **MUST** SET IFAIL C CALL RKNINT AGAIN AND AN ADDITIONAL MAXSTP STEPS WILL C BE ALLOWED. C C IFAIL = 3, THE TEST FOR LIMITING PRECISION IN THE COMPUTER ARITHMET C BEING USED HAS BEEN VIOLATED. IT MAY BE POSSIBLE C TO CONTINUE THE INTEGRATION IF THE ERROR PARAMETERS C TOL, THRES, THRESP ARE INCREASED. C USE RKNSET TO CHANGE TOL, THRES, THRESP. C IF YOU ARE SURE YOU WANT TO CONTINUE WITH RELAXED C ERROR PARAMETERS CALL RKNINT AGAIN WITH THE NEW C VALUES AND YOU **MUST** SET IFAIL=0. C C IFAIL = 4, THE CODE HAS DETECTED INEFFICIENT USE OF THE ROUTINE BY C THE USER'S CALLING (SUB)PROGRAM. OUT OF THE LAST 100 C OR MORE SUCCESSFUL STEPS, MORE THAN TEN PER CENT ARE C STEPS WITH SIZES THAT HAVE HAD TO REDUCED BY A FACTOR C GREATER THAN 0.5 IN ORDER TO HIT THE USER SPECIFIED C POINTS TEND. TO CONTINUE YOU **MUST** SET IFAIL=0 C AND CALL RKNINT AGAIN. C C **** IF RKNINT RETURNS TO THE CALLING PROGRAM WITH IFAIL .GT. 0 C **** YOU MUST TAKE SOME ACTION TO CORRECT THE ERROR BEFORE RKNINT C **** IS CALLED AGAIN. OTHERWISE AN INFINITE LOOP MAY OCCUR. THE C **** ROUTINE WILL ATTEMPT TO DETECT THE CASE WHERE AN ERROR STATE C **** IS NOT CORRECTED. IF TWO FAILURES OCCUR AT THE SAME VALUE OF C **** T (IN A SINGLE INTEGRATION) RKNINT WILL PRINT A MESSAGE C **** AND STOP. THE ROUTINE WILL ALSO STOP IF THE INPUT VALUE OF C **** IFAIL IS NONZERO. C C*********************************************************************** C OPTIONAL INPUTS AND OUTPUTS C*********************************************************************** C C THE SETUP ROUTINE RKNSET CAN BE USED TO PASS OPTIONAL VALUES C TO RKNINT. THE ROUTINE MUST BE CALLED BEFORE THE FIRST CALL TO C RKNINT FOR A PARTICULAR PROBLEM AND CAN BE CALLED BETWEEN C CONTINUATION CALLS IF YOU WISH TO CHANGE ANY OPTIONAL VALUES. C SEE THE DOCUMENTATION FOR RKNSET. C C THE OPTIONAL INPUTS ARE: C C H - STEP SIZE FOR THE NEXT STEP. C C TOL, THRES, THRESP - LOCAL ERROR CONTROL PARAMETERS. C C MAXSTP - MAXIMUM NUMBER OF ATTEMPTED STEPS ALLOWED ON ANY C CALL TO RKNINT. C C ONESTP - MODE OF INTEGRATION. C C HIGH - CHOICE OF FORMULA PAIR. C C C THE DIAGNOSTIC ROUTINE RKNDIA CAN BE CALLED AFTER ANY CALL TO C RKNINT TO RETURN INFORMATION REGARDING THE PERFORMANCE OF THE C INTEGRATOR. C C OPTIONAL OUTPUTS: C C HSTART - THE STEPSIZE ACTUALLY USED BY RKNINT ON THE FIRST STEP. C C HUSED - THE STEPSIZE LAST USED BY RKNINT. C C HNEXT - THE STEPSIZE PREDICTED BY RKNINT FOR USE ON THE NEXT STEP. C C NSUCC - THE NUMBER OF SUCCESSFUL STEPS TAKEN BY RKNINT SINCE LAST C INITIALIZATION. C C NFAIL - THE NUMBER OF FAILED STEPS SINCE LAST INITIALIZATION. C C NATT - THE NUMBER OF ATTEMPTED STEPS TAKEN TO FIND AN APPROPRIATE C INITIAL STEPSIZE. C C*********************************************************************** C SUBROUTINE RKNDIA(NEQ, HNEXT, HUSED, HSTART, NSUCC, NFAIL, NATT, C + THRES, THRESP, RWORK) C C********************************************************************** C ABSTRACT C********************************************************************** C C RKNDIA IS THE DIAGNOSTIC ROUTINE ASSOCIATED WITH THE NYSTROM SOLVER, C RKNINT. AFTER A CALL TO RKNINT, THIS ROUTINE CAN BE CALLED TO C OBTAIN INFORMATION ABOUT THE PERFORMANCE OF RKNINT AND THE CURRENT C SETTINGS OF OPTIONAL VALUES. THIS INFORMATION IS SAVED IN RWORK C BY RKNINT AND EXTRACTED BY RKNDIA. C C THE USER SHOULD BE AWARE THAT AFTER AN UNSUCCESSFUL CALL TO RKNINT, C SOME OF THE INFORMATION RETURNED BY RKNDIA MAY REFER TO THE LAST C SUCCESSFUL STEP. C C ASSOCIATED ROUTINES ARE THE NYSTROM INTEGRATOR, RKNINT, AND THE C SETUP ROUTINE, RKNSET. C C********************************************************************** C TO CALL RKNDIA C********************************************************************** C C NEQ - INTEGER. ON INPUT, NUMBER OF DIFFERENTIAL EQUATIONS BEING SOLVED C UNCHANGED ON OUTPUT. C C HNEXT - REAL. OUTPUT ONLY. THE PREDICTED STEPSIZE CALCULATED BY C RKNINT FOR USE ON THE NEXT STEP. C C HUSED - REAL. OUTPUT ONLY. THE STEPSIZE LAST USED BY RKNINT. C C HSTART - REAL. OUTPUT ONLY. THE STEPSIZE USED BY RKNINT FOR THE C FIRST STEP OF THE CURRENT PROBLEM. IT MAY BE DIFFERENT C THAN THE VALUE INPUT BY THE USER. THIS VALUE COULD BE C USEFUL AS AN INITIAL ESTIMATE IN STARTING INTEGRATION C OF PROBLEMS OF SIMILAR TYPE. C C NSUCC - INTEGER. OUTPUT ONLY. NUMBER OF SUCCESSFUL STEPS SINCE C LAST INITIALIZATION OF RKNINT. C C NFAIL - INTEGER. OUTPUT ONLY. NUMBER OF FAILED STEPS SINCE LAST C INITIALIZATION OF RKNINT. C C NATT - INTEGER. OUTPUT ONLY. NUMBER OF ATTEMPTED STEPS TAKEN TO C FIND AN APPROPRIATE STARTING STEPSIZE (WHEN THE USER DOES C NOT GIVE AN ESTIMATE). C THE COST OF AN ATTEMPTED STEP OF THIS SORT WILL BE C APPROXIMATELY HALF THAT OF A USUAL INTEGRATION STEP. C C NOTE THAT NFAIL DOES NOT INCLUDE NATT. A PARTIAL STARTING STEP DOES C NOT COUNT AS A FAILED STEP. C C THRES, THRESP - REAL ARRAYS OF SIZE .GE. NEQ. OUTPUT ONLY. THRESHOLD C TOLERANCES USED IN THE ERROR CONTROL STRATEGY FOR C Y AND YP, RESPECTIVELY. C C RWORK - REAL ARRAY. MUST BE THE SAME AS C THE WORKING STORAGE ARRAY PASSED TO RKNSET AND RKNINT. THE C INFORMATION PASSED BACK BY RKNDIA IS EXTRACTED FROM RWORK. C C*********************************************************************** C SUBROUTINE EVAL(NEQ, T, Y, YP, NWANT, TWANT, YWANT, C + YPWANT, RWORK, IFAIL) C C*********************************************************************** C ABSTRACT C*********************************************************************** C C THIS IS THE DENSE OUTPUT ROUTINE ASSOCIATED WITH THE NYSTROM C INTEGRATOR, RKNINT. EVAL USES THE THIRD FORMULA OF THE RUNGE-KUTTA- C NYSTROM 6(4)6 TRIPLE TO EVALUATE Y AND Y' AT THE POINT SIGMA * H, C WHERE H IS THE STEPSIZE LAST USED BY RKNINT AND SIGMA SHOULD C USUALLY BE IN THE INTERVAL (0, 1). SIGMA IS COMPUTED AS C SIGMA = (TWANT - TOLD) / H. C C EVAL CALLS THE ROUTINE CFEVAL2. C C*********************************************************************** C OVERVIEW OF ARGUMENTS C*********************************************************************** C C NEQ - NUMBER OF DIFFERENTIAL EQUATIONS. C C T - CURRENT VALUE OF THE DEPENDENT VARIABLE, AS RETURNED BY RKNINT. C C Y, YP - VALUES OF Y AND Y' AT T, AS RETURNED BY RKNINT. C C NWANT - THE USER WISHES INTERMEDIATE VALUES FOR THE FIRST C NWANT COMPONENTS ONLY (NWANT .LE. NEQ). C C TWANT - THE USER WISHES INTERMEDIATE VALUES AT TWANT. C C YWANT, YPWANT - ON OUTPUT, VALUES OF Y AND Y' AT TWANT. C C RWORK - WORKING STORAGE ARRAY USED BY RKNINT. THE INTERMEDIATE C STAGES USED IN THE NYSTROM INTERPOLATION FORMULA ARE C PASSED FROM RKNINT IN THIS ARRAY, ALONG WITH AUXILIARY C INFORMATION. C C IFAIL - STATUS OF ROUTINE. C C*********************************************************************** C TO CALL EVAL C*********************************************************************** C C NEQ - INTEGER. NUMBER OF DIFFERENTIAL EQUATIONS. SAME AS THE VALUE C PASSED TO RKNINT. C C T - REAL. VALUE OF THE INDEPENDENT VARIABLE LAST RETURNED BY RKNINT. C C Y, YP - REAL ARRAYS OF SIZE .GE. NEQ. SHOULD BE PASSED UNCHANGED C FROM THE VALUES RETURNED BY RKNINT. C C NWANT - INTEGER .LE. NEQ. VALUES OF Y AND Y' FOR THE FIRST NWANT C COMPONENTS WILL BE RETURNED. C C TWANT - REAL. INTERPOLATION AT TWANT IS REQUIRED. C C YWANT, YPWANT - REAL ARRAYS OF SIZE .GE. NWANT. USED FOR OUTPUT. C C RWORK - REAL ARRAY. SAME AS THE WORKING C STORAGE ARRAY USED BY RKNINT. SHOULD NOT BE ALTERED BEFORE C A CALL TO EVAL. C C IFAIL - INTEGER. USED FOR OUTPUT ONLY. C C*********************************************************************** C ON EXIT FROM EVAL C*********************************************************************** C C YWANT AND YPWANT WILL CONTAIN APPROXIMATIONS TO Y AND Y' AT TWANT, C PROVIDED IFAIL IS 0 OR 1. OTHERWISE, NO ACTION IS TAKEN. C C IFAIL = 0 == SUCCESS. C C = 1 == EXTRAPOLATION WAS PERFORMED, I.E., THE VALUE OF TWANT C WAS NOT BETWEEN T AND T-H, WHERE H IS THE LAST STEPSIZE C USED BY RKNINT. THE USER SHOULD BE AWARE THAT THE C ERROR IN AN EXTRAPOLATED VALUE MAY BE ** LARGE **. C C = 2 == NWANT .GT. NEQ. EVAL CAN BE CALLED AGAIN DIRECTLY C WHEN NWANT IS RESET. C C = 3 == RKNINT LAST USED THE 12(10) FORMULA PAIR TO INTEGRATE C THE DIFFERENTIAL EQUATIONS. INTERPOLATION IS NOT C AVAILABLE WITH THIS PAIR. C C*********************************************************************** SUBROUTINE RKNSET(NEQ,H,TOL,THRES,THRESP,MAXSTP,START,ONESTP,HIGH, * LWORK,RWORK,IFAIL) C C THIS IS THE SET-UP ROUTINE FOR THE NYSTROM CODE RKNINT. IT C MUST BE CALLED BEFORE THE FIRST CALL TO RKNINT AND CAN BE C USED BETWEEN CONTINUATION CALLS TO CHANGE OPTIONAL INPUTS. C C THE ROUTINE USES THE MACHINE-DEPENDENT CONSTANTS UROUND C AND THRDEF, WHICH IS CURRENTLY SET TO 50*UROUND. C C IF STORAGE IS AT PREMIUM, THEN IN THE CALLING (SUB)PROGRAM THE DUMMY C ARGUMENT 'THRES' MAY BE THE ACTUAL ARGUMENT 'RWORK(15+4*NEQ)' AND C SIMILARLY 'THRESP' MAY BE 'RWORK(15+5*NEQ)'. C C .. Parameters .. DOUBLE PRECISION UROUND, THRDEF, ONE, ZERO, U10 INTEGER NOVHD, NSTMAX PARAMETER (UROUND=3.0D-17,THRDEF=50.0D0*UROUND,ONE=1.0D0, * ZERO=0.0D0,U10=10.0D0*UROUND,NOVHD=14, * NSTMAX=1000) C .. C .. Scalar Arguments .. DOUBLE PRECISION H, TOL INTEGER IFAIL, LWORK, MAXSTP, NEQ LOGICAL HIGH, ONESTP, START C .. C .. Array Arguments .. DOUBLE PRECISION RWORK(LWORK), THRES(NEQ), THRESP(NEQ) C .. C .. Local Scalars .. INTEGER K, NTHR, NTHRP C .. C .. Intrinsic Functions .. INTRINSIC ABS, DBLE C .. C .. Executable Statements .. C IF ((HIGH .AND. LWORK.LT.NOVHD+20*NEQ) * .OR. (LWORK.LT.NOVHD+11*NEQ)) THEN IFAIL = 2 GO TO 100 END IF C C NOTE: ONE == TRUE, ZERO == FALSE C NTHR = NOVHD + 4*NEQ NTHRP = NOVHD + 5*NEQ IF (START) THEN RWORK(7) = ONE RWORK(1) = ABS(H) RWORK(2) = ZERO RWORK(3) = ZERO RWORK(4) = ZERO RWORK(5) = ZERO RWORK(11) = ZERO RWORK(13) = UROUND RWORK(14) = ZERO ELSE RWORK(7) = ZERO IF (H.NE.ZERO) RWORK(3) = ABS(H) END IF IF (ONESTP) THEN RWORK(8) = ONE ELSE RWORK(8) = ZERO END IF IF (MAXSTP.LE.0) THEN RWORK(6) = DBLE(NSTMAX) ELSE RWORK(6) = DBLE(MAXSTP) END IF IF (HIGH) THEN RWORK(10) = ONE ELSE RWORK(10) = ZERO END IF IF ((TOL.GT.ONE) .OR. (TOL.LT.U10)) THEN IFAIL = 3 GO TO 100 ELSE RWORK(12) = TOL END IF C C SET THRES INTO RWORK C IF (THRES(1).LE.ZERO) THEN DO 20 K = 1, NEQ RWORK(NTHR+K) = THRDEF 20 CONTINUE ELSE CVEC CVEC THE FOLLOWING LOOP WAS NOT VECTORIZABLE ON A CDC CYBER 205 DUE TO CVEC THE 'GO TO' STATEMENT. CVEC DO 40 K = 1, NEQ IF (THRES(K).LE.ZERO) THEN IFAIL = 1 GO TO 100 ELSE RWORK(NTHR+K) = THRES(K) END IF 40 CONTINUE END IF C C SET THRESP INTO RWORK C IF (THRESP(1).LE.ZERO) THEN DO 60 K = 1, NEQ RWORK(NTHRP+K) = THRDEF 60 CONTINUE ELSE CVEC CVEC THE FOLLOWING LOOP IS NOT VECTORIZABLE ON A CDC CYBER 205 DUE TO CVEC THE 'GO TO' STATEMENT CVEC DO 80 K = 1, NEQ IF (THRESP(K).LE.ZERO) THEN IFAIL = 1 GO TO 100 ELSE RWORK(NTHRP+K) = THRESP(K) END IF 80 CONTINUE END IF C RWORK(9) = ONE IFAIL = 0 C 100 CONTINUE START = .FALSE. C END SUBROUTINE RKNINT(F,NEQ,T,TEND,Y,YP,YDP,RWORK,IFAIL) C C C PARAMETER STATEMENTS. C THE FIRST SET RELATES TO THE STEP SIZE SELECTION ALGORITHM. C THE SECOND SET PROVIDES INTEGER CONSTANTS AND A UNIT NUMBER C FOR THE OUTPUT MESSAGE WHEN CODE TERMINATES EXECUTION. C THE THIRD SET CONTAINS DOUBLE PRECISION CONSTANTS USED C IN THE CODE. C C C NOTE - C THE LOCAL ARRAYS ARE LARGE ENOUGH TO ACCOMODATE THE COEFFICIENTS C OF THE 12(10) FORMULA PAIR. THIS ADDS CLARITY TO THE PROGRAM BUT C SOME STORAGE COULD BE SAVED IF NECESSARY (BY ADDITIONAL PROGRAMMING C EFFORT). TWO POSSIBLE WAYS OF SAVING STORAGE ARE: C 1) STORE THE MATRIX OF COEFFICIENTS A IN STRICTLY LOWER C TRIANGULAR FORM. THIS WOULD SAVE 153 LOCATIONS. C 2) INCLUDE THE STORAGE REQUIREMENT FOR THE COEFFICIENTS C IN THE WORKSPACE ARRAY RWORK. C C C IF THE USE OF COMMON STATEMENTS IS TO BE AVOIDED, /SAVE1/ CAN C BE REPLACED BY CALLING THE APPROPRIATE COEFFICIENT ROUTINE C (COEFHI OR COEF64) ON EACH ENTRY TO RKNINT AND SETTING THE ARRAY PTR C APPROPRIATELY. SIMILARLY, THE VALUES IN /SAVE3/ COULD BE RECOMPUTED C ON EACH ENTRY. OTHER MEANS WILL HAVE TO BE FOUND TO REPLACE /SAVE2/. C C .. Parameters .. DOUBLE PRECISION R, RMN1, RMN2, RMN3, SCALEL, SCALEH PARAMETER (R=10.0D0,RMN1=1.0D0/R,RMN2=1.0D0/(R*R), * RMN3=1.0D0/(R*R*R),SCALEL=0.8D0,SCALEH=0.75D0) INTEGER IDEV, UROUND, NOVHD, HUNDRD PARAMETER (IDEV=6,UROUND=13,NOVHD=14,HUNDRD=100) DOUBLE PRECISION ONE, ZERO, FIFTH, TWO, HALF, QUAR, ELVNTH, FOUR, * EIGHT, TEN, TWENTY, TWOHUN, TENPC PARAMETER (ONE=1.0D0,ZERO=0.0D0,FIFTH=0.2D0,TWO=2.0D0, * HALF=0.5D0,QUAR=0.25D0,ELVNTH=1.0D0/11.0D0, * FOUR=4.0D0,EIGHT=8.0D0,TEN=1.0D1,TWENTY=2.0D1, * TWOHUN=2.0D2,TENPC=0.1D0) C .. C .. Scalar Arguments .. DOUBLE PRECISION T, TEND INTEGER IFAIL, NEQ C .. C .. Array Arguments .. DOUBLE PRECISION RWORK(*), Y(NEQ), YDP(NEQ), YP(NEQ) C .. C .. Subroutine Arguments .. EXTERNAL F C .. C .. Scalars in Common .. DOUBLE PRECISION DIR, EXPON, H, RANGE, RHIGH, SCALE, STBRAD, * TCHECK, TOL, TOOSML INTEGER EFFCT1, EFFCT2, FINPH2, ISTP, NASTEP, NATT, * NBASE, NEQOLD, NFLSTP, NRT1, NRT2, NRT3, NRT4, * NSTAGE, NSTMAX, NTHR, NTHRP, NWT, NWTP LOGICAL FAILED, FIRST, HGIVEN, LAST, LOOP, OPTRED, * PHASE2, SUCCES C .. C .. Arrays in Common .. DOUBLE PRECISION A(17,17), B(16), BHAT(15), BP(16), BPHAT(17), * C(17) INTEGER PTR(16) C .. C .. Local Scalars .. DOUBLE PRECISION ABSH, ABSHP, ALPHA, BETA, DEL, ERR, ERRP, FDEL, * HSQ, HUSED, SCALE1, SCL, SCLP, SK, SPK, TAU, * THRES, THRESP, TOLD, TRY, WT, WTP, YDPDEL, * YDPNRM, YINTK, YPDEL, YPINTK, YPNRM INTEGER I, J, JSTAGE, K LOGICAL DFLT, DFLTP, INEFFC, WASTE C .. C .. External Subroutines .. EXTERNAL COEF64, COEFHI C .. C .. Intrinsic Functions .. INTRINSIC ABS, DBLE, INT, MAX, MIN, SIGN C .. C .. Common blocks .. COMMON /SAVE1/A, C, BPHAT, B, BP, BHAT, PTR COMMON /SAVE2/DIR, RANGE, TCHECK, H, EXPON, RHIGH, TOL, * SCALE, STBRAD, TOOSML, NEQOLD, NSTAGE, NSTMAX, * NASTEP, ISTP, NFLSTP, NATT, EFFCT1, EFFCT2, * FINPH2, LAST, FAILED, FIRST, HGIVEN, SUCCES, * OPTRED, PHASE2, LOOP COMMON /SAVE3/NRT1, NRT2, NRT3, NRT4, NBASE, NTHR, * NTHRP, NWT, NWTP C .. C .. Save statement .. SAVE /SAVE1/, /SAVE2/, /SAVE3/ C .. C .. Executable Statements .. C C BLOCK A. C CHECK FOR NONZERO ENTRY VALUE OF IFAIL. C IF (IFAIL.NE.0) THEN WRITE (UNIT=IDEV,FMT=99998) IFAIL STOP END IF NASTEP = 0 IFAIL = 0 C IF (RWORK(7).EQ.ONE) THEN C C BLOCK B. C FIRST CALL IN AN INTEGRATION. C EXTRACT DATA FROM SETUP ROUTINE RKNSET. C TOLD = T H = RWORK(1) HUSED = ZERO RWORK(1) = ZERO NSTMAX = INT(RWORK(6)+FIFTH) RWORK(9) = ZERO RWORK(7) = ZERO TOL = RWORK(12) C C INITIALIZATION. C EFFCT1 = 0 EFFCT2 = 0 ISTP = 0 NFLSTP = 0 NATT = 0 FIRST = .TRUE. LOOP = .FALSE. IF (T.EQ.TEND) THEN IFAIL = 1 GO TO 620 END IF NEQOLD = NEQ RANGE = ABS(TEND-T) DIR = SIGN(ONE,TEND-T) C C THE NEXT THREE LINES RELATE TO THE CODE CHOICE OF H0. C OPTRED = .FALSE. HGIVEN = H .NE. ZERO PHASE2 = FIRST .AND. .NOT. HGIVEN C C SET UP POINTERS TO RWORK. C NRT1 = NOVHD NRT2 = NOVHD + NEQ NRT3 = NOVHD + 2*NEQ NRT4 = NOVHD + 3*NEQ NTHR = NOVHD + 4*NEQ NTHRP = NOVHD + 5*NEQ NBASE = NOVHD + 6*NEQ C C PICK UP CONSTANTS DEPENDING ON ORDER OF METHOD SELECTED AND SET C POINTERS. THE LOCAL ARRAY PTR CONTAINS POINTERS TO WHERE THE STAGES C ARE HELD IN THE ARRAY RWORK. THIS IS TO ALLOW SAVINGS IN STORAGE FOR C 12(10) PAIR AND BOTH PAIRS TO BE IMPLEMENTED IN THE SAME BLOCK OF C CODE. C RHIGH = RWORK(10) IF (RWORK(10).EQ.ZERO) THEN EXPON = FIFTH STBRAD = FOUR TOOSML = TWENTY*RWORK(UROUND) NSTAGE = 6 SCALE = SCALEL CALL COEF64(A,C,BPHAT,B,BP) FINPH2 = 4 NWT = NBASE + 3*NEQ NWTP = NBASE + 4*NEQ PTR(1) = 0 PTR(2) = 1 PTR(3) = 2 PTR(4) = 3 PTR(5) = 4 ELSE EXPON = ELVNTH STBRAD = EIGHT TOOSML = TWOHUN*RWORK(UROUND) NSTAGE = 17 SCALE = SCALEH CALL COEFHI(A,C,BHAT,BPHAT,B,BP) FINPH2 = 14 NWT = NBASE + 12*NEQ NWTP = NBASE + 13*NEQ PTR(1) = 0 PTR(2) = 1 PTR(3) = 2 PTR(4) = 0 PTR(5) = 3 PTR(6) = 4 PTR(7) = 5 PTR(8) = 6 PTR(9) = 7 PTR(10) = 8 PTR(11) = 9 PTR(12) = 10 PTR(13) = 11 PTR(14) = 12 PTR(15) = 13 PTR(16) = 3 END IF C CALL F(NEQ,T,Y,YDP) C IF ( .NOT. HGIVEN) THEN C C PHASE 1 OF CALCULATION OF AN INITIAL STEPSIZE H0 BY THE CODE C IF || YP(T0) || .NE. 0 THEN H0 = (TOL**EXPON) / ||YP(T0)|| C ELSE H0 = TEND - T0 C IF || YDP(T0) || .NE. 0 THEN HP0 = (TOL**EXPON) / ||YDP(T0)|| C ELSE HP0 = TEND - T0 C H = MIN(H0, HP0) C NOTE: 1) EXPON = 1/5 FOR THE 6(4) PAIR, C EXPON = 1/11 FOR THE 12(10) PAIR. C 2) THE STRATEGY GIVEN ABOVE HAS BEEN MODIFIED TO TRY TO C AVOID SELECTING TOO SMALL A VALUE FOR H0 WHEN THE C INITIAL VALUES OF Y AND YP ARE SMALL. C YPNRM = ZERO YDPNRM = ZERO DFLT = .FALSE. DFLTP = .FALSE. CVEC CVEC THE FOLLOWING LOOP WAS NOT VECTORIZABLE ON A CRAY 1S. CVEC THE FOLLOWING LOOP WAS NOT VECTORIZABLE ON A CDC CYBER 205. CVEC C C EVALUATE NORMS AS DESCRIBED ABOVE AND INITIALISE THE WEIGHT VECTORS C FOR USE IN PHASE2. C DO 20 K = 1, NEQ THRES = RWORK(NTHR+K) WT = MAX(THRES,ABS(Y(K))) RWORK(NWT+K) = WT TRY = ABS(YP(K)/WT) IF (TRY.GT.YPNRM) THEN YPNRM = TRY DFLT = WT .EQ. THRES END IF THRESP = RWORK(NTHRP+K) WTP = MAX(THRESP,ABS(YP(K))) RWORK(NWTP+K) = WTP TRY = ABS(YDP(K)/WTP) IF (TRY.GT.YDPNRM) THEN YDPNRM = TRY DFLTP = WTP .EQ. THRESP END IF 20 CONTINUE ABSH = RANGE ABSHP = RANGE TAU = TOL**EXPON IF (ABSH*YPNRM.GT.TAU .AND. .NOT. DFLT) ABSH = TAU/YPNRM IF (ABSHP*YDPNRM.GT.TAU .AND. .NOT. DFLTP) * ABSHP = TAU/YDPNRM H = MIN(ABSH,ABSHP) C C END OF PHASE 1. C END IF C C END OF BLOCK B. C ELSE C C BLOCK C. C CONTINUATION CALL. C CHECK DATA FOR CONSISTENCY. C TOLD = RWORK(14) HUSED = RWORK(2) IF (T.EQ.TEND) IFAIL = 1 IF (NEQ.NE.NEQOLD) IFAIL = 1 IF (DIR*(TEND-T).LT.ZERO) IFAIL = 1 IF (IFAIL.EQ.1) GO TO 620 C C CHECK FOR OPTIONAL INPUTS. C IF (RWORK(9).EQ.ONE) THEN H = RWORK(3) NSTMAX = INT(RWORK(6)+FIFTH) RWORK(9) = ZERO TOL = RWORK(12) IF (RHIGH.NE.RWORK(10)) THEN RHIGH = RWORK(10) IF (RWORK(10).EQ.ZERO) THEN EXPON = FIFTH STBRAD = FOUR TOOSML = TWENTY*RWORK(UROUND) NSTAGE = 6 SCALE = SCALEL FINPH2 = 4 NWT = NBASE + 3*NEQ NWTP = NBASE + 4*NEQ PTR(1) = 0 PTR(2) = 1 PTR(3) = 2 PTR(4) = 3 PTR(5) = 4 CALL COEF64(A,C,BPHAT,B,BP) ELSE EXPON = ELVNTH STBRAD = EIGHT TOOSML = TWOHUN*RWORK(UROUND) NSTAGE = 17 SCALE = SCALEH CALL COEFHI(A,C,BHAT,BPHAT,B,BP) FINPH2 = 14 NWT = NBASE + 12*NEQ NWTP = NBASE + 13*NEQ PTR(1) = 0 PTR(2) = 1 PTR(3) = 2 PTR(4) = 0 PTR(5) = 3 PTR(6) = 4 PTR(7) = 5 PTR(8) = 6 PTR(9) = 7 PTR(10) = 8 PTR(11) = 9 PTR(12) = 10 PTR(13) = 11 PTR(14) = 12 PTR(15) = 13 PTR(16) = 3 END IF END IF END IF C C END OF BLOCK C. C END IF C C BLOCK D. C THE START OF A NEW STEP. C H = SIGN(ABS(H),DIR) C 40 CONTINUE SUCCES = .FALSE. 60 LAST = .FALSE. C C THE CODE IS REQUIRED TO HIT TEND EXACTLY WITH A REASONABLE STEPSIZE. C CHECK THAT T+H .LE. TEND OR T+2*H .LE. TEND IN THE DIRECTION C OF INTEGRATION. A CHECK IS INCORPORATED LATER IN THE CODE TO CHECK C WHEN H IS REDUCED TOO OFTEN BY A SIGNIFICANT AMOUNT, HENCE IMPAIRING C EFFICIENCY. C IF (DIR*(T+H-TEND).GE.ZERO) THEN INEFFC = (TEND-T)/H .LT. HALF H = TEND - T LAST = .TRUE. ELSE IF (DIR*(T+TWO*H-TEND).GE.ZERO) THEN H = (TEND-T)/TWO END IF END IF C FAILED = .FALSE. C C RESTART HERE AFTER A STEP FAILURE. C 80 CONTINUE C C THE CHECK FOR LIMITING PRECISION, USING TOOSML, IS BASED ON C THE DISTANCE BETWEEN COEFFICIENTS C(I) IN THE RKN FORMULAS, C AS DESCRIBED IN C SHAMPINE,L.F.,WATTS,H.A. (1979). THE ART OF WRITING A C RUNGE-KUTTA CODE. APPL. MATH. COMPUTATION 5, 93-121. C IF (ABS(H).LT.TOOSML*MAX(ABS(T),ABS(T+H))) THEN C C STEP TOO SMALL. CONTROL RETURNS TO USER THROUGH BOTTOM OF CODE. C IFAIL = 3 GO TO 600 END IF NASTEP = NASTEP + 1 IF (NASTEP.GT.NSTMAX) THEN IFAIL = 2 GO TO 600 END IF C C END OF BLOCK D. C C BLOCK E. C FORM THE STAGES AND STORE THEM IN THE WORKING STORAGE ARRAY, RWORK. C TEMPORARY STORAGE: C RWORK(NRT1+K) CONTAINS THE SUM OF THE STAGES SUM A(I,J)*G(J) C USED IN FORMING THE INTERNAL Y VALUES. AT THE END OF C THE BLOCK IT WILL CONTAIN SUM BHAT(I) * G(I) AND C WILL BE USED IN ERROR ESTIMATION FOR Y. C RWORK(NRT2+K) WILL CONTAIN SUM BPHAT(I) * G(I) AND WILL C BE USED IN ESTIMATING THE ERROR IN YP. C AT THE END OF THE BLOCK, THE NEW APPROXIMATION TO Y(K) WILL BE C IN RWORK(NRT3+K) AND THE APPROXIMATION TO YP(K) WILL BE IN C RWORK(NRT4+K). C C NOTE. C MANY OF THE FOLLOWING LOOPS OVER K=1,NEQ HAVE CONSTANT ARRAY VALUES C INSIDE. THE CODE IS WRITTEN WITH CLARITY IN MIND. ANY OPTIMIZING C COMPILER WILL IDENTIFY THESE OCCURRENCES AND MOVE THESE CONSTANT C VALUES OUTSIDE THE LOOPS. TWO LOOPS CONTAIN CHECKS FOR ZERO C MULTIPLIERS SO AS TO PREVENT NEEDLESS COMPUTATION. C HSQ = H*H DO 100 K = 1, NEQ RWORK(NRT2+K) = BPHAT(1)*YDP(K) 100 CONTINUE DO 260 I = 2, NSTAGE C DO 120 K = 1, NEQ RWORK(NRT1+K) = A(I,1)*YDP(K) 120 CONTINUE C DO 160 J = 2, I - 1 JSTAGE = NBASE + PTR(J-1)*NEQ IF (A(I,J).NE.ZERO) THEN DO 140 K = 1, NEQ RWORK(NRT1+K) = RWORK(NRT1+K) + A(I,J)*RWORK(JSTAGE+K) 140 CONTINUE END IF 160 CONTINUE C DO 180 K = 1, NEQ RWORK(NRT3+K) = Y(K) + C(I)*H*YP(K) + HSQ*RWORK(NRT1+K) 180 CONTINUE C JSTAGE = NBASE + PTR(I-1)*NEQ CALL F(NEQ,T+C(I)*H,RWORK(NRT3+1),RWORK(JSTAGE+1)) C IF (BPHAT(I).NE.ZERO) THEN DO 200 K = 1, NEQ RWORK(NRT2+K) = RWORK(NRT2+K) + BPHAT(I)*RWORK(JSTAGE+K) 200 CONTINUE END IF C C CONTROL USUALLY GOES FROM HERE TO THE START OF THE LOOP (I=2,NSTAGE), C UNLESS IT IS THE FIRST STEP AND THE CODE IS COMPUTING H0, C IN WHICH CASE THE NEXT BLOCK OF CODE IS EXECUTED. C IF (PHASE2) THEN YPNRM = ZERO YDPNRM = ZERO DEL = ZERO FDEL = ZERO YPDEL = ZERO YDPDEL = ZERO CVEC CVEC THE FOLLOWING LOOP WAS NOT VECTORIZABLE ON A CDC CYBER 205. CVEC DO 220 K = 1, NEQ YINTK = RWORK(NRT3+K) YPINTK = YP(K) + H*RWORK(NRT2+K) WT = MAX(RWORK(NWT+K),ABS(YINTK)) WTP = MAX(RWORK(NWTP+K),ABS(YPINTK)) RWORK(NWT+K) = WT RWORK(NWTP+K) = WTP YPNRM = MAX(YPNRM,ABS(YINTK)/WT,ABS(Y(K))/WT) YDPNRM = MAX(YDPNRM,ABS(YPINTK)/WTP,ABS(YP(K))/WTP) DEL = MAX(DEL,ABS(YINTK-Y(K))/WT) FDEL = MAX(FDEL,ABS(YPINTK-YP(K))/WT) YPDEL = MAX(YPDEL,ABS(YPINTK-YP(K))/WTP) YDPDEL = MAX(YDPDEL,ABS(RWORK(JSTAGE+K)-YDP(K))/WTP) 220 CONTINUE DEL = MAX(DEL,ABS(C(I)*H)/RANGE) SCL = ONE SCLP = ONE IF (DEL.GT.TEN*RWORK(UROUND)*YPNRM) THEN IF (ABS(H)*FDEL.GT.STBRAD*DEL) * SCL = STBRAD/R*MAX(DEL/(ABS(H)*FDEL),RMN3) END IF IF (YPDEL.GT.TEN*RWORK(UROUND)*YDPNRM) THEN IF (ABS(H)*YDPDEL.GT.STBRAD*YPDEL) * SCLP = STBRAD/R*MAX(YPDEL/(ABS(H)*YDPDEL),RMN3) END IF SCALE1 = MIN(SCL,SCLP) IF (SCALE1.LT.ONE) THEN NATT = NATT + 1 H = SCALE1*ABS(H) H = DIR*H LAST = .FALSE. C C RESET THE WEIGHT VECTORS C DO 240 K = 1, NEQ RWORK(NWT+K) = MAX(ABS(Y(K)),RWORK(NTHR+K)) RWORK(NWTP+K) = MAX(ABS(YP(K)),RWORK(NTHRP+K)) 240 CONTINUE GO TO 80 END IF C C TO SAVE ON STORAGE THE WEIGHT VECTORS RWORK(NWT+...) AND C RWORK(NWTP+...) ARE IN LOCATIONS WHERE THE STAGES ARE USUALLY KEPT C AND HENCE PHASE2 OPERATES UPTO THE FOURTH STAGE FOR THE 6(4) PAIR AND C THE 14TH STAGE FOR THE 12(10) PAIR. C PHASE2 = I .LT. FINPH2 END IF 260 CONTINUE C DO 280 K = 1, NEQ RWORK(NRT4+K) = YP(K) + H*RWORK(NRT2+K) 280 CONTINUE C C THE 12(10) PAIR DOES NOT USE FSAL (FIRST STAGE ON STEP N+1 IS SAME C AS LAST STAGE ON STEP N), HENCE MUST ADD STAGES UP TO OBTAIN THE C APPROXIMATIONS. UTILIZE FACT THAT C BHAT(2)=BHAT(3)=BHAT(4)=BHAT(5)=BHAT(6)=BHAT(16)=BHAT(17)=0. C IF (RHIGH.EQ.ONE) THEN DO 300 K = 1, NEQ RWORK(NRT1+K) = BHAT(1)*YDP(K) 300 CONTINUE DO 340 J = 7, 15 JSTAGE = NBASE + PTR(J-1)*NEQ DO 320 K = 1, NEQ RWORK(NRT1+K) = RWORK(NRT1+K) + BHAT(J)*RWORK(JSTAGE+K) 320 CONTINUE 340 CONTINUE DO 360 K = 1, NEQ RWORK(NRT3+K) = Y(K) + H*YP(K) + HSQ*RWORK(NRT1+K) 360 CONTINUE END IF C C END OF BLOCK E. C C C BLOCK F. C ERROR ESTIMATION. C USE THE TEMPORARY STORAGE VECTORS RWORK(NRT1+K) AND RWORK(NRT2+K) C TO COMPUTE YHAT - Y AND YPHAT - YP RESPECTIVELY. C SEPARATE SECTIONS ARE USED FOR THE 12(10) AND 6(4) PAIRS C IN ORDER TO TAKE ADVANTAGE OF THE FACT THAT FOR THE 6(4) PAIR C CASE, B(4) = B(5) = B(6) = ZERO, WHILST FOR THE 12(10) PAIR C B(2) = B(3) = B(4) = B(5) = B(6) = B(17) = ZERO AND SIMILARLY FOR C BP. C DO 380 K = 1, NEQ RWORK(NRT1+K) = RWORK(NRT1+K) - B(1)*YDP(K) RWORK(NRT2+K) = RWORK(NRT2+K) - BP(1)*YDP(K) 380 CONTINUE C IF (RHIGH.EQ.ZERO) THEN C C ERROR ESTIMATES FOR THE 6(4) PAIR. C DO 420 I = 2, 3 JSTAGE = NBASE + PTR(I-1)*NEQ DO 400 K = 1, NEQ RWORK(NRT1+K) = RWORK(NRT1+K) - B(I)*RWORK(JSTAGE+K) 400 CONTINUE 420 CONTINUE DO 460 I = 2, NSTAGE JSTAGE = NBASE + PTR(I-1)*NEQ DO 440 K = 1, NEQ RWORK(NRT2+K) = RWORK(NRT2+K) - BP(I)*RWORK(JSTAGE+K) 440 CONTINUE 460 CONTINUE C ELSE C C ERROR ESTIMATES FOR THE 12(10) PAIR. C DO 500 I = 7, 16 JSTAGE = NBASE + PTR(I-1)*NEQ DO 480 K = 1, NEQ RWORK(NRT1+K) = RWORK(NRT1+K) - B(I)*RWORK(JSTAGE+K) RWORK(NRT2+K) = RWORK(NRT2+K) - BP(I)*RWORK(JSTAGE+K) 480 CONTINUE 500 CONTINUE END IF C C FORM ERR, THE MAXIMUM OF THE WEIGHTED MAX NORMS OF THE ESTIMATED C LOCAL ERRORS IN Y AND YP. C ERR = ZERO ERRP = ZERO DO 520 K = 1, NEQ SK = HSQ*RWORK(NRT1+K) SPK = H*RWORK(NRT2+K) WT = HALF*ABS(Y(K)) + QUAR*ABS(RWORK(NRT3+K)) + * QUAR*ABS(RWORK(NRT3+K)-SK) WT = MAX(WT,RWORK(NTHR+K)) WTP = HALF*ABS(YP(K)) + QUAR*ABS(RWORK(NRT4+K)) + * QUAR*ABS(RWORK(NRT4+K)-SPK) WTP = MAX(WTP,RWORK(NTHRP+K)) ERR = MAX(ERR,ABS(SK)/WT) ERRP = MAX(ERRP,ABS(SPK)/WTP) 520 CONTINUE ERR = MAX(ERR,ERRP) C C END OF BLOCK F. C IF (ERR.GT.TOL) THEN C C BLOCK G. C FAILED STEP. C LAST = .FALSE. IF (FIRST .AND. .NOT. HGIVEN) THEN C C PHASE 3 FOR CODE CHOICE OF H0. C NATT = NATT + 1 IF (SUCCES) THEN C C THE CODE HAS DISCARDED AN INITIAL STEP IN FAVOUR OF A MUCH LARGER C PREDICTED STEP BUT THIS LARGER STEP HAS FAILED. THEREFORE CARRY C OUT OPTIMAL REDUCTION. C OPTRED = .TRUE. ALPHA = SCALE*(TOL/ERR)**EXPON ALPHA = MAX(RMN2,ALPHA) ELSE C C NO SUCCESSFUL STEP YET. REDUCE H0 TO H0/R AND RESET WEIGHT VECTORS. C ALPHA = RMN1 PHASE2 = .TRUE. DO 540 K = 1, NEQ RWORK(NWT+K) = MAX(ABS(Y(K)),RWORK(NTHR+K)) RWORK(NWTP+K) = MAX(ABS(YP(K)),RWORK(NTHRP+K)) 540 CONTINUE END IF ELSE C C NOT THE FIRST STEP (OR H0 WAS SUPPLIED BY THE USER), SO USE THE C NORMAL STEP REDUCTION ALGORITHM. C NFLSTP = NFLSTP + 1 IF (FAILED) THEN ALPHA = HALF ELSE FAILED = .TRUE. ALPHA = SCALE*(TOL/ERR)**EXPON ALPHA = MAX(RMN1,ALPHA) END IF END IF H = ALPHA*H C C TRY THE STEP AGAIN. C GO TO 80 C C END OF BLOCK G. C ELSE C C BLOCK H. C SUCCESSFUL STEP. C SUCCES = .TRUE. BETA = ((ERR/TOL)**EXPON)/SCALE IF (FIRST) THEN C C PHASE 3 FOR CODE CHOICE OF H0. C ALPHA = ONE/MAX(RMN3,BETA) ELSE C C NORMAL STEPSIZE PREDICTION. C ALPHA = ONE/MAX(RMN1,BETA) IF (FAILED) ALPHA = MIN(ONE,ALPHA) END IF HUSED = H H = ALPHA*H C IF (FIRST .AND. .NOT. LAST .AND. .NOT. OPTRED .AND. ALPHA.GT.R) * THEN C C PHASE 3 FOR CODE CHOICE OF H0. THE CODE HAS ATTEMPTED THE FIRST STEP, C IS NOT AT THE END OF THE INTEGRATION RANGE, HAS NOT ENCOUNTERED A STEP C FAILURE IN PHASE3, AND THE PREDICTED INCREASE IS LARGER THAN THE C MAXIMUM PERMITTED ON A NORMAL STEP. RETAKE THE STEP. C HUSED = ZERO NATT = NATT + 1 GO TO 60 END IF C IF (FIRST) RWORK(1) = HUSED FIRST = .FALSE. TOLD = T T = T + HUSED EFFCT1 = EFFCT1 + 1 IF (LAST) THEN T = TEND IF (INEFFC) EFFCT2 = EFFCT2 + 1 C C TEST FOR INEFFICIENT USE OF FORMULAS FOR OUTPUT. C IF (EFFCT1.GT.HUNDRD) THEN WASTE = DBLE(EFFCT2)/DBLE(EFFCT1) .GT. TENPC IF (WASTE) THEN IFAIL = 4 EFFCT1 = 0 EFFCT2 = 0 END IF END IF END IF ISTP = ISTP + 1 C C ON A SUCCESSFUL STEP, COPY THE NEW APPROXIMATIONS INTO Y AND YP. C SAVE THE OLD VALUES OF Y, YP AND YDP IN RWORK, C ALONG WITH THE INTERMEDIATE STAGES FOR USE IN CASE DENSE C OUTPUT IS REQUIRED. C DO 560 K = 1, NEQ RWORK(NRT1+K) = Y(K) Y(K) = RWORK(NRT3+K) RWORK(NRT2+K) = YP(K) YP(K) = RWORK(NRT4+K) RWORK(NRT3+K) = YDP(K) 560 CONTINUE C C SINCE THE 12(10) PAIR DOES NOT USE FSAL (SEE ABOVE), A FUNCTION C EVALUATION IS DONE HERE WHICH WILL BE USED ON THE NEXT STEP. FOR THE C 6(4) PAIR SET YDP FROM THE LAST STAGE. C IF (RHIGH.EQ.ONE) THEN CALL F(NEQ,T,Y,YDP) ELSE JSTAGE = NBASE + PTR(5)*NEQ DO 580 K = 1, NEQ YDP(K) = RWORK(JSTAGE+K) 580 CONTINUE END IF C C IF THE CODE HAS NOT REACHED THE END OF THE INTEGRATION RANGE OR IS C OPERATING IN INTERVAL MODE, ATTEMPT THE NEXT STEP. C IF ( .NOT. (LAST .OR. (RWORK(8).EQ.ONE))) GO TO 40 C C END OF BLOCK H. C END IF C C BLOCK I. C SET DIAGNOSTIC INFORMATION AND EXIT. C 600 RWORK(2) = HUSED RWORK(3) = H RWORK(4) = DBLE(ISTP) RWORK(5) = DBLE(NFLSTP) RWORK(11) = DBLE(NATT) RWORK(14) = TOLD 620 CONTINUE C C TEST FOR SUCCESSIVE FAILURES AT THE SAME VALUE OF T. C IF (IFAIL.GT.0) THEN IF ( .NOT. LOOP) THEN LOOP = .TRUE. TCHECK = T ELSE IF (TCHECK.EQ.T) THEN WRITE (UNIT=IDEV,FMT=99999) T STOP ELSE LOOP = .FALSE. END IF END IF ELSE LOOP = .FALSE. END IF C C END OF BLOCK I. C RETURN C C 99999 FORMAT (/' RKNINT STOPPED DUE TO TWO SUCCESSIVE FAILURES',' AT', * ' THE SAME VALUE OF T (= ',1P,D12.5,')',/) 99998 FORMAT (/' RKNINT STOPPED DUE TO NONZERO INPUT VALUE OF IFAIL. ', * 'IFAIL = ',I3) END SUBROUTINE RKNDIA(NEQ,HNEXT,HUSED,HSTART,NSUCC,NFAIL,NATT,THRES, * THRESP,RWORK) C C THIS IS THE DIAGNOSTIC ROUTINE ASSOCIATED WITH THE NYSTROM C SOLVER RKNINT. THE VALUES OF THE PARAMETERS IN THE SUBROUTINE C CALL ARE TAKEN OUT OF RWORK. C C SEE COMMENT IN RKNSET ABOUT STORAGE REDUCTION FOR THRES AND THRESP. C C .. Parameters .. DOUBLE PRECISION XINC INTEGER NOVHD PARAMETER (XINC=0.2D0,NOVHD=14) C .. C .. Scalar Arguments .. DOUBLE PRECISION HNEXT, HSTART, HUSED INTEGER NATT, NEQ, NFAIL, NSUCC C .. C .. Array Arguments .. DOUBLE PRECISION RWORK(*), THRES(NEQ), THRESP(NEQ) C .. C .. Local Scalars .. INTEGER K, NTHR, NTHRP C .. C .. Intrinsic Functions .. INTRINSIC INT C .. C .. Executable Statements .. C NTHR = NOVHD + 4*NEQ NTHRP = NOVHD + 5*NEQ C HSTART = RWORK(1) HUSED = RWORK(2) HNEXT = RWORK(3) NSUCC = INT(RWORK(4)+XINC) NFAIL = INT(RWORK(5)+XINC) NATT = INT(RWORK(11)+XINC) DO 20 K = 1, NEQ THRES(K) = RWORK(NTHR+K) THRESP(K) = RWORK(NTHRP+K) 20 CONTINUE RETURN END SUBROUTINE EVAL(NEQ,T,Y,YP,NWANT,TWANT,YWANT,YPWANT,RWORK,IFAIL) C C THIS ROUTINE USES THE THIRD FORMULA OF THE RKN6(4)6 TRIPLE C TO RETURN THE VALUES OF Y AND YP AT THE POINT TWANT, WHICH C SHOULD BE BETWEEN T-H AND T. C C .. Parameters .. INTEGER NOVHD DOUBLE PRECISION ONE PARAMETER (NOVHD=14,ONE=1.0D0) C .. C .. Scalar Arguments .. DOUBLE PRECISION T, TWANT INTEGER IFAIL, NEQ, NWANT C .. C .. Array Arguments .. DOUBLE PRECISION RWORK(NOVHD+14*NEQ), Y(NEQ), YP(NEQ), * YPWANT(NWANT), YWANT(NWANT) C .. C .. Local Scalars .. DOUBLE PRECISION H, HSIG, SIGMA, SUM, SUMP, TOLD INTEGER I, JSTAGE, K, NBASE, NRT1, NRT2, NRT3, NSTAGE LOGICAL INTERP C .. C .. Local Arrays .. DOUBLE PRECISION BPSTAR(6), BSTAR(6) C .. C .. External Subroutines .. EXTERNAL CFEVAL2 C .. C .. Intrinsic Functions .. INTRINSIC SIGN C .. C .. Executable Statements .. C IF (RWORK(10).EQ.ONE) THEN IFAIL = 3 GO TO 100 END IF NSTAGE = 6 IF (NWANT.GT.NEQ .OR. NWANT.LT.1) THEN IFAIL = 2 GO TO 100 END IF C C RWORK(NRT1+...) CONTAINS THE PREVIOUS VALUE OF Y. C RWORK(NRT2+...) CONTAINS THE PREVIOUS VALUE OF YP. C RWORK(NRT3+...) CONTAINS THE PREVIOUS VALUE OF YDP. C NRT1 = NOVHD NRT2 = NOVHD + NEQ NRT3 = NOVHD + 2*NEQ NBASE = NOVHD + 6*NEQ C C CHECK FOR TRIVIAL CASES C IFAIL = 0 IF (TWANT.EQ.T) THEN DO 20 K = 1, NWANT YWANT(K) = Y(K) YPWANT(K) = YP(K) 20 CONTINUE GO TO 100 END IF H = RWORK(2) TOLD = RWORK(14) IF (TWANT.EQ.TOLD) THEN DO 40 K = 1, NWANT YWANT(K) = RWORK(NRT1+K) YPWANT(K) = RWORK(NRT2+K) 40 CONTINUE GO TO 100 END IF C C CHECK FOR EXTRAPOLATION C IF (SIGN(ONE,H).EQ.ONE) THEN INTERP = (TOLD.LT.TWANT) .AND. (TWANT.LT.T) ELSE INTERP = (TOLD.GT.TWANT) .AND. (TWANT.GT.T) END IF IF ( .NOT. INTERP) IFAIL = 1 C HSIG = TWANT - TOLD SIGMA = HSIG/H CALL CFEVAL2(BSTAR,BPSTAR,SIGMA) C C COMPUTE YWANT AND YPWANT. C THE FACT THAT BSTAR(2) = 0 = BPSTAR(2) IS USED. C DO 80 K = 1, NWANT SUM = BSTAR(1)*RWORK(NRT3+K) SUMP = BPSTAR(1)*RWORK(NRT3+K) DO 60 I = 3, NSTAGE JSTAGE = NBASE + (I-2)*NEQ SUM = SUM + BSTAR(I)*RWORK(JSTAGE+K) SUMP = SUMP + BPSTAR(I)*RWORK(JSTAGE+K) 60 CONTINUE YWANT(K) = RWORK(NRT1+K) + HSIG*(RWORK(NRT2+K)+HSIG*SUM) YPWANT(K) = RWORK(NRT2+K) + HSIG*SUMP 80 CONTINUE C 100 RETURN END SUBROUTINE COEF64(A,C,BPHAT,B,BP) C C THIS ROUTINE RETURNS THE COEFFICIENTS FOR THE C DORMAND/PRINCE RKN6(4)6FD PAIR USED IN THE CODE C TO COMPUTE Y AND YP. THE COEFFICIENTS FOR THE DENSE C OUTPUT FORMULA ARE RETURNED BY THE ROUTINE CFEVAL2 TO C THE INTERPOLATION SUBROUTINE. C C .. Array Arguments .. DOUBLE PRECISION A(17,17), B(16), BP(16), BPHAT(17), C(17) C .. C .. Executable Statements .. C DOUBLE PRECISION R C INTRINSIC SQRT C C R = SQRT(8581.0) C C C(1) = 0. C C(2) = (209.0-R)/900.0 C C(3) = (209.0-R)/450.0 C C(4) = (209.0+R)/450.0 C C(5) = 9./10. C C(6) = 1. C C A(2,1) = (26131.0-209.0*R)/810000.0 C A(3,1) = (26131.0-209.0*R)/607500.0 C A(3,2) = (26131.0-209.0*R)/303750.0 C A(4,1) = (980403512254.0+7781688431.0*R)/11694469921875.0 C A(4,2) = -(1262884486208.0+15385481287.0*R)/11694469921875.0 C A(4,3) = (7166233891441.0+78694563299.0*R)/46777879687500.0 C A(5,1) = -9.0*(329260.0+3181.0*R)/27040000.0 C A(5,2) = 27.0*(35129.0+3331.0*R)/13520000.0 C A(5,3) = -27.0*(554358343.0+31040327.0*R)/464060480000.0 C A(5,4) = 153.0*(8555257.0-67973.0*R)/2745920000.0 C A(6,1) = 329.0/4212.0 C A(6,2) = 0. C A(6,3) = (84119543.0+366727.0*R)/409622616.0 C A(6,4) = (84119543.0-366727.0*R)/409622616.0 C A(6,5) = 200.0/17901.0 C C BPHAT(1) = 329.0/4212.0 C BPHAT(2) = 0. C BPHAT(3) = (389225579.0+96856.0*R)/1024056540.0 C BPHAT(4) = (389225579.0-96856.0*R)/1024056540.0 C BPHAT(5) = 2000.0/17901.0 C BPHAT(6) = 1./20. C C B(1) = (2701.0+23.0*R)/4563.0 C B(2) = -(9829.0+131.0*R)/9126.0 C B(3) = 5.0*(1798.0+17.0*R)/9126.0 C C BP(1) = 115.0/2106.0 C BP(2) = 0. C BP(3) = (84119543.0+366727.0*R)/256014135.0 C BP(4) = (84119543.0-366727.0*R)/256014135.0 C BP(5) = 6950.0/17901.0 C BP(6) = -1./10. C C COEFFICIENTS FOR RKN646FD COMPUTED TO 30D BY MACSYMA C C(1) = 0.0D0 C(2) = 1.29295903136704415288990053209D-1 C(3) = 2.58591806273408830577980106418D-1 C(4) = 6.70297082615480058310908782471D-1 C(5) = 9.0D-1 C(6) = 1.0D0 C A(2,1) = 8.35871528396802532822102346743D-3 A(3,1) = 1.11449537119573671042946979566D-2 A(3,2) = 2.22899074239147342085893959132D-2 A(4,1) = 1.45474742801091785895935232316D-1 A(4,2) = -2.29860640522647473120261456297D-1 A(4,3) = 3.09034987202967536528726080729D-1 A(5,1) = -2.07668262950789954335146205252D-1 A(5,2) = 6.86366784292514312273571851042D-1 A(5,3) = -1.99549277872349252201326358888D-1 A(5,4) = 1.25850756530624894262900713098D-1 A(6,1) = 7.81101614434947768281101614435D-2 A(6,2) = 0.0D0 A(6,3) = 2.882917411897667776841772093D-1 A(6,4) = 1.22425537174570410182422422005D-1 A(6,5) = 1.1172560192168035305290207251D-2 C BPHAT(1) = 7.81101614434947768281101614435D-2 BPHAT(2) = 0.0D0 BPHAT(3) = 3.88843478705982602715952208523D-1 BPHAT(4) = 3.71320757928842267403035557523D-1 BPHAT(5) = 1.1172560192168035305290207251D-1 BPHAT(6) = 5.0D-2 C B(1) = 1.05885926037041827822001005886D0 B(2) = -2.40675137192445205319176777632D0 B(3) = 1.84789211155403377497175771746D0 B(4) = 0.0D0 B(5) = 0.0D0 B(6) = 0.0D0 C BP(1) = 5.46058879392212725546058879392D-2 BP(2) = 0.0D0 BP(3) = 4.6126678590362684429468353488D-1 BP(4) = 1.95880859479312656291875875208D-1 BP(5) = 3.88246466677839226858834701972D-1 BP(6) = -1.0D-1 C RETURN END SUBROUTINE CFEVAL2(B,BP,S) C C THIS ROUTINE RETURNS THE COEFFICIENTS FOR Y AND YP FROM THE C DORMAND/PRINCE DENSE OUTPUT PAIR OF RKN6(4)6FD. C C .. Scalar Arguments .. DOUBLE PRECISION S C .. C .. Array Arguments .. DOUBLE PRECISION B(6), BP(6) C .. C .. Local Scalars .. DOUBLE PRECISION U, V, W, X C .. C .. Executable Statements .. C C R = SQRT(8581.0D0) C B(1) = (((900.D0*S-3819.D0)*S+6386.D0)*S-5244.D0)*S + 2106.D0 B(1) = B(1)/4212.D0 BP(1) = (((5400.D0*S-19095.D0)*S+25544.D0)*S-15732.D0)*S + 4212.D0 BP(1) = BP(1)/4212.D0 B(2) = 0.D0 BP(2) = 0.D0 C C U = 5860823.0D0 - 152228.0D0*R C V = 4929647204.0D0 - 156109769.0D0*R C W = 22190560391.0D0 - 1109665151.0D0*R C X = 81356461.0D0 + 25954829.0D0*R C Y = 22529243880.0D0 C Z = 11264621940.0D0 C B(3) = (((1800.*U*S-6.*V)*S+W)*S+18.*X)*S/Y C BP(3) = (((5400.*U*S-15.*V)*S+2.*W)*S+27.*X)*S/Z C U = -.365774278775981018074400638167D-3 V = -.423066852735158392593161679389D0 W = -.357765287229492094385275612595D1 X = .110329844381694622139549133756D0 B(3) = (((1800.D0*U*S-6.D0*V)*S+W)*S+18.D0*X)*S U = -.731548557551962036148801276334D-3 V = -.846133705470316785186323358781D0 W = -.715530574458984188770551225188D1 X = .220659688763389244279098267512D0 BP(3) = (((5400.D0*U*S-15.D0*V)*S+2.D0*W)*S+27.D0*X)*S C C U = 5860823.0D0 + 152228.0D0*R C V = 4929647204.0D0 + 156109769.0D0*R C W = 22190560391.0D0 + 1109665151.0D0*R C X = 81356461.0D0 - 25954829.0D0*R C B(4) = (((1800.*U*S-6.*V)*S+W)*S+18.*X)*S/Y C BP(4) = (((5400.*U*S-15.*V)*S+2.*W)*S+27.*X)*S/Z C U = .886060093179442480359392334923D-3 V = .860688925651348955842548458187D0 W = .554758675105232911302563118117D1 X = -.103107545982925635135268011828D0 B(4) = (((1800.D0*U*S-6.D0*V)*S+W)*S+18.D0*X)*S U = .177212018635888496071878466985D-2 V = .172137785130269791168509691638D1 W = .110951735021046582260512623623D2 X = -.206215091965851270270536023656D0 BP(4) = (((5400.D0*U*S-15.D0*V)*S+2.D0*W)*S+27.D0*X)*S C B(5) = ((225.D0*S-651.D0)*S+620.D0)*S - 195.D0 B(5) = -200.D0*S*B(5)/17901.D0 BP(5) = ((270.D0*S-651.D0)*S+496.D0)*S - 117.D0 BP(5) = -1000.D0*S*BP(5)/17901.D0 B(6) = S*(S-1.D0)*((300.D0*S-523.D0)*S+234.D0)/220.D0 BP(6) = S*(((1800.D0*S-4115.D0)*S+3028.D0)*S-702.D0)/220.D0 C RETURN END SUBROUTINE COEFHI(A,C,BCAP,BDCAP,B,BD) C C THIS ROUTINE RETURNS THE COEFFICIENTS FOR THE RUNGE-KUTTA C NYSTROM 12(10) PAIR OF DORMAND AND PRINCE. THE COEFFICIENTS C WERE SUPPLIED BY D/P. C COEFFICIENTS OF RKN12(10)17M TO 30 DECIMAL PLACES FROM MACSYMA. C C .. Array Arguments .. DOUBLE PRECISION A(17,17), B(16), BCAP(15), BD(16), BDCAP(17), * C(17) C .. C .. Executable Statements .. C C(1) = 0.0D0 C(2) = 2.0D-2 C(3) = 4.0D-2 C(4) = 1.0D-1 C(5) = 1.33333333333333333333333333333D-1 C(6) = 1.6D-1 C(7) = 5.0D-2 C(8) = 2.0D-1 C(9) = 2.5D-1 C(10) = 3.33333333333333333333333333333D-1 C(11) = 5.0D-1 C(12) = 5.55555555555555555555555555556D-1 C(13) = 7.5D-1 C(14) = 8.57142857142857142857142857143D-1 C(15) = 9.45216222272014340129957427739D-1 C(16) = 1.0D0 C(17) = 1.0D0 C A(2,1) = 2.0D-4 C A(3,1) = 2.66666666666666666666666666667D-4 A(3,2) = 5.33333333333333333333333333333D-4 C A(4,1) = 2.91666666666666666666666666667D-3 A(4,2) = -4.16666666666666666666666666667D-3 A(4,3) = 6.25D-3 C A(5,1) = 1.64609053497942386831275720165D-3 A(5,2) = 0.0D0 A(5,3) = 5.48696844993141289437585733882D-3 A(5,4) = 1.75582990397805212620027434842D-3 C A(6,1) = 1.9456D-3 A(6,2) = 0.0D0 A(6,3) = 7.15174603174603174603174603175D-3 A(6,4) = 2.91271111111111111111111111111D-3 A(6,5) = 7.89942857142857142857142857143D-4 C A(7,1) = 5.6640625D-4 A(7,2) = 0.0D0 A(7,3) = 8.80973048941798941798941798942D-4 A(7,4) = -4.36921296296296296296296296296D-4 A(7,5) = 3.39006696428571428571428571429D-4 A(7,6) = -9.94646990740740740740740740741D-5 C A(8,1) = 3.08333333333333333333333333333D-3 A(8,2) = 0.0D0 A(8,3) = 0.0D0 A(8,4) = 1.77777777777777777777777777778D-3 A(8,5) = 2.7D-3 A(8,6) = 1.57828282828282828282828282828D-3 A(8,7) = 1.08606060606060606060606060606D-2 C A(9,1) = 3.65183937480112971375119150338D-3 A(9,2) = 0.0D0 A(9,3) = 3.96517171407234306617557289807D-3 A(9,4) = 3.19725826293062822350093426091D-3 A(9,5) = 8.22146730685543536968701883401D-3 A(9,6) = -1.31309269595723798362013884863D-3 A(9,7) = 9.77158696806486781562609494147D-3 A(9,8) = 3.75576906923283379487932641079D-3 C A(10,1) = 3.70724106871850081019565530521D-3 A(10,2) = 0.0D0 A(10,3) = 5.08204585455528598076108163479D-3 A(10,4) = 1.17470800217541204473569104943D-3 A(10,5) = -2.11476299151269914996229766362D-2 A(10,6) = 6.01046369810788081222573525136D-2 A(10,7) = 2.01057347685061881846748708777D-2 A(10,8) = -2.83507501229335808430366774368D-2 A(10,9) = 1.48795689185819327555905582479D-2 C A(11,1) = 3.51253765607334415311308293052D-2 A(11,2) = 0.0D0 A(11,3) = -8.61574919513847910340576078545D-3 A(11,4) = -5.79144805100791652167632252471D-3 A(11,5) = 1.94555482378261584239438810411D0 A(11,6) = -3.43512386745651359636787167574D0 A(11,7) = -1.09307011074752217583892572001D-1 A(11,8) = 2.3496383118995166394320161088D0 A(11,9) = -7.56009408687022978027190729778D-1 A(11,10) = 1.09528972221569264246502018618D-1 C A(12,1) = 2.05277925374824966509720571672D-2 A(12,2) = 0.0D0 A(12,3) = -7.28644676448017991778247943149D-3 A(12,4) = -2.11535560796184024069259562549D-3 A(12,5) = 9.27580796872352224256768033235D-1 A(12,6) = -1.65228248442573667907302673325D0 A(12,7) = -2.10795630056865698191914366913D-2 A(12,8) = 1.20653643262078715447708832536D0 A(12,9) = -4.13714477001066141324662463645D-1 A(12,10) = 9.07987398280965375956795739516D-2 A(12,11) = 5.35555260053398504916870658215D-3 C A(13,1) = -1.43240788755455150458921091632D-1 A(13,2) = 0.0D0 A(13,3) = 1.25287037730918172778464480231D-2 A(13,4) = 6.82601916396982712868112411737D-3 A(13,5) = -4.79955539557438726550216254291D0 A(13,6) = 5.69862504395194143379169794156D0 A(13,7) = 7.55343036952364522249444028716D-1 A(13,8) = -1.27554878582810837175400796542D-1 A(13,9) = -1.96059260511173843289133255423D0 A(13,10) = 9.18560905663526240976234285341D-1 A(13,11) = -2.38800855052844310534827013402D-1 A(13,12) = 1.59110813572342155138740170963D-1 C A(14,1) = 8.04501920552048948697230778134D-1 A(14,2) = 0.0D0 A(14,3) = -1.66585270670112451778516268261D-2 A(14,4) = -2.1415834042629734811731437191D-2 A(14,5) = 1.68272359289624658702009353564D1 A(14,6) = -1.11728353571760979267882984241D1 A(14,7) = -3.37715929722632374148856475521D0 A(14,8) = -1.52433266553608456461817682939D1 A(14,9) = 1.71798357382154165620247684026D1 A(14,10) = -5.43771923982399464535413738556D0 A(14,11) = 1.38786716183646557551256778839D0 A(14,12) = -5.92582773265281165347677029181D-1 A(14,13) = 2.96038731712973527961592794552D-2 C A(15,1) = -9.13296766697358082096250482648D-1 A(15,2) = 0.0D0 A(15,3) = 2.41127257578051783924489946102D-3 A(15,4) = 1.76581226938617419820698839226D-2 A(15,5) = -1.48516497797203838246128557088D1 A(15,6) = 2.15897086700457560030782161561D0 A(15,7) = 3.99791558311787990115282754337D0 A(15,8) = 2.84341518002322318984542514988D1 A(15,9) = -2.52593643549415984378843352235D1 A(15,10) = 7.7338785423622373655340014114D0 A(15,11) = -1.8913028948478674610382580129D0 A(15,12) = 1.00148450702247178036685959248D0 A(15,13) = 4.64119959910905190510518247052D-3 A(15,14) = 1.12187550221489570339750499063D-2 C A(16,1) = -2.75196297205593938206065227039D-1 A(16,2) = 0.0D0 A(16,3) = 3.66118887791549201342293285553D-2 A(16,4) = 9.7895196882315626246509967162D-3 A(16,5) = -1.2293062345886210304214726509D1 A(16,6) = 1.42072264539379026942929665966D1 A(16,7) = 1.58664769067895368322481964272D0 A(16,8) = 2.45777353275959454390324346975D0 A(16,9) = -8.93519369440327190552259086374D0 A(16,10) = 4.37367273161340694839327077512D0 A(16,11) = -1.83471817654494916304344410264D0 A(16,12) = 1.15920852890614912078083198373D0 A(16,13) = -1.72902531653839221518003422953D-2 A(16,14) = 1.93259779044607666727649875324D-2 A(16,15) = 5.20444293755499311184926401526D-3 C A(17,1) = 1.30763918474040575879994562983D0 A(17,2) = 0.0D0 A(17,3) = 1.73641091897458418670879991296D-2 A(17,4) = -1.8544456454265795024362115588D-2 A(17,5) = 1.48115220328677268968478356223D1 A(17,6) = 9.38317630848247090787922177126D0 A(17,7) = -5.2284261999445422541474024553D0 A(17,8) = -4.89512805258476508040093482743D1 A(17,9) = 3.82970960343379225625583875836D1 A(17,10) = -1.05873813369759797091619037505D1 A(17,11) = 2.43323043762262763585119618787D0 A(17,12) = -1.04534060425754442848652456513D0 A(17,13) = 7.17732095086725945198184857508D-2 A(17,14) = 2.16221097080827826905505320027D-3 A(17,15) = 7.00959575960251423699282781988D-3 A(17,16) = 0.0D0 C C BCAP(1) = 1.21278685171854149768890395495D-2 BCAP(2) = 0.0D0 BCAP(3) = 0.0D0 BCAP(4) = 0.0D0 BCAP(5) = 0.0D0 BCAP(6) = 0.0D0 BCAP(7) = 8.62974625156887444363792274411D-2 BCAP(8) = 2.52546958118714719432343449316D-1 BCAP(9) = -1.97418679932682303358307954886D-1 BCAP(10) = 2.03186919078972590809261561009D-1 BCAP(11) = -2.07758080777149166121933554691D-2 BCAP(12) = 1.09678048745020136250111237823D-1 BCAP(13) = 3.80651325264665057344878719105D-2 BCAP(14) = 1.16340688043242296440927709215D-2 BCAP(15) = 4.65802970402487868693615238455D-3 C BCAP(16) = 0.0D0 C BCAP(17) = 0.0D0 C BDCAP(1) = 1.21278685171854149768890395495D-2 BDCAP(2) = 0.0D0 BDCAP(3) = 0.0D0 BDCAP(4) = 0.0D0 BDCAP(5) = 0.0D0 BDCAP(6) = 0.0D0 BDCAP(7) = 9.08394342270407836172412920433D-2 BDCAP(8) = 3.15683697648393399290429311645D-1 BDCAP(9) = -2.63224906576909737811077273181D-1 BDCAP(10) = 3.04780378618458886213892341513D-1 BDCAP(11) = -4.15516161554298332243867109382D-2 BDCAP(12) = 2.46775609676295306562750285101D-1 BDCAP(13) = 1.52260530105866022937951487642D-1 BDCAP(14) = 8.14384816302696075086493964505D-2 BDCAP(15) = 8.50257119389081128008018326881D-2 BDCAP(16) = -9.15518963007796287314100251351D-3 BDCAP(17) = 2.5D-2 C B(1) = 1.70087019070069917527544646189D-2 B(2) = 0.0D0 B(3) = 0.0D0 B(4) = 0.0D0 B(5) = 0.0D0 B(6) = 0.0D0 B(7) = 7.22593359308314069488600038463D-2 B(8) = 3.72026177326753045388210502067D-1 B(9) = -4.01821145009303521439340233863D-1 B(10) = 3.35455068301351666696584034896D-1 B(11) = -1.31306501075331808430281840783D-1 B(12) = 1.89431906616048652722659836455D-1 B(13) = 2.68408020400290479053691655806D-2 B(14) = 1.63056656059179238935180933102D-2 B(15) = 3.79998835669659456166597387323D-3 B(16) = 0.0D0 C B(17) = 0.0D0 C BD(1) = 1.70087019070069917527544646189D-2 BD(2) = 0.0D0 BD(3) = 0.0D0 BD(4) = 0.0D0 BD(5) = 0.0D0 BD(6) = 0.0D0 BD(7) = 7.60624588745593757356421093119D-2 BD(8) = 4.65032721658441306735263127583D-1 BD(9) = -5.35761526679071361919120311817D-1 BD(10) = 5.03182602452027500044876052344D-1 BD(11) = -2.62613002150663616860563681567D-1 BD(12) = 4.26221789886109468625984632024D-1 BD(13) = 1.07363208160116191621476662322D-1 BD(14) = 1.14139659241425467254626653171D-1 BD(15) = 6.93633866500486770090602920091D-2 BD(16) = 2.0D-2 C BD(17) = 0.0D0 C RETURN C C C c = 0 C 1 C C C 1 C c = -- C 2 50 C C C 1 C c = -- C 3 25 C C C 1 C c = -- C 4 10 C C C 2 C c = -- C 5 15 C C C 4 C c = -- C 6 25 C C C 1 C c = -- C 7 20 C C C 1 C c = - C 8 5 C C C 1 C c = - C 9 4 C C C 1 C c = - C 10 3 C C C 1 C c = - C 11 2 C C C 5 C c = - C 12 9 C C C 3 C c = - C 13 4 C C C 6 C c = - C 14 7 C C C 8437 C c = ---- C 15 8926 C C C c = 1 C 16 C C C c = 1 C 17 C C C 1 C a = ---- C 2, 1 5000 C C C 1 C a = ---- C 3, 1 3750 C C C 1 C a = ---- C 3, 2 1875 C C C 7 C a = ---- C 4, 1 2400 C C C 1 C a = - --- C 4, 2 240 C C C 1 C a = --- C 4, 3 160 C C C 2 C a = ---- C 5, 1 1215 C C C a = 0 C 5, 2 C C C 4 C a = --- C 5, 3 729 C C C 32 C a = ----- C 5, 4 18225 C C C 152 C a = ----- C 6, 1 78125 C C C a = 0 C 6, 2 C C C 1408 C a = ------ C 6, 3 196875 C C C 2048 C a = ------ C 6, 4 703125 C C C 432 C a = ------ C 6, 5 546875 C C C 29 C a = ----- C 7, 1 51200 C C C a = 0 C 7, 2 C C C 341 C a = ------ C 7, 3 387072 C C C 151 C a = - ------ C 7, 4 345600 C C C 243 C a = ------ C 7, 5 716800 C C C 11 C a = - ------ C 7, 6 110592 C C C 37 C a = ----- C 8, 1 12000 C C C a = 0 C 8, 2 C C C a = 0 C 8, 3 C C C 2 C a = ---- C 8, 4 1125 C C C 27 C a = ----- C 8, 5 10000 C C C 5 C a = ---- C 8, 6 3168 C C C 224 C a = ----- C 8, 7 20625 C C C 100467472123373 C a = ----------------- C 9, 1 27511470744477696 C C C a = 0 C 9, 2 C C C 101066550784375 C a = ----------------- C 9, 3 25488568483854336 C C C 49478218404275 C a = ----------------- C 9, 4 15475202293768704 C C C 21990175014231 C a = ---------------- C 9, 5 2674726322379776 C C C 3576386017671875 C a = - ------------------- C 9, 6 2723635603703291904 C C C 16163228153 C a = ------------- C 9, 7 1654104722787 C C C 38747524076705 C a = ----------------- C 9, 8 10316801529179136 C C C 62178936641284701329 C a = ----------------------- C 10, 1 16772293867250014666848 C C C a = 0 C 10, 2 C C C 46108564356250 C a = ---------------- C 10, 3 9072835168325103 C C C 1522561724950 C a = ---------------- C 10, 4 1296119309760729 C C C 45978886013453735443 C a = - ---------------------- C 10, 5 2174186242050927827184 C C C 299403512366617849203125 C a = ------------------------- C 10, 6 4981371278573254356053856 C C C 15571226634087127616 C a = --------------------- C 10, 7 774466927638876610083 C C C 133736375367792139885 C a = - ---------------------- C 10, 8 4717207650164066625051 C C C 7461389216 C a = ------------ C 10, 9 501451974639 C C C 501256914705531962342417557181 C a = -------------------------------- C 11, 1 14270506505142656332600844507392 C C C a = 0 C 11, 2 C C C 1143766215625 C a = - --------------- C 11, 3 132752960853408 C C C 6864570325 C a = - ------------- C 11, 4 1185294293334 C C C 194348369382310456605879163404183 C a = --------------------------------- C 11, 5 99893545535998594328205911551744 C C C 94634958447010580589908066176109375 C a = - ----------------------------------- C 11, 6 27549212808177898050085930321520256 C C C 17006472665356285286219618514 C a = - ------------------------------ C 11, 7 155584463413110817059022733377 C C C 33530528814694461893884349656345 C a = -------------------------------- C 11, 8 14270506505142656332600844507392 C C C 13439782155791134368 C a = - -------------------- C 11, 9 17777268379678341919 C C C 1441341768767571 C a = ----------------- C 11, 10 13159456712985856 C C C 105854110734231079069010159870911189747853 C a = ------------------------------------------- C 12, 1 5156624149476760916008179453333467046288864 C C C a = 0 C 12, 2 C C C 144579793509250000 C a = - -------------------- C 12, 3 19842290513127000261 C C C 101935644099967250 C a = - -------------------- C 12, 4 48188419817594143491 C C C 1585474394319811696785932424388196965 C a = ------------------------------------- C 12, 5 1709257457318830856936350991091849456 C C C 843499776333774172853009613469456309715703125 C a = - --------------------------------------------- C 12, 6 510505790798199330684809765880013237582597536 C C C 15057703799298260121553794369056896088480 C a = - ------------------------------------------ C 12, 7 714327132646734138085088291809720015274157 C C C 1749840442221344572962864758990584360232600 C a = ------------------------------------------- C 12, 8 1450300542040339007627300471250037606768743 C C C 11255775246405733991656178432768 C a = - -------------------------------- C 12, 9 27206626483067760480757659602193 C C C 669010348769579696 C a = ------------------- C 12, 10 7368057640845834597 C C C 4598083098752 C a = --------------- C 12, 11 858563707934367 C C C 1639758773684715326849438048667467886824967397 C a = - ----------------------------------------------- C 13, 1 11447568726280607813664651120965112496134881280 C C C a = 0 C 13, 2 C C C 3942453384375 C a = --------------- C 13, 3 314673684985856 C C C 11737114158175 C a = ---------------- C 13, 4 1719466921529856 C C C 23710715033675876683332701739887457 C a = - ----------------------------------- C 13, 5 4940189888325748664958546898558976 C C C 498150575499633273684774666731162498301909124515625 C a = --------------------------------------------------- C 13, 6 87415924307623977386706008889913792042985180430336 C C C 64881557768202140428371179540010005713998551 C a = -------------------------------------------- C 13, 7 85896810580242200654071863296887242202224768 C C C 2336309182318568698279006266321563486172654055 C a = - ----------------------------------------------- C 13, 8 18316109962048972501863441793544179993815810048 C C C 493399374030747471036018890494175 C a = - --------------------------------- C 13, 9 251658285736841065236836942273664 C C C 418285003077108927126515545155 C a = ------------------------------ C 13, 10 455369916679568501838710898688 C C C 15171723902781457 C a = - ----------------- C 13, 11 63532954684873728 C C C 1501203688494867 C a = ---------------- C 13, 12 9434957026426880 C C C 34188549803371802849576690267872548602326398788953 C a = -------------------------------------------------- C 14, 1 42496542183406636759747616530102745233754251202880 C C C a = 0 C 14, 2 C C C 18971246281693750 C a = - ------------------- C 14, 3 1138830954584356089 C C C 59230464334542700 C a = - ------------------- C 14, 4 2765732318276293359 C C C 5147939981309774383134903239728881770043 C a = ---------------------------------------- C 14, 5 305929030949718561059100251282184099064 C C C 362572021355026772337065830211467821556305840522907812 C a = - ------------------------------------------------------ C 14, 6 324512095420929759624784749347170583153994213035432256 C C C 60305503318319653518547439098565661266182518307816 C a = - -------------------------------------------------- C 14, 7 17856872599361492097414471889911176856851308259643 C C C 1036461878759982363277481306266144563833492657780645 C a = - ---------------------------------------------------- C 14, 8 67994467493450618815596186448164392374006801924608 C C C 128398681100219349205889126776607047000 C a = --------------------------------------- C 14, 9 7473801441221286756994805323613917077 C C C 49156374556350058671822606102117 C a = - -------------------------------- C 14, 10 9039888303968618912866414995904 C C C 12253036339964386945 C a = -------------------- C 14, 11 8828680926314891943 C C C 647188390508758231059 C a = - ---------------------- C 14, 12 1092148506009694282240 C C C 10915833599872 C a = --------------- C 14, 13 368729913707897 C C Ca = C 15, 1 C C 4939337286263213195547765488387521892799075623007291241961609516532 C - ------------------------------------------------------------------- C 5408250052307451520718178852915698257207815452080611897685945761264 C C C a = 0 C 15, 2 C C C 7588799849596321243074032368290625 C a = ------------------------------------- C 15, 3 3147217749590114939838670370597819616 C C C 16870665568420512953501332587233725 C a = ------------------------------------ C 15, 4 955405388268427749593882076788623812 C C C 80864251591837801485030858227147601466956843757908779606 C a = - -------------------------------------------------------- C 15, 5 54447992506702009927986632715967769032585338753056786562 C C Ca = C 15, 6 C C 4610328329649866588704236006423149172472141907645890762410296050212 C ------------------------------------------------------------------- C 2135428689710103309390449198881479603148467934048051598947383737508 C C Ca = C 15, 7 C C 4159963831215576225909381034291748993887819834160487158570788681 C ---------------------------------------------------------------- C 1040533184037697645660563795162185415624171583014576682740416336 C C Ca = C 15, 8 C C 738139214212435127943380193414870655354213707189052136566460666444958 C --------------------------------------------------------------------- C 259596002510757672994472584939953516345975141699869371088925396540699 C C C 333683433458405281346882867597135977469443722954786270692 C a = - --------------------------------------------------------- C 15, 9 132102862435303266640535426836147775872819092781208127980 C C C 42661937996741208687503901295747546613008142604821349179 C a = -------------------------------------------------------- C 15, 10 55162410119399855550108207148248549410926885937244965785 C C C 630755628691078947314733435975762542732598947 C a = - --------------------------------------------- C 15, 11 333503232300511886435069380727586592765317456 C C C 1522350657470125698997653827133798314909646891 C a = ---------------------------------------------- C 15, 12 1520094067152619944607524353149267399623188480 C C C 305575414262755427083262606101825880 C a = -------------------------------------- C 15, 13 65839748482572312891297405431209259829 C C C 256624643108055110568255672032710477795 C a = ----------------------------------------- C 15, 14 22874609758516552135947898572671559986304 C C C 571597862947184314270186718640978947715678864684269066846 C a = - --------------------------------------------------------- C 16, 1 207705506488030390761613596901272001190776700439774478634 C C C a = 0 C 16, 2 C C C 66981514290625 C a = ---------------- C 16, 3 1829501741761029 C C C 43495576635800 C a = ---------------- C 16, 4 4443075658562499 C C C 127865248353371207265315478623656127 C a = - ------------------------------------ C 16, 5 10401415428935853634424440540325344 C C C 131656514265807573955723157408023481433806699348396032656 C a = --------------------------------------------------------- C 16, 6 92668695535091962564795912774190176478892159517481612467 C C C 3881494143728609118531066904799685950051960514138645179820 C a = ---------------------------------------------------------- C 16, 7 2446349095978358868919950548516272963929118212742344026549 C C C 16292266704968075585259245375842819400619822954470178684291 C a = ----------------------------------------------------------- C 16, 8 66288722243155885736983218667976563740242178853010092663614 C C C 43986024977384568043684084266385512680544563954 C a = - ----------------------------------------------- C 16, 9 4922783599524658241955780540171948284522386185 C C C 285912200202585226675651763671663063668290787 C a = --------------------------------------------- C 16, 10 65371192072964016939690070594254881767827200 C C C 6776815256667778089672518929 C a = - ---------------------------- C 16, 11 3693654613173093729492918708 C C C 398946554885847045598775476868169 C a = --------------------------------- C 16, 12 344154261237450078839899047372800 C C C 76630698033396272 C a = - ------------------- C 16, 13 4432017119727044925 C C C 28401702316003037 C a = ------------------- C 16, 14 1469612686944417840 C C C 66049942462586341419969330578128801 C a = -------------------------------------- C 16, 15 12691068622536592094919763114637498325 C C C 83940754497395557520874219603241359529066454343054832302344735 Ca = -------------------------------------------------------------- C 17, 1 64192596456995578553872477759926464976144474354415663868673233 C C C a = 0 C 17, 2 C C C 892543892035485503125 C a = ----------------------- C 17, 3 51401651664490002607536 C C C 12732238157949399705325 C a = - ------------------------ C 17, 4 686579204375687891972088 C C C 5290376174838819557032232941734928484252549 C a = ------------------------------------------- C 17, 5 357179779572898187570048915214361602000384 C C C 2687322933801750693719999180471745666665021538793817303193221 C a = ------------------------------------------------------------- C 17, 6 2863980005760296740624015421425947092438943496681472214589916 C C Ca = C 17, 7 C C 197649786681880330585741729796159873563741413724149351549277865 C - --------------------------------------------------------------- C 378029217824623393200881653405474359138017953416246216408422692 C C Ca = C 17, 8 C C 100286075630483975704018828319990067604207336241794360144098685695 C - ------------------------------------------------------------------ C 20486915674765670626893195919603679319429068544972409068469849579 C C C 8739866119696575810411768434844068608106287881671139259 C a = ------------------------------------------------------- C 17, 9 2282122412587168891929052689609009868137678763277087160 C C C 7922242431969626895355493632206885458496418610471389 C a = - ---------------------------------------------------- C 17, 10 748272134517487495468365669337985635214015258726400 C C C 2777643183645212014464950387658055285 C a = ------------------------------------- C 17, 11 1141545470045611737197667093465955392 C C C 1372659703515496442825084239977218110461 C a = - ---------------------------------------- C 17, 12 1313121960368535725613950174847107891200 C C C 6144417902699179309851023 C a = -------------------------- C 17, 13 85608793932459282773805825 C C C 140294243355138853053241 C a = -------------------------- C 17, 14 64884622846351585391642880 C C C 168671028523891369934964082754523881107337 C a = -------------------------------------------- C 17, 15 24062875279623260368388427013982199424119600 C C C a = 0 C 17, 16 C C C 63818747 C bc = ---------- C 1 5262156900 C C C bc = 0 C 2 C C C bc = 0 C 3 C C C bc = 0 C 4 C C C bc = 0 C 5 C C C bc = 0 C 6 C C C 22555300000000 C bc = --------------- C 7 261366897038247 C C C 1696514453125 C bc = ------------- C 8 6717619827072 C C C 45359872 C bc = - --------- C 9 229764843 C C C 19174962087 C bc = ----------- C 10 94371046000 C C C 19310468 C bc = - --------- C 11 929468925 C C C 16089185487681 C bc = --------------- C 12 146694672924800 C C C 1592709632 C bc = ----------- C 13 41841694125 C C C 52675701958271 C bc = ---------------- C 14 4527711056573100 C C C 12540904472870916741199505796420811396 C bc = ---------------------------------------- C 15 2692319557780977037279406889319526430375 C C C bc = 0 C 16 C C C bc = 0 C 17 C C C 63818747 C bdc = ---------- C 1 5262156900 C C C bdc = 0 C 2 C C C bdc = 0 C 3 C C C bdc = 0 C 4 C C C bdc = 0 C 5 C C C bdc = 0 C 6 C C C 451106000000000 C bdc = ---------------- C 7 4965971043726693 C C C 8482572265625 C bdc = -------------- C 8 26870479308288 C C C 181439488 C bdc = - --------- C 9 689294529 C C C 57524886261 C bdc = ------------ C 10 188742092000 C C C 38620936 C bdc = - --------- C 11 929468925 C C C 144802669389129 C bdc = --------------- C 12 586778691699200 C C C 6370838528 C bdc = ----------- C 13 41841694125 C C C 368729913707897 C bdc = ---------------- C 14 4527711056573100 C C C 111940113324845802831946788738852162520696 C bdc = ------------------------------------------- C 15 1316544263754897771229629968877248424453375 C C C 113178587 C bdc = - ----------- C 16 12362232960 C C C 1 C bdc = -- C 17 40 C C C 27121957 C b = ---------- C 1 1594593000 C C C b = 0 C 2 C C C b = 0 C 3 C C C b = 0 C 4 C C C b = 0 C 5 C C C b = 0 C 6 C C C 4006163300000 C b = -------------- C 7 55441463008113 C C C 9466403125 C b = ----------- C 8 25445529648 C C C 163199648 C b = - --------- C 9 406149975 C C C 23359833 C b = -------- C 10 69636250 C C C 18491714 C b = - --------- C 11 140828625 C C C 11052304606701 C b = -------------- C 12 58344472186000 C C C 1191129152 C b = ----------- C 13 44377554375 C C C 2033811086741 C b = --------------- C 14 124730332137000 C C C 3616943474975740389660406409450169802 C b = --------------------------------------- C 15 951830146690244407118982233597812374375 C C C b = 0 C 16 C C C b = 0 C 17 C C C 27121957 C bd = ---------- C 1 1594593000 C C C bd = 0 C 2 C C C bd = 0 C 3 C C C bd = 0 C 4 C C C bd = 0 C 5 C C C bd = 0 C 6 C C C 4217014000000 C bd = -------------- C 7 55441463008113 C C C 47332015625 C bd = ------------ C 8 101782118592 C C C 652798592 C bd = - ---------- C 9 1218449925 C C C 70079499 C bd = --------- C 10 139272500 C C C 36983428 C bd = - --------- C 11 140828625 C C C 99470741460309 C bd = --------------- C 12 233377888744000 C C C 4764516608 C bd = ----------- C 13 44377554375 C C C 14236677607187 C bd = --------------- C 14 124730332137000 C C C 198066487470143918516004831967805004004 C bd = ---------------------------------------- C 15 2855490440070733221356946700793437123125 C C C 1 C bd = -- C 16 50 C C C bd = 0 C 17 C END C