SUBROUTINE DINTM (NDIMI,ANSWER,WORK,NWORK,IOPT) c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. c ALL RIGHTS RESERVED. c Based on Government Sponsored Research NAS7-03001. C>> 2009-11-03 DINTM Krogh Initialized lots of variables. C>> 1996-03-31 DINTM Krogh Removed unused variable in common. c>> 1995-11-20 DINTM Krogh Converted from SFTRAN to Fortran 77. C>> 1994-11-14 DINTM Krogh Declared all vars. c>> 1994-10-19 DINTM Krogh Changes to use M77CON c>> 1994-07-07 DINTM Snyder set up for CHGTYP. C>> 1993-05-18 DINTM Krogh -- Changed "END" to "END PROGRAM" C>> 1994-05-02 DINTM Snyder corrected some comments C>> 1991-09-20 DINTM Krogh converted '(1)' dimensioning to '(*)'. C>> 1987-11-19 DINTM Snyder Initial code. c c--D replaces "?": ?INTA,?intc,?intec,?INTF,?INTM,?INTMA,?INTNC,?INTOP C C ****************************************************************** C C THIS SUBROUTINE ATTEMPTS TO CALCULATE THE INTEGRAL OF A FUNCTION C PROVIDED BY THE USER. THE FUNCTION IS PROVIDED BY THE USER VIA A C SUBROUTINE REFERENCED BY CALL DINTF(F,X,IFLAG) OR BY REVERSE C COMMUNICATION. ALL ABSCISSAE ARE WITHIN THE REGION OF C INTEGRATION. C C THE RESULT IS OBTAINED USING QUADRATURE FORMULAE DUE TO C T. N. L. PATTERSON, MATHEMATICS OF COMPUTATION, VOLUME 22, C PAGES 847-856, 1968. C C ***** WARNING *********************************************** C C THE RELIABILITY AND EFFICIENCY OF THIS PROGRAM ARE STRONGLY C INFLUENCED BY DISCONTINUITIES IN THE FUNCTION OR IT'S C DERIVATIVES, INTEGRABLE SINGULARITIES IN THE REGION OF C INTEGRATION, AND NON-INTEGRABLE SINGULARITIES NEAR THIS REGION. C (INCLUDING COMPLEX VALUES). THE EFFICIENCY AND RELIABILITY C OF INTEGRATING SUCH FUNCTIONS MAY BE GREATLY IMPROVED BY MANUALLY C SUBDIVIDING THE REGION OF INTEGRATION AT THE DISCONTINUITY, C SINGULARITY OR CLOSEST POINT TO A POLE, AND SUMMING THE ANSWERS. C A CHANGE OF VARIABLE TO ELIMINATE OR REDUCE THE STRENGTH OF THE C SINGULARITY WILL SIGNIFICANTLY IMPROVE PERFORMANCE. C C ***** FORMAL ARGUMENTS ************************************* C C NDIMI IS THE NUMBER OF DIMENSIONS OF INTEGRATION. C ANSWER IS THE ESTIMATE OF THE INTEGRAL. WHEN REVERSE COMMUNICATION C IS SPECIFIED IT IS USED TO PASS FUNCTION VALUES FROM C THE USER PROGRAM TO DINTA OR DINTMA. DOUBLE PRECISION ANSWER C WORK CONTAINS THE LIMITS, WORKING SPACE AND MAY BE REFERENCED BY C THE OPTION VECTOR (SEE IOPT BELOW). WHEN THE INTEGRATION C IS COMPLETE, WORK(1) CONTAINS THE ESTIMATED ABSOLUTE ERROR. C WHEN REVERSE COMMUNICATION IS SPECIFIED, WORK(1) IS USED TO C PASS ABSCISSAE FROM DINTA OR DINTMA TO THE USER PROGRAM. C WORK(1), ..., WORK(NDIMI) ARE COMPONENTS OF THE VECTOR OF C ABSCISSAE. WORK(NDIMI+1), ..., WORK(2*NDIMI) ARE LOWER C LIMITS. WORK(2*NDIMI+1), ..., WORK(3*NDIMI) ARE UPPER LIMITS. C THE ABSCISSAE AND LIMITS ARE STORED INNERMOST FIRST. DINTF C IS CALLED TO REQUEST THE LIMITS OF EVERY DIMENSION, BUT C CONSTANT LIMITS (THE LIMITS OF THE OUTER DIMENSION ARE ALWAYS C CONSTANT) MAY INSTEAD BE STORED BEFORE INTEGRATION BEGINS. C WORK(3*NDIMI+1), ..., WORK(3*NDIMI+KWORK*(NDIMI-1)) ARE C WORKING STORAGE, WHERE KWORK DEPENDS ON THE MACHINE AND THE C PRECISION OF THE PROGRAM VERSION. KWORK MUST SPECIFY ENOUGH C STORAGE FOR 4 VARIABLES WHICH ARE ALWAYS DOUBLE PRECISION, C 139 VARIABLES WHICH ARE DOUBLE PRECISION IN THE DOUBLE C PRECISION VERSION, 30 VARIABLES WHICH ARE DOUBLE PRECISION C IN THE DOUBLE PRECISION PROGRAM IF THE LARGEST REPRESENTABLE C NUMBERS IN SINGLE AND DOUBLE PRECISION ARE DIFFERENT BY MORE C THAN A FACTOR OF 4, 29 INTEGER VARIABLES AND 11 LOGICAL C VARIABLES. REMEMBER TO CONSIDER THE DATA TYPE OF WORK WHEN C SPECIFYING THE VALUE OF KWORK. EXAMPLES ARE C-- Begin mask code changes C IF PRECISION IS REAL, AND INTEGER AND LOGICAL VARIABLES C OCCUPY THE SAME SPACE AS REAL VARIABLES, KWORK=217. C IF PRECISION IS REAL, AND INTEGER AND LOGICAL VARIABLES C OCCUPY HALF AS MUCH SPACE AS REAL VARIABLES, KWORK=? C IF PRECISION IS PARTLY DOUBLE, AND INTEGER AND LOGICAL C VARIABLES OCCUPY THE SAME SPACE AS REAL VARIABLES, C KWORK=? C IF PRECISION IS PARTLY DOUBLE, AND INTEGER AND LOGICAL C VARIABLES OCCUPY HALF AS MUCH SPACE AS REAL VARIABLES, C KWORK=? C IF PRECISION IS DOUBLE, AND INTEGER AND LOGICAL VARIABLES C OCCUPY THE SAME SPACE AS REAL VARIABLES, KWORK=193. C IF PRECISION IS DOUBLE, AND INTEGER AND LOGICAL VARIABLES C OCCUPY HALF AS MUCH SPACE AS REAL VARIABLES, KWORK=? C FOR PORTABILITY, ASSUME THE WORST CASE. THAT IS, IF THE C PRECISION IS SINGLE, KWORK=217. IF THE PRECISION IS DOUBLE, C KWORK=193. C WORK(3*NDIMI+KWORK*(NDIMI-1)+1), ... ARE AVAILABLE FOR C REFERENCE THROUGH THE OPTION VECTOR. C-- End mask code changes DOUBLE PRECISION WORK(*) C NWORK IS THE NUMBER OF ELEMENTS ALLOCATED FOR WORK. C NWORK MUST BE GREATER THAN 3*NDIMI+KWORK*(NDIMI-1). C IOPT IS A VECTOR OF INTEGERS USED TO RETURN THE STATUS AND TO C SPECIFY OPTIONS. INTEGER IOPT(*) C IOPT(1) IS USED TO RETURN A STATUS INDICATOR TO THE USER. C VALUES OF THIS INDICATOR ARE C -NDIM - NORMAL TERMINATION, EITHER THE ABSOLUTE OR THE C RELATIVE ERROR CRITERIA ARE SATISFIED. C -NDIM-1 - NORMAL TERMINATION, NEITHER THE ABSOLUTE NOR THE C RELATIVE ERROR CRITERIA ARE SATISFIED, BUT THE C ERROR RELATIVE TO THE OBTAINABLE PRECISION C PRECISION CRITERION IS SATISFIED. C -NDIM-2 - NORMAL TERMINATION, BUT NONE OF THE ERROR CRITERIA C ARE SATISFIED. C -NDIM-3 - NOT ENOUGH SPACE IN WORK, NWORK IS TOO SMALL. C -NDIM-4 - BAD IOPT VALUE. C -NDIM-5 - TOO MANY FUNCTION VALUES NEEDED. C -NDIM-KDIM-5 - APPARENT NON-INTEGRABLE SINGULARITY IN C DIMENSION KDIM. ANSWER CONTAINS THE APPROXIMATE C ABSCISSA OF THE SINGULARITY IN THE KDIM-TH C DIMENSION, WORK(KDIM+1) THROUGH WORK(NDIM) C CONTAIN THE ABSCISSAE OF EXTERIOR INTEGRALS. C -2*NDIM-KDIM-5 - INCORRECT INNER INTEGRAL DIMENSIONALITY C SPECIFIED DURING EVALUATION OF INTEGRAL AT C DIMENSION KDIM, THAT IS, IABS(IFLAG(1)) = KDIM. C ENTRIES IN IOPT STARTING WITH IOPT(2) ARE DESCRIBED BELOW. C IOPT(I) ENTRY IN IOPT(I) MEANS C 0 NO MORE OPTIONS, BEGIN INTEGRATION. C 1 IOPT(I+1) CONTAINS THE UNIT NUMBER FOR OUTPUT. C 2 IOPT(I+1) IS AN NDIMI DIGIT INTEGER, WHERE C EACH DIGIT IS THE DIAGNOSTIC PRINT LEVEL FOR C ONE DIMENSION. THE LOW ORDER DIGIT IS THE PRINT C LEVEL FOR THE INNER DIMENSION. C 0 - NO PRINTING C 1 - MINIMUM PRINTING - ERROR MESSAGES (DEFAULT) C 2 - PANEL BOUNDARIES AND ANSWERS C 3 - ERROR ESTIMATES FOR EACH QUADRATURE FORMULA C 4 - DETAILED OUTPUT (DIFFERENCE LINES, ETC). C 3 WORK(IOPT(I+1)) CONTAINS THE ABSOLUTE TOLERANCE, C WORK(IOPT(I+1)+1) CONTAINS TOLERANCE RELATIVE TO C THE LOCALLY OBTAINABLE PRECISION, AND C WORK(IOPT(I+1)+2) CONTAINS TOLERANCE RELATIVE TO C THE VALUE OF THE INTEGRAL. THE TOLERANCE RELATIVE C TO THE LOCALLY OBTAINABLE PRECISION IS SPECIFIED AS C THE FRACTION OF LOCALLY OBTAINABLE DIGITS THAT ARE C PERMITTED TO BE WRONG, AND IS INTERNALLY BOUNDED C BETWEEN 0.0 AND 1.0. IF ANY OF THE ERROR CONTROL C CRITERIA ARE SATISFIED, THE ANSWER IS ACCEPTED. IF C THIS OPTION IS NOT USED, THE ABSOLUTE AND RELATIVE C TOLERANCES ARE SET TO ZERO, AND THE TOLERANCE C RELATIVE TO THE LOCALLY OBTAINABLE PRECISION IS SET C TO 0.25. C 4 WORK(IOPT(I+1)) CONTAINS ABSOLUTE ERROR COMMITTED C COMPUTING F. WORK(IOPT(I+1)) MAY BE CHANGED DURING C THE INTEGRATION. C 5 WORK(IOPT(I+1)) CONTAINS RELATIVE ERROR EXPECTED TO C BE COMMITTED COMPUTING F. CHANGES TO THIS VALUE C DURING INTEGRATION WILL NOT BE DETECTED. C 6 USE REVERSE COMMUNICATION - C CALL DINTM (NDIMI,ANSWER,WORK,NWORK,IOPT) C 1 CALL DINTMA (ANSWER,WORK,IOPT(1)) C IF (IOPT(1)) 3,2,4 C 2 COMPUTE THE INNERMOST INTEGRAND - C ANSWER=F(WORK(K)), K=1,NDIMI C THEN EITHER C GO TO 1 C OR (FOR BETTER EFFICIENCY) C CALL DINTA (ANSWER,WORK,IOPT(1)) C IF (IOPT(1)) 1,2,5 C VALUES OF IOPT(1) PRODUCED BY DINTA AND DINTMA C HAVE DIFFERENT MEANINGS. C 3 IF (IOPT(1)+NDIMI.LE.0) GO TO 6 C THE INTEGRAL OVER THE (-IOPT(1))TH DIMENSION IS C CONTAINED IN ANSWER. IF A TRANSFORMATION OF C THE INTEGRAL IS NOT NECESSARY BEFORE IT IS USED C AS AN INTEGRAND FOR THE (1-IOPT(1))TH DIMENSION C GO TO 1 C ELSE COMPUTE ANSWER=G(ANSWER,WORK(K)), C K=1-IOPT(1),NDIMI AND G IS THE TRANSFORMATION. C COMPUTE ALSO WORK(1)=PARTIAL DERIVATIVE OF G C WITH RESPECT TO ANSWER THEN C GO TO 1 C 4 COMPUTE WORK(NDIMI+IOPT(1))=A(WORK(K)), C K=IOPT(1)+1,NDIMI, WORK(2*NDIMI+IOPT(1))=B(WORK(K)) C K=IOPT(1)+1,NDIM, WHERE A IS THE LOWER LIMIT OF C INTEGRATION FOR THE (IOPT(1))TH DIMENSION, AND B IS C THE UPPER LIMIT. IF A TRANSFORMATION AS DISCUSSED C ABOVE IS TO BE APPLIED WHEN THE INTEGRATION OF THIS C DIMENSION IS COMPLETE, COMPUTE WORK(1)=AN ESTIMATE C OF THE UPPER BOUND OF THE PARTIAL DERIVATIVE OF G C WITH RESPECT TO THE INTEGRAL. IF G IS LINEAR IN C THE INTEGRAL, THIS DERIVATIVE CAN BE COMPUTED C WITHOUT KNOWING THE INTEGRAL, RATHER THAN C ESTIMATED. IF THERE IS A SINGULARITY IN THE C INTEGRAND FOR THIS DIMENSION (PERHAPS INTRODUCED C BY THE BOUNDARY OR THE TRANSFORMATION G), SET C IOPT(1) TO THE LOCATION IN WORK OF THE ABSCISSA IN C THIS DIMENSION OF THE SINGULARITY, AND WITH THE C SIGN SET AS DESCRIBED FOR OPTION 11 BELOW. THE C MAGNITUDE OF IOPT(1) MUST BE GREATER THAN NDIM FOR C THIS REQUEST TO BE DETECTED. IN PARTICULAR, IF C THERE IS A SINGULARITY AT ONE OF THE LIMITS, SET C IOPT(1) TO IOPT(1)+NDIM OR -(IOPT(1)+2*NDIM). THEN C GO TO 1. C 5 IOPT(1)=-(IOPT(1)+NDIM) C 6 CONTINUE C AT THIS POINT, IERR CONTAINS THE STATUS DESCRIBED C BELOW, FOR OPTION 8. C 7 SET MINIMUM INDEX OF QUADRATURE FORMULA TO USE C BEFORE SUBDIVISION TO IOPT(I+1). C 8 NO EFFECT. IOPT(I+1) MAY BE USED TO PASS C INFORMATION TO DINTF. INCREMENT I BY 2. C 9 IOPT(I+1) IS THE MAXIMUM NUMBER OF FUNCTION VALUES C TO USE TO INTEGRATE THE ENTIRE PROBLEM. IF C IOPT(I+1) .LE. 0, THE NUMBER OF FUNCTION VALUES C IS NOT CONTROLLED. C 10 IOPT(I+1) IS USED TO RETURN THE NUMBER OF FUNCTION C VALUES USED TO INTEGRATE THE ENTIRE PROBLEM. C 11 WORK(IABS(IOPT(I+1))) IS THE LOCATION OF A C SINGULARITY OR DISCONTINUITY IN THE OUTERMOST C INTEGRAL. IF THE LOCATION IS INSIDE THE INTERVAL, C THE INTERVAL IS IMMEDIATELY SUBDIVIDED. IF C IOPT(I+1) .GT. 0, A T**2 SUBSTITUTION BASED AT C WORK(...) WILL BE USED. IF IOPT(I+1) .LT. 0, C A T**4 SUBSTITUTION WILL BE USED. IF IOPT(I+1) C .EQ. 0, NO SUBSTITUTION WILL BE USED. TO NOTIFY C THE PROGRAM OF SINGULARITIES OR DISCONTINUITIES IN C INTERIOR INTEGRANDS, SET IFLAG NEGATIVE WHEN C ASKED TO COMPUTE LIMITS. SEE THE DESCRIPTION C OF REVERSE COMMUNICATION ABOVE, OR THE DESCRIPTION C OF DINTF BELOW. C 12 WORK(IOPT(I+1)) CONTAINS THE ABSOLUTE ERROR IN THE C LOWER LIMIT OF THE CURRENT DIMENSION. THE ERROR IN C THE UPPER LIMIT IS IN WORK(IOPT(I+1)+1). C 13 IOPT(I+1) IS THE LOCATION IN IOPT IN WHICH THE USER C IS TO STORE NON-STANDARD CHANGES TO THE DIMENSION OF C INNER INTEGRALS. THE DEFAULT VALUE IS 1, ALLOWING C SUCH NONSTANDARD CHANGES TO BE STORED IN THE CELL OF C IFLAG (=IOPT) USUALLY USED FOR COMMUNICATION WITH C DINTF, IN THE CASE WHEN IT IS NOT SIMULTANEOUSLY C NECESSARY TO INDICATE A SINGULARITY. SEE THE C DESCRIPTION OF THE INTERFACE TO DINTF BELOW. C C ALL OPTIONS ARE SET TO THEIR NOMINAL VALUES BEFORE THE OPTION C VECTOR IS PROCESSED. C C ***** DINTF *********************************************** C C DINTF IS REFERENCED VIA CALL DINTF (ANSWER,WORK,IFLAG). C WORK AND ANSWER ARE AS DESCRIBED ABOVE. VALUES OF IFLAG ARE C IFLAG .LT. 0 C THE INTEGRAL OVER THE (-IFLAG)TH DIMENSION IS C CONTAINED IN ANSWER. IF A TRANSFORMATION OF C THE INTEGRAL IS NOT NECESSARY BEFORE IT IS USED C AS AN INTEGRAND FOR THE (1-IFLAG)TH DIMENSION C RETURN C ELSE IF THE INTEGRAND DEPENDS ON ONLY ONE INNER C INTEGRAL COMPUTE ANSWER=G(ANSWER,WORK(K)), C ALSO MULTIPLY WORK(1) BY THE PARTIAL DERIVATIVE C OF G WITH RESPECT TO ANSWER THEN C RETURN C ELSE IF THE INTEGRAND DEPENDS ON MORE INTEGRALS C SAVE ANSWER AND WORK(1), AND SET IFLAG(IXKDIM) TO C THE DIMENSION OF INTEGRAL NEXT TO BE COMPUTED C (IXKDIM IS THE VALUE SPECIFIED BY OPTION 12, OR 1 C IF OPTION 12 IS NOT SELECTED), THEN C RETURN C ELSE (ALL INNER INTEGRALS UPON WHICH THE INNER C INTEGRAND DEPENDS HAVE BEEN EVALUATED) CALCULATE C ANSWER=G(ANSWER,WORK(K),SAVED INTEGRALS) C K=1-IFLAG, NDIMI (G IS THE TRANSFORMATION), C MULTIPLY WORK(1) BY THE PARTIAL DERIVATIVE OF G C WITH RESPECT TO THE LAST INTEGRAL, AND ADD THE C PRODUCT OF EACH OF THE SAVED VALUES OF WORK(1) AND C PARTIAL DERIVATIVE OF G WITH RESPECT TO THE C CORRESPONDING SAVED INTEGRAL ONTO WORK(1). THAT C IS, WORK(1) IS THE TOTAL ERROR IN G, CALCULATED AS C THE INNER PRODUCT OF THE ERRORS IN THE INDIVIDUAL C ARGUMENTS, AND THE PARTIAL DERIVATIVES OF G WITH C RESPECT TO THOSE ARGUMENTS. C IFLAG .EQ. 0 C COMPUTE THE INNERMOST INTEGRAND - C ANSWER=F(WORK(K)), K=1,NDIMI C AND C RETURN C IFLAG .GT. 0 C COMPUTE WORK(NDIMI+IFLAG)=A(WORK(K)), C K=IFLAG+1,NDIMI, WORK(2*NDIMI+IFLAG)=B(WORK(K)), C K=IFLAG+1,NDIM, WHERE A IS THE LOWER LIMIT OF C INTEGRATION FOR THE (IFLAG)TH DIMENSION, AND B IS C THE UPPER LIMIT. IF A TRANSFORMATION AS DISCUSSED C ABOVE IS TO BE APPLIED WHEN THE INTEGRATION OF THIS C DIMENSION IS COMPLETE, COMPUTE WORK(1)=AN ESTIMATE C OF THE UPPER BOUND OF THE PARTIAL DERIVATIVE C OF G WITH RESPECT TO THE INTEGRAL. IF G IS C LINEAR IN THE INTEGRAL, THIS DERIVATIVE MAY BE C COMPUTED WITHOUT KNOWING THE INTEGRAL, RATHER C THAN ESTIMATED. IF THERE IS A SINGULARITY IN C THE INTEGRAND FOR THIS DIMENSION (PERHAPS C INTRODUCED BY THE BOUNDARY OR THE TRANSFORMATION), C SET IFLAG TO THE LOCATION IN WORK OF THE ABSCISSA C IN THIS DIMENSION OF THE SINGULARITY, AND SET THE C SIGN AS FOR OPTION 11 ABOVE. NOTE THAT THE C MAGNITUDE OF IFLAG MUST BE GREATER THAN NDIM FOR C THIS REQUEST TO BE DETECTED. IN PARTICULAR, IF C THERE IS A SINGULARITY AT ONE OF THE LIMITS, SET C THE MAGNITUDE OF IFLAG TO IFLAG+NDIM OR C IFLAG+2*NDIM, WITH THE SIGN AS DESCRIBED FOR C OPTION 11. THEN C RETURN. C IF THE FIRST INTEGRAL TO BE COMPUTED DOES NOT HAVE C DIMENSION IFLAG(1)-1, SET IFLAG(IXKDIM) TO THE C DIMENSION OF INTEGRAL TO BE NEXT COMPUTED. SEE C OPTION 12 AND THE DISCUSSION OF IFLAG .LT. 0 ABOVE C FOR INFORMATION ON IXKDIM. C C ***** EXTERNAL REFERENCES ********************************* C C DINTMA TO DO THE INTEGRATION. C DINTOP TO SET OPTIONS. C C ***** LOCAL VARIABLES ************************************* C INTEGER NDIMI, NWORK C C ***** COMMON STORAGE ************************************** C C COMMON /DINTNC/ CONTAINS VARIABLES NOT SEPARATELY SAVED FOR C EACH DIMENSION OF A MULTIPLE QUADRATURE. COMMON /DINTC/ C CONTAINS VARIABLES THAT MUST BE SAVED FOR EACH DIMENSION OF THE C QUADRATURE. THE VARIABLES IN EACH COMMON BLOCK ARE STORED IN THE C ORDER - ALWAYS DOUBLE, DOUBLE IF DOUBLE PRECISION PROGRAM, DOUBLE C IF DOUBLE PRECISION PROGRAM AND EXPONENT RANGE OF DOUBLE AND C SINGLE VERY DIFFERENT, SINGLE, INTEGER, LOGICAL. A PAD OF LOGICAL C VARIABLES IS INCLUDED AT THE END OF /DINTC/. THE DIMENSION OF C THE PAD MAY NEED TO BE VARIED SO THAT NO VARIABLES BEYOND THE END C OF THE COMMON BLOCK ARE ALTERED. C C DECLARATIONS OF COMMON /DINTNC/ VARIABLES. C DOUBLE PRECISION AINIT, BINIT, FNCVAL, S, TP DOUBLE PRECISION FER, FER1, RELOBT, TPS, XJ, XJP INTEGER FEA, FEA1, INC, INC2, IPRINT, 1 ISTOP(2,2),JPRINT, KDIM, KK, KMAXF, NDIM, 2 NFINDX, NFMAX, NFMAXM, RELTOL, REVERM, REVERS, 3 WHEREM LOGICAL NEEDH C C DECLARATIONS OF COMMON /DINTC/ VARIABLES. C c--D Next line special: S => D, X => Q, D => D, P => D DOUBLE PRECISION ACUM, PACUM, RESULT(2) C 139 $.TYPE.$ VARIABLES DOUBLE PRECISION 1 AACUM, ABSCIS, DELMIN, DELTA, DIFF, DISCX(2), 2 END(2), ERRINA, ERRINB, FAT(2), FSAVE, 3 FUNCT(24), F1, F2, LOCAL(4), PAACUM, PF1, 4 PF2, PHISUM, PHTSUM, PX, SPACE(6), 5 STEP(2), START(2), SUM, T, TA, TASAVE, 6 TB, TEND, WORRY(2), X, X1, 7 X2, XT(17), FT(17), PHI(34) c Note XT, FT, and PHI above are last, because they must be in adjacent c locations in DINTC. C 30 $DSTYP$ VARIABLES DOUBLE PRECISION 1 ABSDIF, COUNT, EDUE2A, EDUE2B, EP, EPNOIZ, 2 EPS, EPSMAX, EPSMIN, EPSO, EPSR, EPSS, 3 ERR, ERRAT(2), ERRC, ERRF, ERRI, ERRT(2), 4 ESOLD, EXTRA, PEPSMN, RE, RELEPS, REP, 5 REPROD, RNDC, TLEN, XJUMP C 29 INTEGER VARIABLES INTEGER DISCF, DISCHK, ENDPTS, I, INEW, 1 IOLD, IP, IXKDIM, J, J1, J1OLD, 2 J2, J2OLD, K, KAIMT, KMAX, KMIN, 3 L, LENDT, NFEVAL, NFJUMP, NSUB, NSUBSV, 4 NXKDIM, PART, SEARCH, TALOC, WHERE, WHERE2 C 11 TO 18 LOGICALS (7 ARE PADDING). LOGICAL DID1, FAIL, FATS(2), FSAVED, HAVDIF, 1 IEND, INIT, ROUNDF, XCDOBT(2), PAD(7) C C THE COMMON BLOCKS. C COMMON /DINTNC/ c 1 2 3 4 5 6 7 8 W AINIT, BINIT, FNCVAL, S, TP, FER, FER1, RELOBT, c 9 10 11 12 13 1 2 3 X TPS, XJ, XJP, FEA, FEA1, KDIM, INC, INC2, c 4 (2,2) 8 9 10 11 12 13 14 Y ISTOP, JPRINT, IPRINT, KK, KMAXF, NDIM, NFINDX, NFMAX, c 15 16 17 18 19 20 Z NFMAXM, RELTOL, REVERM, REVERS, WHEREM, NEEDH COMMON /DINTC/ 1 ACUM, PACUM, RESULT COMMON /DINTC/ c 1 2 (4) 6 7 8 9 10 11 (2) 1 AACUM, LOCAL, ABSCIS, TA, DELTA, DELMIN, DIFF, DISCX, c 13 (2) 15 16 17 (2) 19 20 (24) 44 2 END, ERRINA, ERRINB, FAT, FSAVE, FUNCT, F2, c 45 46 47 48 49 50 51 (6) 3 PAACUM, PF1, PF2, PHISUM, PHTSUM, PX, SPACE, c 57 (2) 59 (2) 61 62 63 64 65 4 STEP, START, SUM, T, TASAVE, TB, TEND, c 66 (2) 68 69 70 71 72 5 WORRY, X1, X2, X, F1, COUNT, c 73 (17) 90 (17) 107 (34) 6 XT, FT, PHI COMMON /DINTC/ c 141 142 143 144 145 146 1 ABSDIF, EDUE2A, EDUE2B, EP, EPNOIZ, EPSMAX, c 147 148 149 150 (2) 152 153 2 EPSO, EPSR, EPSS, ERRAT, ERRC, ERRF, c 154 (2) 156 157 158 159 160 3 ERRT, ESOLD, EXTRA, PEPSMN, RELEPS, REP, c 161 162 163 4 RNDC, TLEN, XJUMP, c 164 165 166 167 168 169 5 ERRI, ERR, EPSMIN, EPS, RE, REPROD COMMON /DINTC/ c 170 171 172 1 DISCF, DISCHK, ENDPTS, INEW, IOLD, IP, IXKDIM, 2 J, J1, J1OLD, J2, J2OLD, KMAX, KMIN, 3 L, LENDT, NFJUMP, NSUBSV, NXKDIM, TALOC, WHERE2, c 1 2 3 4 5 6 7 8 4 I, K, KAIMT, NSUB, PART, SEARCH, WHERE, NFEVAL COMMON /DINTC/ 1 DID1, FAIL, FATS, FSAVED, HAVDIF, IEND, INIT, ROUNDF, 2 XCDOBT, PAD SAVE /DINTNC/, /DINTC/ C C THE VARIABLES HERE DEFINE THE MACHINE ENVIRONMENT. ALL ARE SET C IN DINTOP. THE MEANING ATTACHED TO THESE VARIABLES CAN BE C FOUND BY LOOKING AT THE DEFINITIONS IN DINTOP. DOUBLE PRECISION 1 EMEPS, EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF, 2 ESMALL, ENZER, EDELM1, ENINF COMMON /DINTEC/ 1 EMEPS, EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF, 2 ESMALL, ENZER, EDELM1, ENINF SAVE /DINTEC/ c For initialization of common blocks integer KREAL, KINT, KLOG, ITEMP parameter (KINT=29, KREAL=169, KLOG=11) logical LMOVE(KLOG) integer IMOVE(KINT) double precision RMOVE(KREAL) equivalence (AACUM, RMOVE) equivalence (DISCF, IMOVE) equivalence (DID1, LMOVE) C C C ***** PROCEDURES ****************************************** C c Initialize common blocks to avoid references to undefined variables. do 5 ITEMP = 1, KLOG LMOVE(ITEMP) = .false. 5 continue do 10 ITEMP = 1, KREAL RMOVE(ITEMP) = 0.D0 10 continue do 15 ITEMP = 1, KINT IMOVE(ITEMP) = 0 15 continue WHEREM=0 NFEVAL=0 NDIM=NDIMI KDIM=NWORK X = 0.D0 C KDIM IS TEMPORARILY USED IN DINTMA TO CHECK THE DIMENSION OF C WORK. CALL DINTOP (IOPT,WORK) IF (IOPT(1).EQ.0) THEN C C CALL DINTMA TO DO THE INTEGRATION. C IF (REVERM.EQ.0) CALL DINTMA (ANSWER,WORK,IOPT) ELSE IOPT(1)=-NDIM-4 END IF RETURN END .