1 SUBROUTINE DINTM (NDIMI,ANSWER,WORK,NWORK,IOPT) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5C>> 2009-11-03 DINTM Krogh Initialized lots of variables. 6C>> 1996-03-31 DINTM Krogh Removed unused variable in common. 7c>> 1995-11-20 DINTM Krogh Converted from SFTRAN to Fortran 77. 8C>> 1994-11-14 DINTM Krogh Declared all vars. 9c>> 1994-10-19 DINTM Krogh Changes to use M77CON 10c>> 1994-07-07 DINTM Snyder set up for CHGTYP. 11C>> 1993-05-18 DINTM Krogh -- Changed "END" to "END PROGRAM" 12C>> 1994-05-02 DINTM Snyder corrected some comments 13C>> 1991-09-20 DINTM Krogh converted '(1)' dimensioning to '(*)'. 14C>> 1987-11-19 DINTM Snyder Initial code. 15c 16c--D replaces "?": ?INTA,?intc,?intec,?INTF,?INTM,?INTMA,?INTNC,?INTOP 17C 18C ****************************************************************** 19C 20C THIS SUBROUTINE ATTEMPTS TO CALCULATE THE INTEGRAL OF A FUNCTION 21C PROVIDED BY THE USER. THE FUNCTION IS PROVIDED BY THE USER VIA A 22C SUBROUTINE REFERENCED BY CALL DINTF(F,X,IFLAG) OR BY REVERSE 23C COMMUNICATION. ALL ABSCISSAE ARE WITHIN THE REGION OF 24C INTEGRATION. 25C 26C THE RESULT IS OBTAINED USING QUADRATURE FORMULAE DUE TO 27C T. N. L. PATTERSON, MATHEMATICS OF COMPUTATION, VOLUME 22, 28C PAGES 847-856, 1968. 29C 30C ***** WARNING *********************************************** 31C 32C THE RELIABILITY AND EFFICIENCY OF THIS PROGRAM ARE STRONGLY 33C INFLUENCED BY DISCONTINUITIES IN THE FUNCTION OR IT'S 34C DERIVATIVES, INTEGRABLE SINGULARITIES IN THE REGION OF 35C INTEGRATION, AND NON-INTEGRABLE SINGULARITIES NEAR THIS REGION. 36C (INCLUDING COMPLEX VALUES). THE EFFICIENCY AND RELIABILITY 37C OF INTEGRATING SUCH FUNCTIONS MAY BE GREATLY IMPROVED BY MANUALLY 38C SUBDIVIDING THE REGION OF INTEGRATION AT THE DISCONTINUITY, 39C SINGULARITY OR CLOSEST POINT TO A POLE, AND SUMMING THE ANSWERS. 40C A CHANGE OF VARIABLE TO ELIMINATE OR REDUCE THE STRENGTH OF THE 41C SINGULARITY WILL SIGNIFICANTLY IMPROVE PERFORMANCE. 42C 43C ***** FORMAL ARGUMENTS ************************************* 44C 45C NDIMI IS THE NUMBER OF DIMENSIONS OF INTEGRATION. 46C ANSWER IS THE ESTIMATE OF THE INTEGRAL. WHEN REVERSE COMMUNICATION 47C IS SPECIFIED IT IS USED TO PASS FUNCTION VALUES FROM 48C THE USER PROGRAM TO DINTA OR DINTMA. 49 DOUBLE PRECISION ANSWER 50C WORK CONTAINS THE LIMITS, WORKING SPACE AND MAY BE REFERENCED BY 51C THE OPTION VECTOR (SEE IOPT BELOW). WHEN THE INTEGRATION 52C IS COMPLETE, WORK(1) CONTAINS THE ESTIMATED ABSOLUTE ERROR. 53C WHEN REVERSE COMMUNICATION IS SPECIFIED, WORK(1) IS USED TO 54C PASS ABSCISSAE FROM DINTA OR DINTMA TO THE USER PROGRAM. 55C WORK(1), ..., WORK(NDIMI) ARE COMPONENTS OF THE VECTOR OF 56C ABSCISSAE. WORK(NDIMI+1), ..., WORK(2*NDIMI) ARE LOWER 57C LIMITS. WORK(2*NDIMI+1), ..., WORK(3*NDIMI) ARE UPPER LIMITS. 58C THE ABSCISSAE AND LIMITS ARE STORED INNERMOST FIRST. DINTF 59C IS CALLED TO REQUEST THE LIMITS OF EVERY DIMENSION, BUT 60C CONSTANT LIMITS (THE LIMITS OF THE OUTER DIMENSION ARE ALWAYS 61C CONSTANT) MAY INSTEAD BE STORED BEFORE INTEGRATION BEGINS. 62C WORK(3*NDIMI+1), ..., WORK(3*NDIMI+KWORK*(NDIMI-1)) ARE 63C WORKING STORAGE, WHERE KWORK DEPENDS ON THE MACHINE AND THE 64C PRECISION OF THE PROGRAM VERSION. KWORK MUST SPECIFY ENOUGH 65C STORAGE FOR 4 VARIABLES WHICH ARE ALWAYS DOUBLE PRECISION, 66C 139 VARIABLES WHICH ARE DOUBLE PRECISION IN THE DOUBLE 67C PRECISION VERSION, 30 VARIABLES WHICH ARE DOUBLE PRECISION 68C IN THE DOUBLE PRECISION PROGRAM IF THE LARGEST REPRESENTABLE 69C NUMBERS IN SINGLE AND DOUBLE PRECISION ARE DIFFERENT BY MORE 70C THAN A FACTOR OF 4, 29 INTEGER VARIABLES AND 11 LOGICAL 71C VARIABLES. REMEMBER TO CONSIDER THE DATA TYPE OF WORK WHEN 72C SPECIFYING THE VALUE OF KWORK. EXAMPLES ARE 73C-- Begin mask code changes 74C IF PRECISION IS REAL, AND INTEGER AND LOGICAL VARIABLES 75C OCCUPY THE SAME SPACE AS REAL VARIABLES, KWORK=217. 76C IF PRECISION IS REAL, AND INTEGER AND LOGICAL VARIABLES 77C OCCUPY HALF AS MUCH SPACE AS REAL VARIABLES, KWORK=? 78C IF PRECISION IS PARTLY DOUBLE, AND INTEGER AND LOGICAL 79C VARIABLES OCCUPY THE SAME SPACE AS REAL VARIABLES, 80C KWORK=? 81C IF PRECISION IS PARTLY DOUBLE, AND INTEGER AND LOGICAL 82C VARIABLES OCCUPY HALF AS MUCH SPACE AS REAL VARIABLES, 83C KWORK=? 84C IF PRECISION IS DOUBLE, AND INTEGER AND LOGICAL VARIABLES 85C OCCUPY THE SAME SPACE AS REAL VARIABLES, KWORK=193. 86C IF PRECISION IS DOUBLE, AND INTEGER AND LOGICAL VARIABLES 87C OCCUPY HALF AS MUCH SPACE AS REAL VARIABLES, KWORK=? 88C FOR PORTABILITY, ASSUME THE WORST CASE. THAT IS, IF THE 89C PRECISION IS SINGLE, KWORK=217. IF THE PRECISION IS DOUBLE, 90C KWORK=193. 91C WORK(3*NDIMI+KWORK*(NDIMI-1)+1), ... ARE AVAILABLE FOR 92C REFERENCE THROUGH THE OPTION VECTOR. 93C-- End mask code changes 94 DOUBLE PRECISION WORK(*) 95C NWORK IS THE NUMBER OF ELEMENTS ALLOCATED FOR WORK. 96C NWORK MUST BE GREATER THAN 3*NDIMI+KWORK*(NDIMI-1). 97C IOPT IS A VECTOR OF INTEGERS USED TO RETURN THE STATUS AND TO 98C SPECIFY OPTIONS. 99 INTEGER IOPT(*) 100C IOPT(1) IS USED TO RETURN A STATUS INDICATOR TO THE USER. 101C VALUES OF THIS INDICATOR ARE 102C -NDIM - NORMAL TERMINATION, EITHER THE ABSOLUTE OR THE 103C RELATIVE ERROR CRITERIA ARE SATISFIED. 104C -NDIM-1 - NORMAL TERMINATION, NEITHER THE ABSOLUTE NOR THE 105C RELATIVE ERROR CRITERIA ARE SATISFIED, BUT THE 106C ERROR RELATIVE TO THE OBTAINABLE PRECISION 107C PRECISION CRITERION IS SATISFIED. 108C -NDIM-2 - NORMAL TERMINATION, BUT NONE OF THE ERROR CRITERIA 109C ARE SATISFIED. 110C -NDIM-3 - NOT ENOUGH SPACE IN WORK, NWORK IS TOO SMALL. 111C -NDIM-4 - BAD IOPT VALUE. 112C -NDIM-5 - TOO MANY FUNCTION VALUES NEEDED. 113C -NDIM-KDIM-5 - APPARENT NON-INTEGRABLE SINGULARITY IN 114C DIMENSION KDIM. ANSWER CONTAINS THE APPROXIMATE 115C ABSCISSA OF THE SINGULARITY IN THE KDIM-TH 116C DIMENSION, WORK(KDIM+1) THROUGH WORK(NDIM) 117C CONTAIN THE ABSCISSAE OF EXTERIOR INTEGRALS. 118C -2*NDIM-KDIM-5 - INCORRECT INNER INTEGRAL DIMENSIONALITY 119C SPECIFIED DURING EVALUATION OF INTEGRAL AT 120C DIMENSION KDIM, THAT IS, IABS(IFLAG(1)) = KDIM. 121C ENTRIES IN IOPT STARTING WITH IOPT(2) ARE DESCRIBED BELOW. 122C IOPT(I) ENTRY IN IOPT(I) MEANS 123C 0 NO MORE OPTIONS, BEGIN INTEGRATION. 124C 1 IOPT(I+1) CONTAINS THE UNIT NUMBER FOR OUTPUT. 125C 2 IOPT(I+1) IS AN NDIMI DIGIT INTEGER, WHERE 126C EACH DIGIT IS THE DIAGNOSTIC PRINT LEVEL FOR 127C ONE DIMENSION. THE LOW ORDER DIGIT IS THE PRINT 128C LEVEL FOR THE INNER DIMENSION. 129C 0 - NO PRINTING 130C 1 - MINIMUM PRINTING - ERROR MESSAGES (DEFAULT) 131C 2 - PANEL BOUNDARIES AND ANSWERS 132C 3 - ERROR ESTIMATES FOR EACH QUADRATURE FORMULA 133C 4 - DETAILED OUTPUT (DIFFERENCE LINES, ETC). 134C 3 WORK(IOPT(I+1)) CONTAINS THE ABSOLUTE TOLERANCE, 135C WORK(IOPT(I+1)+1) CONTAINS TOLERANCE RELATIVE TO 136C THE LOCALLY OBTAINABLE PRECISION, AND 137C WORK(IOPT(I+1)+2) CONTAINS TOLERANCE RELATIVE TO 138C THE VALUE OF THE INTEGRAL. THE TOLERANCE RELATIVE 139C TO THE LOCALLY OBTAINABLE PRECISION IS SPECIFIED AS 140C THE FRACTION OF LOCALLY OBTAINABLE DIGITS THAT ARE 141C PERMITTED TO BE WRONG, AND IS INTERNALLY BOUNDED 142C BETWEEN 0.0 AND 1.0. IF ANY OF THE ERROR CONTROL 143C CRITERIA ARE SATISFIED, THE ANSWER IS ACCEPTED. IF 144C THIS OPTION IS NOT USED, THE ABSOLUTE AND RELATIVE 145C TOLERANCES ARE SET TO ZERO, AND THE TOLERANCE 146C RELATIVE TO THE LOCALLY OBTAINABLE PRECISION IS SET 147C TO 0.25. 148C 4 WORK(IOPT(I+1)) CONTAINS ABSOLUTE ERROR COMMITTED 149C COMPUTING F. WORK(IOPT(I+1)) MAY BE CHANGED DURING 150C THE INTEGRATION. 151C 5 WORK(IOPT(I+1)) CONTAINS RELATIVE ERROR EXPECTED TO 152C BE COMMITTED COMPUTING F. CHANGES TO THIS VALUE 153C DURING INTEGRATION WILL NOT BE DETECTED. 154C 6 USE REVERSE COMMUNICATION - 155C CALL DINTM (NDIMI,ANSWER,WORK,NWORK,IOPT) 156C 1 CALL DINTMA (ANSWER,WORK,IOPT(1)) 157C IF (IOPT(1)) 3,2,4 158C 2 COMPUTE THE INNERMOST INTEGRAND - 159C ANSWER=F(WORK(K)), K=1,NDIMI 160C THEN EITHER 161C GO TO 1 162C OR (FOR BETTER EFFICIENCY) 163C CALL DINTA (ANSWER,WORK,IOPT(1)) 164C IF (IOPT(1)) 1,2,5 165C VALUES OF IOPT(1) PRODUCED BY DINTA AND DINTMA 166C HAVE DIFFERENT MEANINGS. 167C 3 IF (IOPT(1)+NDIMI.LE.0) GO TO 6 168C THE INTEGRAL OVER THE (-IOPT(1))TH DIMENSION IS 169C CONTAINED IN ANSWER. IF A TRANSFORMATION OF 170C THE INTEGRAL IS NOT NECESSARY BEFORE IT IS USED 171C AS AN INTEGRAND FOR THE (1-IOPT(1))TH DIMENSION 172C GO TO 1 173C ELSE COMPUTE ANSWER=G(ANSWER,WORK(K)), 174C K=1-IOPT(1),NDIMI AND G IS THE TRANSFORMATION. 175C COMPUTE ALSO WORK(1)=PARTIAL DERIVATIVE OF G 176C WITH RESPECT TO ANSWER THEN 177C GO TO 1 178C 4 COMPUTE WORK(NDIMI+IOPT(1))=A(WORK(K)), 179C K=IOPT(1)+1,NDIMI, WORK(2*NDIMI+IOPT(1))=B(WORK(K)) 180C K=IOPT(1)+1,NDIM, WHERE A IS THE LOWER LIMIT OF 181C INTEGRATION FOR THE (IOPT(1))TH DIMENSION, AND B IS 182C THE UPPER LIMIT. IF A TRANSFORMATION AS DISCUSSED 183C ABOVE IS TO BE APPLIED WHEN THE INTEGRATION OF THIS 184C DIMENSION IS COMPLETE, COMPUTE WORK(1)=AN ESTIMATE 185C OF THE UPPER BOUND OF THE PARTIAL DERIVATIVE OF G 186C WITH RESPECT TO THE INTEGRAL. IF G IS LINEAR IN 187C THE INTEGRAL, THIS DERIVATIVE CAN BE COMPUTED 188C WITHOUT KNOWING THE INTEGRAL, RATHER THAN 189C ESTIMATED. IF THERE IS A SINGULARITY IN THE 190C INTEGRAND FOR THIS DIMENSION (PERHAPS INTRODUCED 191C BY THE BOUNDARY OR THE TRANSFORMATION G), SET 192C IOPT(1) TO THE LOCATION IN WORK OF THE ABSCISSA IN 193C THIS DIMENSION OF THE SINGULARITY, AND WITH THE 194C SIGN SET AS DESCRIBED FOR OPTION 11 BELOW. THE 195C MAGNITUDE OF IOPT(1) MUST BE GREATER THAN NDIM FOR 196C THIS REQUEST TO BE DETECTED. IN PARTICULAR, IF 197C THERE IS A SINGULARITY AT ONE OF THE LIMITS, SET 198C IOPT(1) TO IOPT(1)+NDIM OR -(IOPT(1)+2*NDIM). THEN 199C GO TO 1. 200C 5 IOPT(1)=-(IOPT(1)+NDIM) 201C 6 CONTINUE 202C AT THIS POINT, IERR CONTAINS THE STATUS DESCRIBED 203C BELOW, FOR OPTION 8. 204C 7 SET MINIMUM INDEX OF QUADRATURE FORMULA TO USE 205C BEFORE SUBDIVISION TO IOPT(I+1). 206C 8 NO EFFECT. IOPT(I+1) MAY BE USED TO PASS 207C INFORMATION TO DINTF. INCREMENT I BY 2. 208C 9 IOPT(I+1) IS THE MAXIMUM NUMBER OF FUNCTION VALUES 209C TO USE TO INTEGRATE THE ENTIRE PROBLEM. IF 210C IOPT(I+1) .LE. 0, THE NUMBER OF FUNCTION VALUES 211C IS NOT CONTROLLED. 212C 10 IOPT(I+1) IS USED TO RETURN THE NUMBER OF FUNCTION 213C VALUES USED TO INTEGRATE THE ENTIRE PROBLEM. 214C 11 WORK(IABS(IOPT(I+1))) IS THE LOCATION OF A 215C SINGULARITY OR DISCONTINUITY IN THE OUTERMOST 216C INTEGRAL. IF THE LOCATION IS INSIDE THE INTERVAL, 217C THE INTERVAL IS IMMEDIATELY SUBDIVIDED. IF 218C IOPT(I+1) .GT. 0, A T**2 SUBSTITUTION BASED AT 219C WORK(...) WILL BE USED. IF IOPT(I+1) .LT. 0, 220C A T**4 SUBSTITUTION WILL BE USED. IF IOPT(I+1) 221C .EQ. 0, NO SUBSTITUTION WILL BE USED. TO NOTIFY 222C THE PROGRAM OF SINGULARITIES OR DISCONTINUITIES IN 223C INTERIOR INTEGRANDS, SET IFLAG NEGATIVE WHEN 224C ASKED TO COMPUTE LIMITS. SEE THE DESCRIPTION 225C OF REVERSE COMMUNICATION ABOVE, OR THE DESCRIPTION 226C OF DINTF BELOW. 227C 12 WORK(IOPT(I+1)) CONTAINS THE ABSOLUTE ERROR IN THE 228C LOWER LIMIT OF THE CURRENT DIMENSION. THE ERROR IN 229C THE UPPER LIMIT IS IN WORK(IOPT(I+1)+1). 230C 13 IOPT(I+1) IS THE LOCATION IN IOPT IN WHICH THE USER 231C IS TO STORE NON-STANDARD CHANGES TO THE DIMENSION OF 232C INNER INTEGRALS. THE DEFAULT VALUE IS 1, ALLOWING 233C SUCH NONSTANDARD CHANGES TO BE STORED IN THE CELL OF 234C IFLAG (=IOPT) USUALLY USED FOR COMMUNICATION WITH 235C DINTF, IN THE CASE WHEN IT IS NOT SIMULTANEOUSLY 236C NECESSARY TO INDICATE A SINGULARITY. SEE THE 237C DESCRIPTION OF THE INTERFACE TO DINTF BELOW. 238C 239C ALL OPTIONS ARE SET TO THEIR NOMINAL VALUES BEFORE THE OPTION 240C VECTOR IS PROCESSED. 241C 242C ***** DINTF *********************************************** 243C 244C DINTF IS REFERENCED VIA CALL DINTF (ANSWER,WORK,IFLAG). 245C WORK AND ANSWER ARE AS DESCRIBED ABOVE. VALUES OF IFLAG ARE 246C IFLAG .LT. 0 247C THE INTEGRAL OVER THE (-IFLAG)TH DIMENSION IS 248C CONTAINED IN ANSWER. IF A TRANSFORMATION OF 249C THE INTEGRAL IS NOT NECESSARY BEFORE IT IS USED 250C AS AN INTEGRAND FOR THE (1-IFLAG)TH DIMENSION 251C RETURN 252C ELSE IF THE INTEGRAND DEPENDS ON ONLY ONE INNER 253C INTEGRAL COMPUTE ANSWER=G(ANSWER,WORK(K)), 254C ALSO MULTIPLY WORK(1) BY THE PARTIAL DERIVATIVE 255C OF G WITH RESPECT TO ANSWER THEN 256C RETURN 257C ELSE IF THE INTEGRAND DEPENDS ON MORE INTEGRALS 258C SAVE ANSWER AND WORK(1), AND SET IFLAG(IXKDIM) TO 259C THE DIMENSION OF INTEGRAL NEXT TO BE COMPUTED 260C (IXKDIM IS THE VALUE SPECIFIED BY OPTION 12, OR 1 261C IF OPTION 12 IS NOT SELECTED), THEN 262C RETURN 263C ELSE (ALL INNER INTEGRALS UPON WHICH THE INNER 264C INTEGRAND DEPENDS HAVE BEEN EVALUATED) CALCULATE 265C ANSWER=G(ANSWER,WORK(K),SAVED INTEGRALS) 266C K=1-IFLAG, NDIMI (G IS THE TRANSFORMATION), 267C MULTIPLY WORK(1) BY THE PARTIAL DERIVATIVE OF G 268C WITH RESPECT TO THE LAST INTEGRAL, AND ADD THE 269C PRODUCT OF EACH OF THE SAVED VALUES OF WORK(1) AND 270C PARTIAL DERIVATIVE OF G WITH RESPECT TO THE 271C CORRESPONDING SAVED INTEGRAL ONTO WORK(1). THAT 272C IS, WORK(1) IS THE TOTAL ERROR IN G, CALCULATED AS 273C THE INNER PRODUCT OF THE ERRORS IN THE INDIVIDUAL 274C ARGUMENTS, AND THE PARTIAL DERIVATIVES OF G WITH 275C RESPECT TO THOSE ARGUMENTS. 276C IFLAG .EQ. 0 277C COMPUTE THE INNERMOST INTEGRAND - 278C ANSWER=F(WORK(K)), K=1,NDIMI 279C AND 280C RETURN 281C IFLAG .GT. 0 282C COMPUTE WORK(NDIMI+IFLAG)=A(WORK(K)), 283C K=IFLAG+1,NDIMI, WORK(2*NDIMI+IFLAG)=B(WORK(K)), 284C K=IFLAG+1,NDIM, WHERE A IS THE LOWER LIMIT OF 285C INTEGRATION FOR THE (IFLAG)TH DIMENSION, AND B IS 286C THE UPPER LIMIT. IF A TRANSFORMATION AS DISCUSSED 287C ABOVE IS TO BE APPLIED WHEN THE INTEGRATION OF THIS 288C DIMENSION IS COMPLETE, COMPUTE WORK(1)=AN ESTIMATE 289C OF THE UPPER BOUND OF THE PARTIAL DERIVATIVE 290C OF G WITH RESPECT TO THE INTEGRAL. IF G IS 291C LINEAR IN THE INTEGRAL, THIS DERIVATIVE MAY BE 292C COMPUTED WITHOUT KNOWING THE INTEGRAL, RATHER 293C THAN ESTIMATED. IF THERE IS A SINGULARITY IN 294C THE INTEGRAND FOR THIS DIMENSION (PERHAPS 295C INTRODUCED BY THE BOUNDARY OR THE TRANSFORMATION), 296C SET IFLAG TO THE LOCATION IN WORK OF THE ABSCISSA 297C IN THIS DIMENSION OF THE SINGULARITY, AND SET THE 298C SIGN AS FOR OPTION 11 ABOVE. NOTE THAT THE 299C MAGNITUDE OF IFLAG MUST BE GREATER THAN NDIM FOR 300C THIS REQUEST TO BE DETECTED. IN PARTICULAR, IF 301C THERE IS A SINGULARITY AT ONE OF THE LIMITS, SET 302C THE MAGNITUDE OF IFLAG TO IFLAG+NDIM OR 303C IFLAG+2*NDIM, WITH THE SIGN AS DESCRIBED FOR 304C OPTION 11. THEN 305C RETURN. 306C IF THE FIRST INTEGRAL TO BE COMPUTED DOES NOT HAVE 307C DIMENSION IFLAG(1)-1, SET IFLAG(IXKDIM) TO THE 308C DIMENSION OF INTEGRAL TO BE NEXT COMPUTED. SEE 309C OPTION 12 AND THE DISCUSSION OF IFLAG .LT. 0 ABOVE 310C FOR INFORMATION ON IXKDIM. 311C 312C ***** EXTERNAL REFERENCES ********************************* 313C 314C DINTMA TO DO THE INTEGRATION. 315C DINTOP TO SET OPTIONS. 316C 317C ***** LOCAL VARIABLES ************************************* 318C 319 INTEGER NDIMI, NWORK 320C 321C ***** COMMON STORAGE ************************************** 322C 323C COMMON /DINTNC/ CONTAINS VARIABLES NOT SEPARATELY SAVED FOR 324C EACH DIMENSION OF A MULTIPLE QUADRATURE. COMMON /DINTC/ 325C CONTAINS VARIABLES THAT MUST BE SAVED FOR EACH DIMENSION OF THE 326C QUADRATURE. THE VARIABLES IN EACH COMMON BLOCK ARE STORED IN THE 327C ORDER - ALWAYS DOUBLE, DOUBLE IF DOUBLE PRECISION PROGRAM, DOUBLE 328C IF DOUBLE PRECISION PROGRAM AND EXPONENT RANGE OF DOUBLE AND 329C SINGLE VERY DIFFERENT, SINGLE, INTEGER, LOGICAL. A PAD OF LOGICAL 330C VARIABLES IS INCLUDED AT THE END OF /DINTC/. THE DIMENSION OF 331C THE PAD MAY NEED TO BE VARIED SO THAT NO VARIABLES BEYOND THE END 332C OF THE COMMON BLOCK ARE ALTERED. 333C 334C DECLARATIONS OF COMMON /DINTNC/ VARIABLES. 335C 336 DOUBLE PRECISION AINIT, BINIT, FNCVAL, S, TP 337 DOUBLE PRECISION FER, FER1, RELOBT, TPS, XJ, XJP 338 INTEGER FEA, FEA1, INC, INC2, IPRINT, 339 1 ISTOP(2,2),JPRINT, KDIM, KK, KMAXF, NDIM, 340 2 NFINDX, NFMAX, NFMAXM, RELTOL, REVERM, REVERS, 341 3 WHEREM 342 LOGICAL NEEDH 343C 344C DECLARATIONS OF COMMON /DINTC/ VARIABLES. 345C 346c--D Next line special: S => D, X => Q, D => D, P => D 347 DOUBLE PRECISION ACUM, PACUM, RESULT(2) 348C 139 $.TYPE.$ VARIABLES 349 DOUBLE PRECISION 350 1 AACUM, ABSCIS, DELMIN, DELTA, DIFF, DISCX(2), 351 2 END(2), ERRINA, ERRINB, FAT(2), FSAVE, 352 3 FUNCT(24), F1, F2, LOCAL(4), PAACUM, PF1, 353 4 PF2, PHISUM, PHTSUM, PX, SPACE(6), 354 5 STEP(2), START(2), SUM, T, TA, TASAVE, 355 6 TB, TEND, WORRY(2), X, X1, 356 7 X2, XT(17), FT(17), PHI(34) 357c Note XT, FT, and PHI above are last, because they must be in adjacent 358c locations in DINTC. 359C 30 $DSTYP$ VARIABLES 360 DOUBLE PRECISION 361 1 ABSDIF, COUNT, EDUE2A, EDUE2B, EP, EPNOIZ, 362 2 EPS, EPSMAX, EPSMIN, EPSO, EPSR, EPSS, 363 3 ERR, ERRAT(2), ERRC, ERRF, ERRI, ERRT(2), 364 4 ESOLD, EXTRA, PEPSMN, RE, RELEPS, REP, 365 5 REPROD, RNDC, TLEN, XJUMP 366C 29 INTEGER VARIABLES 367 INTEGER DISCF, DISCHK, ENDPTS, I, INEW, 368 1 IOLD, IP, IXKDIM, J, J1, J1OLD, 369 2 J2, J2OLD, K, KAIMT, KMAX, KMIN, 370 3 L, LENDT, NFEVAL, NFJUMP, NSUB, NSUBSV, 371 4 NXKDIM, PART, SEARCH, TALOC, WHERE, WHERE2 372C 11 TO 18 LOGICALS (7 ARE PADDING). 373 LOGICAL DID1, FAIL, FATS(2), FSAVED, HAVDIF, 374 1 IEND, INIT, ROUNDF, XCDOBT(2), PAD(7) 375C 376C THE COMMON BLOCKS. 377C 378 COMMON /DINTNC/ 379c 1 2 3 4 5 6 7 8 380 W AINIT, BINIT, FNCVAL, S, TP, FER, FER1, RELOBT, 381c 9 10 11 12 13 1 2 3 382 X TPS, XJ, XJP, FEA, FEA1, KDIM, INC, INC2, 383c 4 (2,2) 8 9 10 11 12 13 14 384 Y ISTOP, JPRINT, IPRINT, KK, KMAXF, NDIM, NFINDX, NFMAX, 385c 15 16 17 18 19 20 386 Z NFMAXM, RELTOL, REVERM, REVERS, WHEREM, NEEDH 387 COMMON /DINTC/ 388 1 ACUM, PACUM, RESULT 389 COMMON /DINTC/ 390c 1 2 (4) 6 7 8 9 10 11 (2) 391 1 AACUM, LOCAL, ABSCIS, TA, DELTA, DELMIN, DIFF, DISCX, 392c 13 (2) 15 16 17 (2) 19 20 (24) 44 393 2 END, ERRINA, ERRINB, FAT, FSAVE, FUNCT, F2, 394c 45 46 47 48 49 50 51 (6) 395 3 PAACUM, PF1, PF2, PHISUM, PHTSUM, PX, SPACE, 396c 57 (2) 59 (2) 61 62 63 64 65 397 4 STEP, START, SUM, T, TASAVE, TB, TEND, 398c 66 (2) 68 69 70 71 72 399 5 WORRY, X1, X2, X, F1, COUNT, 400c 73 (17) 90 (17) 107 (34) 401 6 XT, FT, PHI 402 COMMON /DINTC/ 403c 141 142 143 144 145 146 404 1 ABSDIF, EDUE2A, EDUE2B, EP, EPNOIZ, EPSMAX, 405c 147 148 149 150 (2) 152 153 406 2 EPSO, EPSR, EPSS, ERRAT, ERRC, ERRF, 407c 154 (2) 156 157 158 159 160 408 3 ERRT, ESOLD, EXTRA, PEPSMN, RELEPS, REP, 409c 161 162 163 410 4 RNDC, TLEN, XJUMP, 411c 164 165 166 167 168 169 412 5 ERRI, ERR, EPSMIN, EPS, RE, REPROD 413 COMMON /DINTC/ 414c 170 171 172 415 1 DISCF, DISCHK, ENDPTS, INEW, IOLD, IP, IXKDIM, 416 2 J, J1, J1OLD, J2, J2OLD, KMAX, KMIN, 417 3 L, LENDT, NFJUMP, NSUBSV, NXKDIM, TALOC, WHERE2, 418c 1 2 3 4 5 6 7 8 419 4 I, K, KAIMT, NSUB, PART, SEARCH, WHERE, NFEVAL 420 COMMON /DINTC/ 421 1 DID1, FAIL, FATS, FSAVED, HAVDIF, IEND, INIT, ROUNDF, 422 2 XCDOBT, PAD 423 SAVE /DINTNC/, /DINTC/ 424C 425C THE VARIABLES HERE DEFINE THE MACHINE ENVIRONMENT. ALL ARE SET 426C IN DINTOP. THE MEANING ATTACHED TO THESE VARIABLES CAN BE 427C FOUND BY LOOKING AT THE DEFINITIONS IN DINTOP. 428 DOUBLE PRECISION 429 1 EMEPS, EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF, 430 2 ESMALL, ENZER, EDELM1, ENINF 431 COMMON /DINTEC/ 432 1 EMEPS, EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF, 433 2 ESMALL, ENZER, EDELM1, ENINF 434 SAVE /DINTEC/ 435 436c For initialization of common blocks 437 integer KREAL, KINT, KLOG, ITEMP 438 parameter (KINT=29, KREAL=169, KLOG=11) 439 logical LMOVE(KLOG) 440 integer IMOVE(KINT) 441 double precision RMOVE(KREAL) 442 equivalence (AACUM, RMOVE) 443 equivalence (DISCF, IMOVE) 444 equivalence (DID1, LMOVE) 445C 446C 447C ***** PROCEDURES ****************************************** 448C 449 450c Initialize common blocks to avoid references to undefined variables. 451 452 do 5 ITEMP = 1, KLOG 453 LMOVE(ITEMP) = .false. 454 5 continue 455 do 10 ITEMP = 1, KREAL 456 RMOVE(ITEMP) = 0.D0 457 10 continue 458 do 15 ITEMP = 1, KINT 459 IMOVE(ITEMP) = 0 460 15 continue 461 462 WHEREM=0 463 NFEVAL=0 464 NDIM=NDIMI 465 KDIM=NWORK 466 X = 0.D0 467 468 469C KDIM IS TEMPORARILY USED IN DINTMA TO CHECK THE DIMENSION OF 470C WORK. 471 CALL DINTOP (IOPT,WORK) 472 IF (IOPT(1).EQ.0) THEN 473C 474C CALL DINTMA TO DO THE INTEGRATION. 475C 476 IF (REVERM.EQ.0) CALL DINTMA (ANSWER,WORK,IOPT) 477 ELSE 478 IOPT(1)=-NDIM-4 479 END IF 480 RETURN 481 END 482