1 SUBROUTINE DINTDU 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5C>> 2007-03-28 DINTDU Snyder Don't look at XT(0) or FT(0) 6C>> 2007-03-28 DINTDU Krogh l .le. 0 changed to l .le. 1 7C>> 1996-03-31 DINTDU Krogh Removed unused variable in common. 8c>> 1995-11-20 DINTDU Krogh Converted from SFTRAN to Fortran 77. 9c>> 1994-10-19 DINTDU Krogh Changes to use M77CON 10c>> 1994-08-19 DINTDU Snyder correct "middle" that's really at alocal 11c>> 1994-07-07 DINTDU Snyder set up for CHGTYP. 12C>> 1993-05-18 DINTDU Krogh -- Changed "END" to "END PROGRAM" 13C>> 1987-11-20 DINTDU Snyder Initial code. 14C 15C THIS SUBROUTINE UPDATES DIFFERENCE LINES FOR DINTA DURING 16C THE SEARCHES. 17c 18c--D replaces "?": ?INTA, ?intc, ?INTDU, ?intec, ?intnc 19C 20C ***** INTERNAL AND COMMON VARIABLES ************************ 21C 22C EPSCOR IS A CORRECTION TO BE ADDED ONTO EPSMIN. 23 DOUBLE PRECISION EPSCOR 24C FATA THE FUNCTION VALUE AT THE ALOCAL END OF THE INTERVAL. 25C FATB THE FUNCTION VALUE AT THE BLOCAL END OF THE INTERVAL. 26 DOUBLE PRECISION FATA,FATB 27C PHIT IS THE BACKWARD DIFFERENCE LINE. 28 DOUBLE PRECISION PHIT(17) 29C 30C ***** COMMON STORAGE ****************************************** 31C 32C COMMON /DINTNC/ CONTAINS VARIABLES NOT SEPARATELY SAVED FOR 33C EACH DIMENSION OF A MULTIPLE QUADRATURE. COMMON /DINTC/ 34C CONTAINS VARIABLES THAT MUST BE SAVED FOR EACH DIMENSION OF THE 35C QUADRATURE. THE VARIABLES IN EACH COMMON BLOCK ARE STORED IN THE 36C ORDER - ALWAYS DOUBLE, DOUBLE IF DOUBLE PRECISION PROGRAM, DOUBLE 37C IF DOUBLE PRECISION PROGRAM AND EXPONENT RANGE OF DOUBLE AND 38C SINGLE VERY DIFFERENT, SINGLE, INTEGER, LOGICAL. A PAD OF LOGICAL 39C VARIABLES IS INCLUDED AT THE END OF /DINTC/. THE DIMENSION OF 40C THE PAD MAY NEED TO BE VARIED SO THAT NO VARIABLES BEYOND THE END 41C OF THE COMMON BLOCK ARE ALTERED. 42C 43C DECLARATIONS OF COMMON /DINTNC/ VARIABLES. 44C 45 DOUBLE PRECISION AINIT, BINIT, FNCVAL, S, TP 46 DOUBLE PRECISION FER, FER1, RELOBT, TPS, XJ, XJP 47 INTEGER FEA, FEA1, INC, INC2, IPRINT, 48 1 ISTOP(2,2),JPRINT, KDIM, KK, KMAXF, NDIM, 49 2 NFINDX, NFMAX, NFMAXM, RELTOL, REVERM, REVERS, 50 3 WHEREM 51 LOGICAL NEEDH 52C 53C DECLARATIONS OF COMMON /DINTC/ VARIABLES. 54C 55c--D Next line special: S => D, X => Q, D => D, P => D 56 DOUBLE PRECISION ACUM, PACUM, RESULT(2) 57C 139 $.TYPE.$ VARIABLES 58 DOUBLE PRECISION 59 1 AACUM, ABSCIS, DELMIN, DELTA, DIFF, DISCX(2), 60 2 END(2), ERRINA, ERRINB, FAT(2), FSAVE, 61 3 FUNCT(24), F1, F2, LOCAL(4), PAACUM, PF1, 62 4 PF2, PHISUM, PHTSUM, PX, SPACE(6), 63 5 STEP(2), START(2), SUM, T, TA, TASAVE, 64 6 TB, TEND, WORRY(2), X, X1, 65 7 X2, XT(17), FT(17), PHI(34) 66c Note XT, FT, and PHI above are last, because they must be in adjacent 67c locations in DINTC. 68C 30 $DSTYP$ VARIABLES 69 DOUBLE PRECISION 70 1 ABSDIF, COUNT, EDUE2A, EDUE2B, EP, EPNOIZ, 71 2 EPS, EPSMAX, EPSMIN, EPSO, EPSR, EPSS, 72 3 ERR, ERRAT(2), ERRC, ERRF, ERRI, ERRT(2), 73 4 ESOLD, EXTRA, PEPSMN, RE, RELEPS, REP, 74 5 REPROD, RNDC, TLEN, XJUMP 75C 29 INTEGER VARIABLES 76 INTEGER DISCF, DISCHK, ENDPTS, I, INEW, 77 1 IOLD, IP, IXKDIM, J, J1, J1OLD, 78 2 J2, J2OLD, K, KAIMT, KMAX, KMIN, 79 3 L, LENDT, NFEVAL, NFJUMP, NSUB, NSUBSV, 80 4 NXKDIM, PART, SEARCH, TALOC, WHERE, WHERE2 81C 11 TO 18 LOGICALS (7 ARE PADDING). 82 LOGICAL DID1, FAIL, FATS(2), FSAVED, HAVDIF, 83 1 IEND, INIT, ROUNDF, XCDOBT(2), PAD(7) 84C 85C THE COMMON BLOCKS. 86C 87 COMMON /DINTNC/ 88c 1 2 3 4 5 6 7 8 89 W AINIT, BINIT, FNCVAL, S, TP, FER, FER1, RELOBT, 90c 9 10 11 12 13 1 2 3 91 X TPS, XJ, XJP, FEA, FEA1, KDIM, INC, INC2, 92c 4 (2,2) 8 9 10 11 12 13 14 93 Y ISTOP, JPRINT, IPRINT, KK, KMAXF, NDIM, NFINDX, NFMAX, 94c 15 16 17 18 19 20 95 Z NFMAXM, RELTOL, REVERM, REVERS, WHEREM, NEEDH 96 COMMON /DINTC/ 97 1 ACUM, PACUM, RESULT 98 COMMON /DINTC/ 99c 1 2 (4) 6 7 8 9 10 11 (2) 100 1 AACUM, LOCAL, ABSCIS, TA, DELTA, DELMIN, DIFF, DISCX, 101c 13 (2) 15 16 17 (2) 19 20 (24) 44 102 2 END, ERRINA, ERRINB, FAT, FSAVE, FUNCT, F2, 103c 45 46 47 48 49 50 51 (6) 104 3 PAACUM, PF1, PF2, PHISUM, PHTSUM, PX, SPACE, 105c 57 (2) 59 (2) 61 62 63 64 65 106 4 STEP, START, SUM, T, TASAVE, TB, TEND, 107c 66 (2) 68 69 70 71 72 108 5 WORRY, X1, X2, X, F1, COUNT, 109c 73 (17) 90 (17) 107 (34) 110 6 XT, FT, PHI 111 COMMON /DINTC/ 112c 141 142 143 144 145 146 113 1 ABSDIF, EDUE2A, EDUE2B, EP, EPNOIZ, EPSMAX, 114c 147 148 149 150 (2) 152 153 115 2 EPSO, EPSR, EPSS, ERRAT, ERRC, ERRF, 116c 154 (2) 156 157 158 159 160 117 3 ERRT, ESOLD, EXTRA, PEPSMN, RELEPS, REP, 118c 161 162 163 119 4 RNDC, TLEN, XJUMP, 120c 164 165 166 167 168 169 121 5 ERRI, ERR, EPSMIN, EPS, RE, REPROD 122 COMMON /DINTC/ 123c 170 171 172 124 1 DISCF, DISCHK, ENDPTS, INEW, IOLD, IP, IXKDIM, 125 2 J, J1, J1OLD, J2, J2OLD, KMAX, KMIN, 126 3 L, LENDT, NFJUMP, NSUBSV, NXKDIM, TALOC, WHERE2, 127c 1 2 3 4 5 6 7 8 128 4 I, K, KAIMT, NSUB, PART, SEARCH, WHERE, NFEVAL 129 COMMON /DINTC/ 130 1 DID1, FAIL, FATS, FSAVED, HAVDIF, IEND, INIT, ROUNDF, 131 2 XCDOBT, PAD 132 SAVE /DINTNC/, /DINTC/ 133C 134C THE VARIABLES HERE DEFINE THE MACHINE ENVIRONMENT. ALL ARE SET 135C IN DINTOP. THE MEANING ATTACHED TO THESE VARIABLES CAN BE 136C FOUND BY LOOKING AT THE DEFINITIONS IN DINTOP. 137 DOUBLE PRECISION 138 1 EMEPS, EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF, 139 2 ESMALL, ENZER, EDELM1, ENINF 140 COMMON /DINTEC/ 141 1 EMEPS, EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF, 142 2 ESMALL, ENZER, EDELM1, ENINF 143 SAVE /DINTEC/ 144C 145C ***** EQUIVALENCE STATEMENTS ******************************* 146C 147 EQUIVALENCE (PHI(18),PHIT) 148 EQUIVALENCE (FAT(1),FATA), (FAT(2),FATB) 149C 150C ***** PROCEDURES ****************************************** 151C 152 IF (WHERE-5) 200,40,10 153C 154C UPDATE BY ADDING A FUNCTION VALUE IN THE MIDDLE. 155C 15610 HAVDIF=.FALSE. 157 IF (NFEVAL.GT.NFJUMP+6) THEN 158 if (l .le. 1) GO TO 70 159 WHERE=0 160 L=MIN(L,LENDT+1) 161 IF (L.GE.LENDT+1) GO TO 180 162 I=LENDT 163 EPSCOR=0.5d0*(-ABS(XT(L-1)*(FT(L)-FT(L-1))) 164 1 +ABS(XT(L-1)*(FNCVAL-FT(L-1)))+ABS(X*(FT(L)-FNCVAL))) 165 EPSCOR=EPSCOR*EMEPS 166 EPSCOR=EPSCOR+ABS(FNCVAL*RNDC*(XT(L)-XT(L-1))) 167 IF (FEA.NE.0) EPSCOR=EPSCOR+ABS(0.5d0*ERRF*(XT(L)-XT(L-1))) 168 EPSMIN=EPSMIN+MAX(EPSCOR,0.0D0) 169C DO FOREVER 17020 CONTINUE 171 XT(I+1)=XT(I) 172 FT(I+1)=FT(I) 173 IF (I.EQ.L) GO TO 180 174 I=I-1 175 GO TO 20 176C END FOREVER 177 END IF 178 IF (L-LENDT-1) 80,50,50 179C 180C UPDATE BY ADDING A FUNCTION VALUE ON THE BLOCAL END. 181C 18240 IF (WHERE2.EQ.1) GO TO 70 183 FNCVAL=FATB 184 L=LENDT+1 185C ADD ONE AT THE BLOCAL END 18650 PHIT(L)=FNCVAL 187 TP=1.0d0 188 PHIT(LENDT)=FNCVAL-PHIT(LENDT) 189 I=LENDT 19060 TP=TP*((X-XT(I))/(XT(LENDT)-XT(I-1))) 191 PHIT(I-1)=PHIT(I)-TP*PHIT(I-1) 192 I=I-1 193 IF (I.GE.2) GO TO 60 194 GO TO 140 195C 196C UPDATE BY ADDING A FUNCTION VALUE ON THE ALOCAL END. 197C 19870 FNCVAL=FATA 199 L=1 200C ADD ONE IN THE MIDDLE OR AT THE ALOCAL END. 20180 I=LENDT 202 TP=PHIT(I)-FNCVAL 203 S=XT(I)-X 204C DO FOREVER 20590 CONTINUE 206 XT(I+1)=XT(I) 207 FT(I+1)=FT(I) 208 PHIT(I+1)=PHIT(I) 209 PHI(I+1)=PHI(I) 210 I=I-1 211 IF (I.LT.L) GO TO 100 212 TP=TP+(S/(X-XT(I)))*(TP-PHIT(I)) 213 GO TO 90 214C END FOREVER 215100 CONTINUE 216 PHIT(L)=TP 217 IF (L.EQ.1) THEN 218C ADD ONE AT THE ALOCAL END. 219 PHI(1)=FNCVAL 220 TP=1.0d0 221 PHI(2)=FNCVAL-PHI(2) 222 I=2 223110 TP=TP*((X-XT(I))/(XT(1)-XT(I+1))) 224 PHI(I+1)=PHI(I)-TP*PHI(I+1) 225 I=I+1 226 IF (I.LE.LENDT) GO TO 110 227 GO TO 180 228 END IF 229C UPDATE PHIT FOR ADDING ONE IN THE INTERIOR. 230 I=L-1 231130 PHIT(I)=PHIT(I+1)+(S/(X-XT(I)))*(PHIT(I+1)-PHIT(I)) 232 I=I-1 233 IF (I.GT.0) GO TO 130 234C UPDATE PHI FOR ADDING ONE IN THE INTERIOR OR AT THE BLOCAL END. 235140 TP=PHI(1)-FNCVAL 236 S=XT(1)-X 237 I=2 238 IF (L.NE.2) THEN 239150 TP=TP+(S/(X-XT(I)))*(TP-PHI(I)) 240 I=I+1 241 IF (I.LT.L) GO TO 150 242 END IF 243 PHI(L)=TP 244 IF (L.NE.LENDT+1) THEN 245C I = L AT THIS TIME. 246170 PHI(I+1)=PHI(I)+(S/(X-XT(I+1)))*(PHI(I)-PHI(I+1)) 247 I=I+1 248 IF (I.LE.LENDT) GO TO 170 249 END IF 250180 LENDT=LENDT+1 251 XT(L)=X 252 FT(L)=FNCVAL 253 IF (J1OLD.NE.18) THEN 254 IF (J1OLD.GE.L) J1OLD=J1OLD+1 255 END IF 256 IF (J2OLD.GE.L) J2OLD=J2OLD+1 257 IF (WHERE.NE.0) GO TO 230 258C 259C REFORM THE DIFFERENCE LINES. 260C 261200 NFJUMP=NFEVAL 262 PHI(1)=FT(1) 263 PHIT(1)=FT(2)-FT(1) 264 PHI(2)=-PHIT(1) 265 PHIT(2)=FT(2) 266 DO 220 J=3,LENDT 267 TP=1.0d0 268 S=1.0d0 269 PHIT(J)=FT(J) 270 DO 210 I=3,J 271 PHIT(J-I+2)=PHIT(J-I+3)-TP*PHIT(J-I+2) 272 TP=TP*((XT(J)-XT(J-I+2))/(XT(J-1)-XT(J-I+1))) 273 S=S*((XT(1)-XT(J-I+2))/(XT(J)-XT(J-I+2))) 274210 CONTINUE 275 PHIT(1)=PHIT(2)-TP*PHIT(1) 276 PHI(J)=-S*PHIT(1) 277220 CONTINUE 278C 279230 CONTINUE 280 RETURN 281C 282 END 283