1C 2C DRIVER FOR TESTING CMLIB ROUTINES 3C DQAG DQAGI DQAGP DQAGS DQAWC DQAWF DQAWO DQAWS DQNG 4C 5C ONE INPUT DATA CARD IS REQUIRED 6C READ(LIN,1) KPRINT,TIMES 7C 1 FORMAT(I1,E10.0) 8C 9C KPRINT = 0 NO PRINTING 10C 1 NO PRINTING FOR PASSED TESTS, SHORT MESSAGE 11C FOR FAILED TESTS 12C 2 PRINT SHORT MESSAGE FOR PASSED TESTS, FULLER 13C INFORMATION FOR FAILED TESTS 14C 3 PRINT COMPLETE QUICK-CHECK RESULTS 15C 16C ***IMPORTANT NOTE*** 17C ALL QUICK CHECKS USE ROUTINES R2MACH AND D2MACH 18C TO SET THE ERROR TOLERANCES. 19C TIMES IS A CONSTANT MULTIPLIER THAT CAN BE USED TO SCALE THE 20C VALUES OF R1MACH AND D1MACH SO THAT 21C R2MACH(I) = R1MACH(I) * TIMES FOR I=3,4,5 22C D2MACH(I) = D1MACH(I) * TIMES FOR I=3,4,5 23C THIS MAKES IT EASILY POSSIBLE TO CHANGE THE ERROR TOLERANCES 24C USED IN THE QUICK CHECKS. 25C IF TIMES .LE. 0.0 THEN TIMES IS DEFAULTED TO 1.0 26C 27C ***END NOTE*** 28C 29 COMMON/UNIT/LUN 30 COMMON/MSG/ICNT,JTEST(38) 31 COMMON/XXMULT/TIMES 32 LUN=I1MACH(2) 33 LIN=I1MACH(1) 34 ITEST=1 35C 36C READ KPRINT,TIMES PARAMETERS FROM DATA CARD.. 37C 38 READ(LIN,1) KPRINT,TIMES 391 FORMAT(I1,E10.0) 40 IF(TIMES.LE.0.) TIMES=1. 41 CALL XSETUN(LUN) 42 CALL XSETF(-1) 43 CALL XERMAX(1000) 44 IF(KPRINT.LE.2) CALL XERMAX(0) 45C TEST DOUBLE PRECISION QUADPACK ROUTINES 46 CALL DQADCK(KPRINT,IPASS) 47 ITEST=ITEST*IPASS 48C 49 IF(KPRINT.GE.1.AND.ITEST.NE.1) WRITE(LUN,2) 502 FORMAT(/' ***** WARNING -- AT LEAST ONE TEST FOR SUBLIBRARY, 51 1 QUADPKD FAILED ***** ') 52 IF(KPRINT.GE.1.AND.ITEST.EQ.1) WRITE(LUN,3) 533 FORMAT(/' ----- QUADPKD SUBLIBRARY PASSED ALL TESTS ----- ') 54 END 55 DOUBLE PRECISION FUNCTION D2MACH(I) 56 DOUBLE PRECISION D1MACH 57 COMMON/XXMULT/TIMES 58 D2MACH=D1MACH(I) 59 IF(I.EQ.1.OR. I.EQ.2) RETURN 60 D2MACH = D2MACH * DBLE(TIMES) 61 RETURN 62 END 63 REAL FUNCTION R2MACH(I) 64 COMMON/XXMULT/TIMES 65 R2MACH=R1MACH(I) 66 IF(I.EQ.1.OR. I.EQ.2) RETURN 67 R2MACH = R2MACH * TIMES 68 RETURN 69 END 70 SUBROUTINE DQADCK(KPRINT,IPASS) 71C***BEGIN PROLOGUE DQADCK 72C** CALLING PROGRAM FOR QUICK CHECK ROUTINES CDQNG, CDQAG, 73C CDQAGS, CDQAGP, CDQAGI, CDQAWO, CDQAWF AND CDQAWS 74C***AUTHOR ELSIE M. CLARK 75C NATIONAL BUREAU OF STANDARDS 76C***DESCRIPTION 77C.......................................................... 78C 79C THIS PROGRAM CALLS THE QUICK CHECK ROUTINES FOR QUADPACK 80C INTEGRATORS QDNG, QDAG, QDAGS, QDAGP, QDAGI, QDAWO, QDAWF 81C AND QDAWS 82C KPRINT SPECIFIES THE AMOUNT OF PRINTING TO BE DONE BY 83C EACH QUICK CHECK ROUTINE 84C IPASS IS OUTPUT FROM EACH QUICK CHECK ROUTINE AND IS 1 85C WHEN ALL THE TESTS IN THE QUICK CHECK ROUTINE PASSED, 0 86C OTHERWISE. 87C FLAG IS 1 WHEN ALL THE TESTS PASS FOR EACH QUICK 88C ROUTINE, 0 OTHERWISE. 89C 90C.......................................................... 91C***END PROLOGUE DQADCK 92C 93 COMMON/MSG/ICNT,ITEST(38) 94 COMMON/UNIT/LUN 95 INTEGER IPASS,LUN,KPRINT,FLAG 96C LUN DETERMINES OUTPUT UNIT 97 FLAG = 1 98 IF(KPRINT.GT.2) WRITE (LUN,5) 99 5 FORMAT(' ********',/,' NOTE- (QUADPACK DOUBLE PRECISION TESTS)', 100 X /,' THIS ROUTINE TESTS VARIOUS QUADPACK OPTIONS AND PRINTS', 101 X /,' SOME ERROR DIAGNOSTICS IN SUBROUTINE XERROR. THERE ',/, 102 X ' SHOULD BE 38 ABNORMAL RETURN MESSAGES. THE LAST LINE OF ', 103 X /,' OUTPUT SHOULD BE-ALL QUADPACK QUICK CHECK ROUTINES PASSED-', 104 X /,' (IMMEDIATELY FOLLOWING THE LAST ABNORMAL RETURN MESSAGE).', 105 X /,' ANY OTHER FINAL LINE INDICATES AN ERROR.',/,'********') 106 CALL CDQNG(LUN,KPRINT,IPASS) 107 IF (IPASS.NE.0) GO TO 10 108 IF(KPRINT.GE.1)WRITE(LUN,15) 109 15 FORMAT(' SOME TEST(S) IN CDQNG FAILED ') 110 FLAG = 0 111 10 CONTINUE 112 CALL CDQAG(LUN,KPRINT,IPASS) 113 IF (IPASS.NE.0) GO TO 20 114 IF(KPRINT.GE.1)WRITE(LUN,25) 115 25 FORMAT(' SOME TEST(S) IN CDQAG FAILED ') 116 FLAG = 0 117 20 CONTINUE 118 CALL CDQAGS(LUN,KPRINT,IPASS) 119 IF (IPASS.NE.0) GO TO 30 120 IF(KPRINT.GE.1)WRITE(LUN,35) 121 35 FORMAT(' SOME TEST(S) IN CDQAGS FAILED ') 122 FLAG = 0 123 30 CONTINUE 124 CALL CDQAGP(LUN,KPRINT,IPASS) 125 IF (IPASS.NE.0) GO TO 40 126 IF(KPRINT.GE.1)WRITE(LUN,45) 127 45 FORMAT(' SOME TEST(S) IN CDQAGP FAILED ') 128 FLAG = 0 129 40 CONTINUE 130 CALL CDQAGI(LUN,KPRINT,IPASS) 131 IF (IPASS.NE.0) GO TO 50 132 IF(KPRINT.GE.1)WRITE(LUN,55) 133 55 FORMAT(' SOME TEST(S) IN CDQAGI FAILED ') 134 FLAG = 0 135 50 CONTINUE 136 CALL CDQAWO(LUN,KPRINT,IPASS) 137 IF (IPASS.NE.0) GO TO 60 138 IF(KPRINT.GE.1)WRITE(LUN,65) 139 65 FORMAT(' SOME TEST(S) IN CDQAWO FAILED ') 140 FLAG = 0 141 60 CONTINUE 142 CALL CDQAWF(LUN,KPRINT,IPASS) 143 IF (IPASS.NE.0) GO TO 70 144 IF(KPRINT.GE.1)WRITE(LUN,75) 145 75 FORMAT(' SOME TEST(S) IN CDQAWF FAILED ') 146 FLAG = 0 147 70 CONTINUE 148 CALL CDQAWS(LUN,KPRINT,IPASS) 149 IF (IPASS.NE.0) GO TO 80 150 IF(KPRINT.GE.1)WRITE(LUN,85) 151 85 FORMAT(' SOME TEST(S) IN CDQAWS FAILED ') 152 FLAG = 0 153 80 CONTINUE 154 CALL CDQAWC(LUN,KPRINT,IPASS) 155 IF (IPASS.NE.0) GO TO 90 156 IF(KPRINT.GE.1)WRITE(LUN,95) 157 95 FORMAT(' SOME TEST(S) IN CDQAWC FAILED ') 158 FLAG = 0 159 90 CONTINUE 160 IPASS=FLAG 161 IF(KPRINT.GT.1.AND.IPASS.EQ.1) WRITE(LUN,104) 162 IF(KPRINT.NE.0.AND.IPASS.EQ.0) WRITE(LUN,105) 163 104 FORMAT(' ***** ALL DP QUADPACK ROUTINES PASSED ***** ') 164 105 FORMAT(' ***** SOME DP QUADPACK TESTS FAILED ***** ') 165 RETURN 166 END 167 SUBROUTINE CDQAG(KO,KPRINT,IPASS) 168C 169C SUBROUTINES OR FUNCTIONS NEEDED- D2MACH,DQAG,DF1G,DF2G,DF3G,DPRIN 170C 171C FOR FURTHER DOCUMENTATION SEE ROUTINE CQPDOC 172C 173 DOUBLE PRECISION A,ABSERR,B,D2MACH,EPMACH,EPSABS,EPSREL,ERROR, 174 *EXACT1,EXACT2,EXACT3,DF1G,DF2G,DF3G,PI,RESULT,UFLOW,WORK 175 INTEGER IER,IP,IPASS,IWORK,KEY,KPRINT,LAST,LENW,LIMIT, 176 * NEVAL 177 DIMENSION IERV(2),IWORK(100),WORK(400) 178 EXTERNAL DF1G,DF2G,DF3G 179 DATA PI/0.31415926535897932D+01/ 180 DATA EXACT1/0.1154700538379252D+01/ 181 DATA EXACT2/0.11780972450996172D+00/ 182 DATA EXACT3/0.1855802D+02/ 183C 184C TEST ON IER = 0 185C 186 IPASS = 1 187 LIMIT = 100 188 LENW = LIMIT*4 189 EPSABS = 0.0D+00 190 EPMACH = D2MACH(4) 191 KEY = 6 192 EPSREL = DMAX1(DSQRT(EPMACH),0.1D-07) 193 A = 0.0D+00 194 B = 0.1D+01 195 CALL DQAG(DF1G,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,NEVAL,IER, 196 *LIMIT,LENW,LAST,IWORK,WORK) 197 IERV(1) = IER 198 IP = 0 199 ERROR = DABS(EXACT1-RESULT) 200 IF(IER.EQ.0.AND.ERROR.LE.ABSERR.AND.ABSERR.LE.EPSREL*DABS 201 *(EXACT1))IP = 1 202 IF(IP.EQ.0) IPASS = 0 203 CALL DPRIN(KO,0,KPRINT,IP,EXACT1,RESULT,ABSERR,NEVAL,IERV,1) 204C 205C TEST ON IER = 1 206C 207 LIMIT = 1 208 LENW = LIMIT*4 209 B = PI*0.2D+01 210 CALL DQAG(DF2G,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,NEVAL,IER, 211 * LIMIT,LENW,LAST,IWORK,WORK) 212 IERV(1) = IER 213 IP = 0 214 IF(IER.EQ.1) IP = 1 215 IF(IP.EQ.0) IPASS = 0 216 CALL DPRIN(KO,1,KPRINT,IP,EXACT2,RESULT,ABSERR,NEVAL,IERV,1) 217C 218C TEST ON IER = 2 OR 1 219C 220 UFLOW = D2MACH(1) 221 LIMIT = 100 222 LENW = LIMIT*4 223 CALL DQAG(DF2G,A,B,UFLOW,0.0D+00,KEY,RESULT,ABSERR,NEVAL,IER, 224 * LIMIT,LENW,LAST,IWORK,WORK) 225 IERV(1) = IER 226 IERV(2) = 1 227 IP = 0 228 IF(IER.EQ.2.OR.IER.EQ.1) IP = 1 229 IF(IP.EQ.0) IPASS = 0 230 CALL DPRIN(KO,2,KPRINT,IP,EXACT2,RESULT,ABSERR,NEVAL,IERV,2) 231C 232C TEST ON IER = 3 OR 1 233C 234 B = 0.1D+01 235 CALL DQAG(DF3G,A,B,EPSABS,EPSREL,1,RESULT,ABSERR,NEVAL,IER, 236 * LIMIT,LENW,LAST,IWORK,WORK) 237 IERV(1) = IER 238 IERV(2) = 1 239 IP = 0 240 IF(IER.EQ.3.OR.IER.EQ.1) IP = 1 241 IF(IP.EQ.0) IPASS = 0 242 CALL DPRIN(KO,3,KPRINT,IP,EXACT3,RESULT,ABSERR,NEVAL,IERV,2) 243C 244C TEST ON IER = 6 245C 246 LENW = 1 247 CALL DQAG(DF1G,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,NEVAL,IER, 248 * LIMIT,LENW,LAST,IWORK,WORK) 249 IERV(1) = IER 250 IP = 0 251 IF(IER.EQ.6.AND.RESULT.EQ.0.0D+00.AND.ABSERR.EQ.0.0D+00.AND. 252 * NEVAL.EQ.0.AND.LAST.EQ.0) IP = 1 253 IF(IP.EQ.0) IPASS = 0 254 CALL DPRIN(KO,6,KPRINT,IP,EXACT1,RESULT,ABSERR,NEVAL,IERV,1) 255 RETURN 256 END 257 SUBROUTINE CDQAGI(KO,KPRINT,IPASS) 258C 259C SUBROUTINES OR FUNCTIONS NEEDED- D2MACH,DQAGI,DT0,DT1,DT2,DT3,DT4 260C DT5,DPRIN,DF0S,DF1S,DF2S,DF3S,DF4S,DF5S 261C 262C FOR FURTHER DOCUMENTATION SEE ROUTINE CQPDOC 263C 264 DOUBLE PRECISION ABSERR,BOUND,D2MACH,EPMACH,EPSABS, 265 * EPSREL,ERROR,EXACT0,EXACT1,EXACT2,EXACT3,EXACT4, 266 * OFLOW,RESULT,DT0,DT1,DT2,DT3,DT4,DT5,UFLOW,WORK 267 INTEGER IER,IP,IPASS,IWORK,KO,KPRINT,LAST,LENW,LIMIT,NEVAL 268 DIMENSION WORK(800),IWORK(200),IERV(4) 269 EXTERNAL DT0,DT1,DT2,DT3,DT4,DT5 270 DATA EXACT0/2.0D+00/,EXACT1/0.115470066904D1/ 271 DATA EXACT2/0.909864525656D-02/ 272 DATA EXACT3/0.31415926535897932D+01/ 273 DATA EXACT4/0.19984914554328673D+04/ 274C 275C TEST ON IER = 0 276C 277 IPASS = 1 278 LIMIT = 200 279 LENW = LIMIT*4 280 EPSABS = 0.0D+00 281 EPMACH = D2MACH(4) 282 EPSREL = DMAX1(DSQRT(EPMACH),0.1D-07) 283 BOUND = 0.0D+00 284 INF = 1 285 CALL DQAGI(DT0,BOUND,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER, 286 * LIMIT,LENW,LAST,IWORK,WORK) 287 ERROR = DABS(RESULT-EXACT0) 288 IERV(1) = IER 289 IP = 0 290 IF(IER.EQ.0.AND.ERROR.LE.EPSREL*DABS(EXACT0)) 291 * IP = 1 292 IF(IP.EQ.0) IPASS = 0 293 CALL DPRIN(KO,0,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 294C 295C TEST ON IER = 1 296C 297 CALL DQAGI(DT1,BOUND,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER, 298 * 1,4,LAST,IWORK,WORK) 299 IERV(1) = IER 300 IP = 0 301 IF(IER.EQ.1) IP = 1 302 IF(IP.EQ.0) IPASS = 0 303 CALL DPRIN(KO,1,KPRINT,IP,EXACT1,RESULT,ABSERR,NEVAL,IERV,1) 304C 305C TEST ON IER = 2 OR 4 OR 1 306C 307 UFLOW = D2MACH(1) 308 A = 0.1D+00 309 CALL DQAGI(DT2,BOUND,INF,UFLOW,0.0D+00,RESULT,ABSERR,NEVAL,IER, 310 * LIMIT,LENW,LAST,IWORK,WORK) 311 IERV(1) = IER 312 IERV(2) = 4 313 IERV(3) = 1 314 IP = 0 315 IF(IER.EQ.2.OR.IER.EQ.4.OR.IER.EQ.1) IP = 1 316 IF(IP.EQ.0) IPASS = 0 317 CALL DPRIN(KO,2,KPRINT,IP,EXACT2,RESULT,ABSERR,NEVAL,IERV,3) 318C 319C TEST ON IER = 3 OR 4 OR 1 OR 2 320C 321 A = 0.0D+00 322 B = 0.5D+01 323 CALL DQAGI(DT3,BOUND,INF,UFLOW,0.0D+00,RESULT,ABSERR,NEVAL,IER, 324 * LIMIT,LENW,LAST,IWORK,WORK) 325 IERV(1) = IER 326 IERV(2) = 4 327 IERV(3) = 1 328 IERV(4) = 2 329 IP = 0 330 IF(IER.EQ.3.OR.IER.EQ.4.OR.IER.EQ.1.OR.IER.EQ.2) IP = 1 331 IF(IP.EQ.0) IPASS = 0 332 CALL DPRIN(KO,3,KPRINT,IP,EXACT3,RESULT,ABSERR,NEVAL,IERV,4) 333C 334C TEST ON IER = 4 OR 3 OR 1 OR 0 335C 336 B = 0.1D+01 337 CALL DQAGI(DT4,BOUND,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER, 338 * LIMIT,LENW,LAST,IWORK,WORK) 339 IERV(1) = IER 340 IERV(2) = 3 341 IERV(3) = 1 342 IERV(4) = 0 343 IP = 0 344 IF(IER.EQ.4.OR.IER.EQ.3.OR.IER.EQ.1.OR.IER.EQ.0) IP = 1 345 IF(IP.EQ.0) IPASS = 0 346 CALL DPRIN(KO,4,KPRINT,IP,EXACT4,RESULT,ABSERR,NEVAL,IERV,4) 347C 348C TEST ON IER = 5 349C 350 OFLOW = D2MACH(2) 351 CALL DQAGI(DT5,BOUND,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER, 352 * LIMIT,LENW,LAST,IWORK,WORK) 353 IERV(1) = IER 354 IP = 0 355 IF(IER.EQ.5) IP = 1 356 IF(IP.EQ.0) IPASS = 0 357 CALL DPRIN(KO,5,KPRINT,IP,OFLOW,RESULT,ABSERR,NEVAL,IERV,1) 358C 359C TEST ON IER = 6 360C 361 CALL DQAGI(DT1,BOUND,INF,EPSABS,0.0D+00,RESULT,ABSERR,NEVAL,IER, 362 * LIMIT,LENW,LAST,IWORK,WORK) 363 IERV(1) = IER 364 IP = 0 365 IF(IER.EQ.6.AND.RESULT.EQ.0.0D+00.AND.ABSERR.EQ.0.0D+00.AND. 366 * NEVAL.EQ.0.AND.LAST.EQ.0) IP = 1 367 IF(IP.EQ.0) IPASS = 0 368 CALL DPRIN(KO,6,KPRINT,IP,EXACT1,RESULT,ABSERR,NEVAL,IERV,1) 369 RETURN 370 END 371 SUBROUTINE CDQAGP(KO,KPRINT,IPASS) 372C 373C SUBROUTINES OR FUNCTIONS NEEDED- D2MACH,DQAGP,DF1P,DF2P,DF3P 374C DF4P,DPRIN 375C 376C FOR FURTHER DOCUMENTATION SEE ROUTINE CQPDOC 377C 378 DOUBLE PRECISION A,ABSERR,B,D2MACH,EPMACH,EPSABS,EPSREL,ERROR, 379 * EXACT1, 380 * EXACT2,EXACT3,DF1P,DF2P,DF3P,DF4P,OFLOW,POINTS,P1,P2,RESULT, 381 * UFLOW,WORK 382 INTEGER IER,IP,IPASS,IWORK,KO,KPRINT,LAST,LENIW,LENW,LIMIT, 383 * NEVAL,NPTS2 384 DIMENSION IERV(4),IWORK(205),POINTS(5),WORK(405) 385 EXTERNAL DF1P,DF2P,DF3P,DF4P 386 DATA EXACT1/0.4285277667368085D+01/ 387 DATA EXACT2/0.909864525656D-2/ 388 DATA EXACT3/0.31415926535897932D+01/ 389 DATA P1/0.1428571428571428D+00/ 390 DATA P2/0.6666666666666667D+00/ 391C 392C TEST ON IER = 0 393C 394 IPASS = 1 395 NPTS2 = 4 396 LIMIT = 100 397 LENIW = LIMIT*2+NPTS2 398 LENW = LIMIT*4+NPTS2 399 EPSABS = 0.0D+00 400 EPMACH = D2MACH(4) 401 EPSREL=DMAX1(DSQRT(EPMACH),0.1D-07) 402 A = 0.0D+00 403 B = 0.1D+01 404 POINTS(1) = P1 405 POINTS(2) = P2 406 CALL DQAGP(DF1P,A,B,NPTS2,POINTS,EPSABS,EPSREL,RESULT,ABSERR, 407 * NEVAL,IER,LENIW,LENW,LAST,IWORK,WORK) 408 ERROR = DABS(RESULT-EXACT1) 409 IERV(1) = IER 410 IP=0 411 IF(IER.EQ.0.AND.ERROR.LE.EPSREL*DABS(EXACT1)) IP = 1 412 IF(IP.EQ.0) IPASS = 0 413 CALL DPRIN(KO,0,KPRINT,IP,EXACT1,RESULT,ABSERR,NEVAL,IERV,1) 414C 415C TEST ON IER = 1 416C 417 LENIW = 10 418 LENW = LENIW*2-NPTS2 419 CALL DQAGP(DF1P,A,B,NPTS2,POINTS,EPSABS,EPSREL,RESULT,ABSERR, 420 * NEVAL,IER,LENIW,LENW,LAST,IWORK,WORK) 421 IERV(1) = IER 422 IP = 0 423 IF(IER.EQ.1) IP = 1 424 IF(IP.EQ.0) IPASS = 0 425 CALL DPRIN(KO,1,KPRINT,IP,EXACT1,RESULT,ABSERR,NEVAL,IERV,1) 426C 427C TEST ON IER = 2, 4, 1 OR 3 428C 429 NPTS2 = 3 430 POINTS(1) = 0.1D+00 431 LENIW = LIMIT*2+NPTS2 432 LENW = LIMIT*4+NPTS2 433 UFLOW = D2MACH(1) 434 A = 0.1D+00 435 CALL DQAGP(DF2P,A,B,NPTS2,POINTS,UFLOW,0.0D+00,RESULT,ABSERR, 436 * NEVAL,IER,LENIW,LENW,LAST,IWORK,WORK) 437 IERV(1) = IER 438 IERV(2) = 4 439 IERV(3) = 1 440 IERV(4) = 3 441 IP = 0 442 IF(IER.EQ.2.OR.IER.EQ.4.OR.IER.EQ.1.OR.IER.EQ.3) IP = 1 443 IF(IP.EQ.0) IPASS = 0 444 CALL DPRIN(KO,2,KPRINT,IP,EXACT2,RESULT,ABSERR,NEVAL,IERV,4) 445C 446C TEST ON IER = 3 OR 4 OR 1 OR 2 447C 448 NPTS2 = 2 449 LENIW = LIMIT*2+NPTS2 450 LENW = LIMIT*4+NPTS2 451 A = 0.0D+00 452 B = 0.5D+01 453 CALL DQAGP(DF3P,A,B,NPTS2,POINTS,UFLOW,0.0D+00,RESULT,ABSERR, 454 * NEVAL,IER,LENIW,LENW,LAST,IWORK,WORK) 455 IERV(1) = IER 456 IERV(2) = 4 457 IERV(3) = 1 458 IERV(4) = 2 459 IP = 0 460 IF(IER.EQ.3.OR.IER.EQ.4.OR.IER.EQ.1.OR.IER.EQ.2) IP = 1 461 IF(IP.EQ.0) IPASS = 0 462 CALL DPRIN(KO,3,KPRINT,IP,EXACT3,RESULT,ABSERR,NEVAL,IERV,4) 463C 464C TEST ON IER = 5 465C 466 B = 0.1D+01 467 CALL DQAGP(DF4P,A,B,NPTS2,POINTS,EPSABS,EPSREL,RESULT,ABSERR, 468 * NEVAL,IER,LENIW,LENW,LAST,IWORK,WORK) 469 IERV(1) = IER 470 IP = 0 471 IF(IER.EQ.5) IP = 1 472 IF(IP.EQ.0) IPASS = 0 473 OFLOW = D2MACH(2) 474 CALL DPRIN(KO,5,KPRINT,IP,OFLOW,RESULT,ABSERR,NEVAL,IERV,1) 475C 476C TEST ON IER = 6 477C 478 NPTS2 = 5 479 LENIW = LIMIT*2+NPTS2 480 LENW = LIMIT*4+NPTS2 481 POINTS(1) = P1 482 POINTS(2) = P2 483 POINTS(3) = 0.3D+01 484 CALL DQAGP(DF1P,A,B,NPTS2,POINTS,EPSABS,EPSREL,RESULT,ABSERR, 485 * NEVAL,IER,LENIW,LENW,LAST,IWORK,WORK) 486 IERV(1) = IER 487 IP = 0 488 IF(IER.EQ.6.AND.RESULT.EQ.0.0D+00.AND.ABSERR.EQ.0.0D+00.AND. 489 * NEVAL.EQ.0.AND.LAST.EQ.0) IP = 1 490 IF(IP.EQ.0) IPASS = 0 491 CALL DPRIN(KO,6,KPRINT,IP,EXACT1,RESULT,ABSERR,NEVAL,IERV,1) 492 RETURN 493 END 494 SUBROUTINE CDQAGS(KO,KPRINT,IPASS) 495C 496C SUBROUTINES OR FUNCTIONS NEEDED- D2MACH,DQAGS,DF0S,DF1S,DF2S,DF3S 497C DF45,DF5S,DPRIN 498C 499C FOR FURTHER DOCUMENTATION SEE ROUTINE CQPDOC 500C 501 DOUBLE PRECISION A,ABSERR,B,D2MACH,EPMACH,EPSABS, 502 * EPSREL,ERROR,EXACT0,EXACT1,EXACT2,EXACT3,EXACT4, 503 * DF0S,DF1S,DF2S,DF3S,DF4S,DF5S,OFLOW,RESULT,UFLOW,WORK 504 INTEGER IER,IP,IPASS,IWORK,KPRINT,LAST,LENW,LIMIT,NEVAL 505 DIMENSION IERV(4),IWORK(200),WORK(800) 506 EXTERNAL DF0S,DF1S,DF2S,DF3S,DF4S,DF5S 507 DATA EXACT0/0.2D+01/ 508 DATA EXACT1/0.115470066904D+01/ 509 DATA EXACT2/0.909864525656D-02/ 510 DATA EXACT3/0.31415926535897932D+01/ 511 DATA EXACT4/0.19984914554328673D+04/ 512C 513C TEST ON IER = 0 514C 515 IPASS = 1 516 LIMIT = 200 517 LENW = LIMIT*4 518 EPSABS = 0.0D+00 519 EPMACH = D2MACH(4) 520 EPSREL = DMAX1(DSQRT(EPMACH),0.1D-07) 521 A = 0.0D+00 522 B = 0.1D+01 523 CALL DQAGS(DF0S,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER, 524 * LIMIT,LENW,LAST,IWORK,WORK) 525 ERROR = DABS(RESULT-EXACT0) 526 IERV(1) = IER 527 IP = 0 528 IF(IER.EQ.0.AND.ERROR.LE.EPSREL*DABS(EXACT0)) 529 * IP = 1 530 IF(IP.EQ.0) IPASS = 0 531 CALL DPRIN(KO,0,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 532C 533C TEST ON IER = 1 534C 535 CALL DQAGS(DF1S,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER, 536 * 1,4,LAST,IWORK,WORK) 537 IERV(1) = IER 538 IP = 0 539 IF(IER.EQ.1)IP = 1 540 IF(IP.EQ.0) IPASS = 0 541 CALL DPRIN(KO,1,KPRINT,IP,EXACT1,RESULT,ABSERR,NEVAL,IERV,1) 542C 543C TEST ON IER = 2 OR 4 OR 1 544C 545 UFLOW = D2MACH(1) 546 A = 0.1D+00 547 CALL DQAGS(DF2S,A,B,UFLOW,0.0D+00,RESULT,ABSERR,NEVAL,IER, 548 * LIMIT,LENW,LAST,IWORK,WORK) 549 IERV(1) = IER 550 IERV(2) = 4 551 IERV(3) = 1 552 IP = 0 553 IF(IER.EQ.2.OR.IER.EQ.4.OR.IER.EQ.1) IP = 1 554 IF(IP.EQ.0) IPASS = 0 555 CALL DPRIN(KO,2,KPRINT,IP,EXACT2,RESULT,ABSERR,NEVAL,IERV,3) 556C 557C TEST ON IER = 3 OR 4 OR 1 OR 2 558C 559 A = 0.0D+00 560 B = 0.5D+01 561 CALL DQAGS(DF3S,A,B,UFLOW,0.0D+00,RESULT,ABSERR,NEVAL,IER, 562 * LIMIT,LENW,LAST,IWORK,WORK) 563 IERV(1) = IER 564 IERV(2) = 4 565 IERV(3) = 1 566 IERV(4) = 2 567 IP = 0 568 IF(IER.EQ.3.OR.IER.EQ.4.OR.IER.EQ.1.OR.IER.EQ.2) IP = 1 569 IF(IP.EQ.0) IPASS = 0 570 CALL DPRIN(KO,3,KPRINT,IP,EXACT3,RESULT,ABSERR,NEVAL,IERV,4) 571C 572C TEST ON IER = 4 OR 3 OR 1 OR 0 573C 574 B = 0.1D+01 575 CALL DQAGS(DF4S,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER, 576 * LIMIT,LENW,LAST,IWORK,WORK) 577 IERV(1) = IER 578 IERV(2) = 3 579 IERV(3) = 1 580 IERV(4) = 0 581 IP = 0 582 IF(IER.EQ.4.OR.IER.EQ.3.OR.IER.EQ.1.OR.IER.EQ.0) IP = 1 583 IF(IP.EQ.0) IPASS = 0 584 CALL DPRIN(KO,4,KPRINT,IP,EXACT4,RESULT,ABSERR,NEVAL,IERV,4) 585C 586C TEST ON IER = 5 587C 588 OFLOW = D2MACH(2) 589 CALL DQAGS(DF5S,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER, 590 * LIMIT,LENW,LAST,IWORK,WORK) 591 IERV(1) = IER 592 IP = 0 593 IF(IER.EQ.5) IP = 1 594 IF(IP.EQ.0) IPASS = 0 595 CALL DPRIN(KO,5,KPRINT,IP,OFLOW,RESULT,ABSERR,NEVAL,IERV,1) 596C 597C TEST ON IER = 6 598C 599 CALL DQAGS(DF1S,A,B,EPSABS,0.0D+00,RESULT,ABSERR,NEVAL,IER, 600 * LIMIT,LENW,LAST,IWORK,WORK) 601 IERV(1) = IER 602 IP = 0 603 IF(IER.EQ.6.AND.RESULT.EQ.0.0D+00.AND.ABSERR.EQ.0.0D+00.AND. 604 * NEVAL.EQ.0.AND.LAST.EQ.0) IP = 1 605 IF(IP.EQ.0) IPASS = 0 606 CALL DPRIN(KO,6,KPRINT,IP,EXACT1,RESULT,ABSERR,NEVAL,IERV,1) 607 RETURN 608 END 609 SUBROUTINE CDQAWC(KO,KPRINT,IPASS) 610C 611C SUBROUTINES OR FUNCTIONS NEEDED- D2MACH,DQAWC,DF0O,DPRIN 612C FOR FURTHER DOCUMENTATION SEE ROUTINE CQPDOC 613C 614 DOUBLE PRECISION A,ABSERR,B,D2MACH,EPMACH,EPSABS, 615 * EPSREL,ERROR,EXACT0,EXACT1,DF0C,DF1C,C, 616 * RESULT,UFLOW,WORK 617 INTEGER IER,IP,IPASS,IWORK,KPRINT,LAST,LENW,LIMIT,NEVAL 618 DIMENSION WORK(800),IWORK(200),IERV(2) 619 EXTERNAL DF0C,DF1C 620 DATA EXACT0/-0.6284617285065624D+03/ 621 DATA EXACT1/0.1855802D+01/ 622C 623C TEST ON IER = 0 624C 625 IPASS = 1 626 C = 0.5D+00 627 A = -1.0D+00 628 B = 1.0D+00 629 LIMIT = 200 630 LENW = LIMIT*4 631 EPSABS = 0.0D+00 632 EPMACH = D2MACH(4) 633 EPSREL = DMAX1(DSQRT(EPMACH),0.1D-07) 634 CALL DQAWC(DF0C,A,B,C,EPSABS,EPSREL,RESULT,ABSERR, 635 * NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK) 636 IERV(1) = IER 637 IP = 0 638 ERROR = DABS(EXACT0-RESULT) 639 IF(IER.EQ.0.AND.ERROR.LE.ABSERR.AND.ABSERR.LE.EPSREL*DABS(EXACT0)) 640 * IP = 1 641 IF(IP.EQ.0) IPASS = 0 642 CALL DPRIN(KO,0,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 643C 644C TEST ON IER = 1 645C 646 CALL DQAWC(DF0C,A,B,C,EPSABS,EPSREL,RESULT,ABSERR, 647 * NEVAL,IER,1,4,LAST,IWORK,WORK) 648 IERV(1) = IER 649 IP = 0 650 IF(IER.EQ.1) IP = 1 651 IF(IP.EQ.0) IPASS = 0 652 CALL DPRIN(KO,1,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 653C 654C TEST ON IER = 2 OR 1 655C 656 UFLOW = D2MACH(1) 657 CALL DQAWC(DF0C,A,B,C,UFLOW,0.0D+00,RESULT,ABSERR, 658 * NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK) 659 IERV(1) = IER 660 IERV(2) = 1 661 IP = 0 662 IF(IER.EQ.2.OR.IER.EQ.1) IP = 1 663 IF(IP.EQ.0) IPASS = 0 664 CALL DPRIN(KO,2,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,2) 665C 666C TEST ON IER = 3 OR 1 667C 668 CALL DQAWC(DF1C,0.0D+00,B,C,UFLOW,0.0D+00,RESULT,ABSERR, 669 * NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK) 670 IERV(1) = IER 671 IERV(2) = 1 672 IP = 0 673 IF(IER.EQ.3.OR.IER.EQ.1) IP = 1 674 IF(IP.EQ.0) IPASS = 0 675 CALL DPRIN(KO,3,KPRINT,IP,EXACT1,RESULT,ABSERR,NEVAL,IERV,2) 676C 677C TEST ON IER = 6 678C 679 EPSABS = 0.0D+00 680 EPSREL = 0.0D+00 681 CALL DQAWC(DF0C,A,B,C,EPSABS,EPSREL,RESULT,ABSERR, 682 * NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK) 683 IERV(1) = IER 684 IP = 0 685 IF(IER.EQ.6) IP = 1 686 IF(IP.EQ.0) IPASS = 0 687 CALL DPRIN(KO,6,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 688 RETURN 689 END 690 SUBROUTINE CDQAWF(KO,KPRINT,IPASS) 691C 692C SUBROUTINES OR FUNCTIONS NEEDED- D2MACH,DQAWF,DF0F,DF1F,DPRIN 693C 694C FOR FURTHER DOCUMENTATION SEE ROUTINE CQPDOC 695C 696 DOUBLE PRECISION A,ABSERR,D2MACH,EPSABS,EPMACH, 697 * ERROR,EXACT0,DF0F,DF1F,OMEGA,PI,RESULT,UFLOW,WORK 698 INTEGER IER,IP,IPASS,KPRINT,LENW,LIMIT,LIMLST,LST,NEVAL 699 DIMENSION IERV(4),IWORK(450),WORK(1425) 700 EXTERNAL DF0F,DF1F 701 DATA EXACT0/0.1422552162575912D+01/ 702 DATA PI/0.31415926535897932D+01/ 703C 704C TEST ON IER = 0 705C 706 IPASS = 1 707 MAXP1 = 21 708 LIMLST = 50 709 LIMIT = 200 710 LENIW = LIMIT*2+LIMLST 711 LENW = LENIW*2+MAXP1*25 712 EPMACH = D2MACH(4) 713 EPSABS = DMAX1(DSQRT(EPMACH),0.1D-02) 714 A = 0.0D+00 715 OMEGA = 0.8D+01 716 INTEGR = 2 717 CALL DQAWF(DF0F,A,OMEGA,INTEGR,EPSABS,RESULT,ABSERR,NEVAL, 718 * IER,LIMLST,LST,LENIW,MAXP1,LENW,IWORK,WORK) 719 IERV(1) = IER 720 IP = 0 721 ERROR = DABS(EXACT0-RESULT) 722 IF(IER.EQ.0.AND.ERROR.LE.ABSERR.AND.ABSERR.LE.EPSABS) 723 * IP = 1 724 IF(IP.EQ.0) IPASS = 0 725 CALL DPRIN(KO,0,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 726C 727C TEST ON IER = 1 728C 729 LIMLST = 3 730 LENIW = 403 731 LENW = LENIW*2+MAXP1*25 732 CALL DQAWF(DF0F,A,OMEGA,INTEGR,EPSABS,RESULT,ABSERR,NEVAL, 733 * IER,LIMLST,LST,LENIW,MAXP1,LENW,IWORK,WORK) 734 IERV(1) = IER 735 IP = 0 736 IF(IER.EQ.1) IP = 1 737 IF(IP.EQ.0) IPASS = 0 738 CALL DPRIN(KO,1,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 739C 740C TEST ON IER = 3 OR 4 OR 1 OR 2 741C 742 LIMLST = 50 743 LENIW = LIMIT*2+LIMLST 744 LENW = LENIW*2+MAXP1*25 745 UFLOW = D2MACH(1) 746 CALL DQAWF(DF1F,A,0.0D+00,1,UFLOW,RESULT,ABSERR,NEVAL, 747 * IER,LIMLST,LST,LENIW,MAXP1,LENW,IWORK,WORK) 748 IERV(1) = IER 749 IERV(2) = 4 750 IERV(3) = 1 751 IERV(4) = 2 752 IP = 0 753 IF(IER.EQ.3.OR.IER.EQ.4.OR.IER.EQ.1.OR.IER.EQ.2) IP = 1 754 IF(IP.EQ.0) IPASS = 0 755 CALL DPRIN(KO,3,KPRINT,IP,PI,RESULT,ABSERR,NEVAL,IERV,4) 756C 757C TEST ON IER = 6 758C 759 LIMLST = 50 760 LENIW = 20 761 CALL DQAWF(DF0F,A,OMEGA,INTEGR,EPSABS,RESULT,ABSERR,NEVAL, 762 * IER,LIMLST,LST,LENIW,MAXP1,LENW,IWORK,WORK) 763 IERV(1) = IER 764 IP = 0 765 IF(IER.EQ.6) IP = 1 766 IF(IP.EQ.0) IPASS = 0 767 CALL DPRIN(KO,6,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 768C 769C TEST ON IER = 7 770C 771 LIMLST = 50 772 LENIW = 52 773 LENW = LENIW*2+MAXP1*25 774 CALL DQAWF(DF0F,A,OMEGA,INTEGR,EPSABS,RESULT,ABSERR,NEVAL, 775 * IER,LIMLST,LST,LENIW,MAXP1,LENW,IWORK,WORK) 776 IERV(1) = IER 777 IP = 0 778 IF(IER.EQ.7) IP = 1 779 IF(IP.EQ.0) IPASS = 0 780 CALL DPRIN(KO,7,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 781 RETURN 782 END 783 SUBROUTINE CDQAWO(KO,KPRINT,IPASS) 784C 785C SUBROUTINES OR FUNCTIONS NEEDED- D2MACH,DQAWO,DF0O,DF1O,DF2O,DPRIN 786C 787C FOR FURTHER DOCUMENTATION SEE ROUTINE CQPDOC 788C 789 DOUBLE PRECISION A,ABSERR,B,EPMACH,EPSABS, 790 * EPSREL,ERROR,EXACT0,DF0O,DF1O,DF2O, 791 * OFLOW,OMEGA,PI,RESULT,D2MACH,UFLOW,WORK 792 INTEGER IER,IERV,INTEGR,IP,IPASS,IWORK,KO,KPRINT,LAST,LENW, 793 * MAXP1,NEVAL 794 DIMENSION WORK(1325),IWORK(400),IERV(4) 795 EXTERNAL DF0O,DF1O,DF2O 796 DATA EXACT0/0.1042872789432789D+05/ 797 DATA PI/0.31415926535897932D+01/ 798C 799C TEST ON IER = 0 800C 801 IPASS = 1 802 MAXP1 = 21 803 LENIW = 400 804 LENW = LENIW*2+MAXP1*25 805 EPSABS = 0.0D+00 806 EPMACH = D2MACH(4) 807 EPSREL = DMAX1(DSQRT(EPMACH),0.1D-07) 808 A = 0.0D+00 809 B = PI 810 OMEGA = 0.1D+01 811 INTEGR = 2 812 CALL DQAWO(DF0O,A,B,OMEGA,INTEGR,EPSABS,EPSREL,RESULT,ABSERR, 813 * NEVAL,IER,LENIW,MAXP1,LENW,LAST,IWORK,WORK) 814 IERV(1) = IER 815 IP = 0 816 ERROR = DABS(EXACT0-RESULT) 817 IF(IER.EQ.0.AND.ERROR.LE.ABSERR.AND.ABSERR.LE.EPSREL*DABS(EXACT0)) 818 * IP = 1 819 IF(IP.EQ.0) IPASS = 0 820 CALL DPRIN(KO,0,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 821C 822C TEST ON IER = 1 823C 824 LENIW = 2 825 LENW = LENIW*2+MAXP1*25 826 CALL DQAWO(DF0O,A,B,OMEGA,INTEGR,EPSABS,EPSREL,RESULT,ABSERR, 827 * NEVAL,IER,LENIW,MAXP1,LENW,LAST,IWORK,WORK) 828 IERV(1) = IER 829 IP = 0 830 IF(IER.EQ.1) IP = 1 831 IF(IP.EQ.0) IPASS = 0 832 CALL DPRIN(KO,1,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 833C 834C TEST ON IER = 2 OR 4 OR 1 835C 836 UFLOW = D2MACH(1) 837 LENIW = 400 838 LENW = LENIW*2+MAXP1*25 839 CALL DQAWO(DF0O,A,B,OMEGA,INTEGR,UFLOW,0.0D+00,RESULT,ABSERR, 840 * NEVAL,IER,LENIW,MAXP1,LENW,LAST,IWORK,WORK) 841 IERV(1) = IER 842 IERV(2) = 4 843 IERV(3) = 1 844 IP = 0 845 IF(IER.EQ.2.OR.IER.EQ.4.OR.IER.EQ.1) IP = 1 846 IF(IP.EQ.0) IPASS = 0 847 CALL DPRIN(KO,2,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,3) 848C 849C TEST ON IER = 3 OR 4 OR 1 OR 2 850C 851 B = 0.5D+01 852 OMEGA = 0.0D+00 853 INTEGR = 1 854 CALL DQAWO(DF1O,A,B,OMEGA,INTEGR,UFLOW,0.0D+00,RESULT,ABSERR, 855 * NEVAL,IER,LENIW,MAXP1,LENW,LAST,IWORK,WORK) 856 IERV(1) = IER 857 IERV(2) = 4 858 IERV(3) = 1 859 IERV(4) = 2 860 IP = 0 861 IF(IER.EQ.3.OR.IER.EQ.4.OR.IER.EQ.1.OR.IER.EQ.2) IP = 1 862 IF(IP.EQ.0) IPASS = 0 863 CALL DPRIN(KO,3,KPRINT,IP,PI,RESULT,ABSERR,NEVAL,IERV,4) 864C 865C TEST ON IER = 5 866C 867 B = 0.1D+01 868 OFLOW = D2MACH(2) 869 CALL DQAWO(DF2O,A,B,OMEGA,INTEGR,EPSABS,EPSREL,RESULT,ABSERR, 870 * NEVAL,IER,LENIW,MAXP1,LENW,LAST,IWORK,WORK) 871 IERV(1) = IER 872 IP = 0 873 IF(IER.EQ.5) IP = 1 874 IF(IP.EQ.0) IPASS = 0 875 CALL DPRIN(KO,5,KPRINT,IP,OFLOW,RESULT,ABSERR,NEVAL,IERV,1) 876C 877C TEST ON IER = 6 878C 879 INTEGR = 3 880 CALL DQAWO(DF0O,A,B,OMEGA,INTEGR,EPSABS,EPSREL,RESULT,ABSERR, 881 * NEVAL,IER,LENIW,MAXP1,LENW,LAST,IWORK,WORK) 882 IERV(1) = IER 883 IP = 0 884 IF(IER.EQ.6.AND.RESULT.EQ.0.0D+00.AND.ABSERR.EQ.0.0D+00.AND. 885 * NEVAL.EQ.0.AND.LAST.EQ.0) IP = 1 886 IF(IP.EQ.0) IPASS = 0 887 CALL DPRIN(KO,6,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 888 RETURN 889 END 890 SUBROUTINE CDQAWS(KO,KPRINT,IPASS) 891C 892C SUBROUTINES OR FUNCTIONS NEEDED- D2MACH,DQAWS,DF0WS,DF1WS,DF2WS,DPRIN 893C FOR FUTHER DOCUMENTATION SEE ROUTINE CQPDOC 894C 895 DOUBLE PRECISION A,ABSERR,B,D2MACH,EPMACH,EPSABS, 896 * EPSREL,ERROR,EXACT0,EXACT1,DF0WS,DF1WS,ALFA,BETA, 897 * RESULT,UFLOW,WORK 898 INTEGER IER,IP,IPASS,IWORK,KPRINT,LAST,LENW,LIMIT,NEVAL,INTEGR 899 DIMENSION WORK(800),IWORK(200),IERV(2) 900 EXTERNAL DF0WS,DF1WS 901 DATA EXACT0/0.5350190569223644D+00/ 902 DATA EXACT1/0.1998491554328673D+04/ 903C 904C TEST ON IER = 0 905C 906 IPASS = 1 907 ALFA = -0.5D+00 908 BETA = -0.5D+00 909 INTEGR = 1 910 A = 0.0D+00 911 B = 0.1D+01 912 LIMIT = 200 913 LENW = LIMIT*4 914 EPSABS = 0.0D+00 915 EPMACH = D2MACH(4) 916 EPSREL = DMAX1(DSQRT(EPMACH),0.1D-07) 917 CALL DQAWS(DF0WS,A,B,ALFA,BETA,INTEGR,EPSABS,EPSREL,RESULT,ABSERR, 918 * NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK) 919 IERV(1) = IER 920 IP = 0 921 ERROR = DABS(EXACT0-RESULT) 922 IF(IER.EQ.0.AND.ERROR.LE.EPSREL*DABS(EXACT0)) 923 * IP = 1 924 IF(IP.EQ.0) IPASS = 0 925 CALL DPRIN(KO,0,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 926C 927C TEST ON IER = 1 928C 929 CALL DQAWS(DF0WS,A,B,ALFA,BETA,INTEGR,EPSABS,EPSREL,RESULT,ABSERR, 930 * NEVAL,IER,2,8,LAST,IWORK,WORK) 931 IERV(1) = IER 932 IP = 0 933 IF(IER.EQ.1) IP = 1 934 IF(IP.EQ.0) IPASS = 0 935 CALL DPRIN(KO,1,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 936C 937C TEST ON IER = 2 OR 1 938C 939 UFLOW = D2MACH(1) 940 CALL DQAWS(DF0WS,A,B,ALFA,BETA,INTEGR,UFLOW,0.0D+00,RESULT,ABSERR, 941 * NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK) 942 IERV(1) = IER 943 IERV(2) = 1 944 IP = 0 945 IF(IER.EQ.2.OR.IER.EQ.1) IP = 1 946 IF(IP.EQ.0) IPASS = 0 947 CALL DPRIN(KO,2,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,2) 948C 949C TEST ON IER = 3 OR 1 950C 951 CALL DQAWS(DF1WS,A,B,ALFA,BETA,INTEGR,EPSABS,EPSREL,RESULT,ABSERR, 952 * NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK) 953 IERV(1) = IER 954 IERV(2) = 1 955 IP = 0 956 IF(IER.EQ.3.OR.IER.EQ.1) IP = 1 957 IF(IP.EQ.0) IPASS = 0 958 CALL DPRIN(KO,3,KPRINT,IP,EXACT1,RESULT,ABSERR,NEVAL,IERV,2) 959C 960C TEST ON IER = 6 961C 962 INTEGR = 0 963 CALL DQAWS(DF1WS,A,B,ALFA,BETA,INTEGR,EPSABS,EPSREL,RESULT,ABSERR, 964 * NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK) 965 IERV(1) = IER 966 IP = 0 967 IF(IER.EQ.6) IP = 1 968 IF(IP.EQ.0) IPASS = 0 969 CALL DPRIN(KO,6,KPRINT,IP,EXACT0,RESULT,ABSERR,NEVAL,IERV,1) 970 RETURN 971 END 972 SUBROUTINE CDQNG(KO,KPRINT,IPASS) 973C 974C SUBROUTINES OR FUNCTIONS NEEDED* D2MACH,DQNG,DF1N,DF2N,DPRIN 975C 976C FOR FURTHER DOCUMENTATION SEE ROUTINE CQPDOC 977C 978 DOUBLE PRECISION A,ABSERR,B,D2MACH,EPMACH,EPSABS,EPSREL,EXACT1, 979 * ERROR,EXACT2,DF1N,DF2N,RESULT,UFLOW 980 INTEGER IER,IERV,IP,IPASS,KPRINT,NEVAL 981 DIMENSION IERV(1) 982 EXTERNAL DF1N,DF2N 983 DATA EXACT1/0.7281029132255818D+00/ 984 DATA EXACT2/0.1D+02/ 985C 986C TEST ON IER = 0 987C 988 IPASS = 1 989 EPSABS = 0.0D+00 990 EPMACH = D2MACH(4) 991 UFLOW = D2MACH(1) 992 EPSREL = DMAX1(DSQRT(EPMACH),0.1D-07) 993 A=0.0D+00 994 B=0.1D+01 995 CALL DQNG(DF1N,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER) 996 CALL DQNG(DF1N,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER) 997 IERV(1)=IER 998 IP = 0 999 ERROR = DABS(EXACT1-RESULT) 1000 IF(IER.EQ.0.AND.ERROR.LE.ABSERR.AND.ABSERR.LE.EPSREL*DABS(EXACT1)) 1001 * IP = 1 1002 IF(IP.EQ.0) IPASS = 0 1003 IF(KPRINT.NE.0) CALL DPRIN(KO,0,KPRINT,IP,EXACT1,RESULT,ABSERR, 1004 * NEVAL,IERV,1) 1005C 1006C TEST ON IER = 1 1007C 1008 CALL DQNG(DF2N,A,B,UFLOW,0.0D+00,RESULT,ABSERR,NEVAL,IER) 1009 IERV(1) = IER 1010 IP=0 1011 IF(IER.EQ.1) IP = 1 1012 IF(IP.EQ.0) IPASS = 0 1013 IF(KPRINT.NE.0) CALL DPRIN(KO,1,KPRINT,IP,EXACT2,RESULT,ABSERR, 1014 * NEVAL,IERV,1) 1015C 1016C TEST ON IER = 6 1017C 1018 EPSABS = 0.0D+00 1019 EPSREL = 0.0D+00 1020 CALL DQNG(DF1N,A,B,EPSABS,0.0D+00,RESULT,ABSERR,NEVAL,IER) 1021 IERV(1) = IER 1022 IP = 0 1023 IF(IER.EQ.6.AND.RESULT.EQ.0.0D+00.AND.ABSERR.EQ.0.0D+00.AND. 1024 * NEVAL.EQ.0) IP = 1 1025 IF(IP.EQ.0) IPASS = 0 1026 IF(KPRINT.NE.0) CALL DPRIN(KO,6,KPRINT,IP,EXACT1,RESULT,ABSERR, 1027 * NEVAL,IERV,1) 1028 RETURN 1029 END 1030 DOUBLE PRECISION FUNCTION DF0C(X) 1031 DOUBLE PRECISION X 1032 DF0C = 1.D0/(X*X+1.D-4) 1033 RETURN 1034 END 1035 DOUBLE PRECISION FUNCTION DF0F(X) 1036 DOUBLE PRECISION X 1037 DF0F = 0.0D+00 1038 IF(X.NE.0.0D+00) DF0F = DSIN(0.5D+02*X)/(X*DSQRT(X)) 1039 RETURN 1040 END 1041 DOUBLE PRECISION FUNCTION DF0O(X) 1042 DOUBLE PRECISION X 1043 DF0O = (0.2D+01*DSIN(X))**14 1044 RETURN 1045 END 1046 DOUBLE PRECISION FUNCTION DF0S(X) 1047 DOUBLE PRECISION X 1048 DF0S = 0.0D+00 1049 IF(X.NE.0.0D+00) DF0S = 0.1D+01/DSQRT(X) 1050 RETURN 1051 END 1052 DOUBLE PRECISION FUNCTION DF0WS(X) 1053 DOUBLE PRECISION X 1054 DF0WS = DSIN(0.1D+02*X) 1055 RETURN 1056 END 1057 DOUBLE PRECISION FUNCTION DF1C(X) 1058 DOUBLE PRECISION X 1059 DF1C = 0.0D+00 1060 IF(X.NE.0.33D+00) DF1C = (X-0.5D+00)*DABS(X-0.33D+00)**(-0.9D+00) 1061 RETURN 1062 END 1063 DOUBLE PRECISION FUNCTION DF1F(X) 1064 DOUBLE PRECISION X,X1,Y 1065 X1 = X+0.1D+01 1066 DF1F = 0.5D+01/X1/X1 1067 Y = 0.5D+01/X1 1068 IF(Y.GT.0.31415926535897932D+01) DF1F = 0.0D+00 1069 RETURN 1070 END 1071 DOUBLE PRECISION FUNCTION DF1G(X) 1072 DOUBLE PRECISION PI,X 1073 DATA PI/0.31415926535897932D+01/ 1074 DF1G = 0.2D+01/(0.2D+01+DSIN(0.1D+02*PI*X)) 1075 RETURN 1076 END 1077 DOUBLE PRECISION FUNCTION DF1N(X) 1078 DOUBLE PRECISION X 1079 DF1N=0.1D+01/(X**4+X**2+0.1D+01) 1080 RETURN 1081 END 1082 DOUBLE PRECISION FUNCTION DF1O(X) 1083 DOUBLE PRECISION X 1084 DF1O = 0.1D+01 1085 IF(X.GT.0.31415926535897932D+01) DF1O = 0.0D+00 1086 RETURN 1087 END 1088 DOUBLE PRECISION FUNCTION DF1P(X) 1089 DOUBLE PRECISION ALFA1,ALFA2,P1,P2,X,D1,D2 1090C P1 = 1/7, P2 = 2/3 1091 DATA P1/0.1428571428571428D+00/ 1092 DATA P2/0.6666666666666667D+00/ 1093 ALFA1 = -0.25D0 1094 ALFA2 = -0.5D0 1095 D1=DABS(X-P1) 1096 D2=DABS(X-P2) 1097 DF1P = 0.0D+00 1098 IF(D1.NE.0.0D+00.AND.D2.NE.0.0D+00) DF1P = D1**ALFA1+D2**ALFA2 1099 RETURN 1100 END 1101 DOUBLE PRECISION FUNCTION DF1S(X) 1102 DOUBLE PRECISION X 1103 DF1S = 0.2D+01/(0.2D+01+DSIN(0.314159D+02*X)) 1104 RETURN 1105 END 1106 DOUBLE PRECISION FUNCTION DF1WS(X) 1107 DOUBLE PRECISION X 1108 DF1WS = 0.00D+00 1109 IF(X-0.33D+00 .NE. 0.00D+00) DF1WS=DABS(X-0.33D+00)**(-0.999D+00) 1110 RETURN 1111 END 1112 DOUBLE PRECISION FUNCTION DF2G(X) 1113 DOUBLE PRECISION X 1114 DF2G = X*DSIN(0.3D+02*X)*DCOS(0.5D+02*X) 1115 RETURN 1116 END 1117 DOUBLE PRECISION FUNCTION DF2N(X) 1118 DOUBLE PRECISION X 1119 DF2N=X**(-0.9D+00) 1120 RETURN 1121 END 1122 DOUBLE PRECISION FUNCTION DF2O(X) 1123 DOUBLE PRECISION X 1124 DF2O = 0.0D+00 1125 IF(X.NE.0.0D+00) DF2O = 0.1D+01/(X*X*DSQRT(X)) 1126 RETURN 1127 END 1128 DOUBLE PRECISION FUNCTION DF2P(X) 1129 DOUBLE PRECISION X 1130 DF2P = DSIN(0.314159D+03*X)/(0.314159D+01*X) 1131 RETURN 1132 END 1133 DOUBLE PRECISION FUNCTION DF2S(X) 1134 DOUBLE PRECISION X 1135 DF2S = 0.1D+03 1136 IF(X.NE.0.0D+00) DF2S = DSIN(0.314159D+03*X)/(0.314159D+01*X) 1137 RETURN 1138 END 1139 DOUBLE PRECISION FUNCTION DF3G(X) 1140 DOUBLE PRECISION X 1141 DF3G=DABS(X-0.33D+00)**(-.90D+00) 1142 RETURN 1143 END 1144 DOUBLE PRECISION FUNCTION DF3P(X) 1145 DOUBLE PRECISION X 1146 DF3P = 0.1D+01 1147 IF(X.GT.0.31415926535897932D+01) DF3P = 0.0D+00 1148 RETURN 1149 END 1150 DOUBLE PRECISION FUNCTION DF3S(X) 1151 DOUBLE PRECISION X 1152 DF3S = 0.1D+01 1153 IF(X.GT.0.31415926535897932D+01) DF3S = 0.0D+00 1154 RETURN 1155 END 1156 DOUBLE PRECISION FUNCTION DF4P(X) 1157 DOUBLE PRECISION X 1158 DF4P = 0.0D+00 1159 IF(X.GT.0.0D+00) DF4P = 0.1D+01/(X*DSQRT(X)) 1160 RETURN 1161 END 1162 DOUBLE PRECISION FUNCTION DF4S(X) 1163 DOUBLE PRECISION X 1164 DF4S = 0.00D+00 1165 IF(X-0.33D+00 .NE. 0.00D+00) DF4S=DABS(X-0.33D+00)**(-0.999D+00) 1166 RETURN 1167 END 1168 DOUBLE PRECISION FUNCTION DF5S(X) 1169 DOUBLE PRECISION X 1170 DF5S = 0.0D+00 1171 IF(X.NE.0.0D+00) DF5S = 1.0D+00/(X*DSQRT(X)) 1172 RETURN 1173 END 1174 SUBROUTINE DPRIN(KO,NUM1,KPRINT,IP,EXACT,RESULT,ABSERR,NEVAL, 1175 * IERV,LIERV) 1176C 1177C***BEGIN PROLOGUE DPRIN 1178C***DPRIN ROUTINE FOR QUADPACK QUICK CHECKS, 10.27.81 1179C***AUTHORS ROBERT PIESSENS AND ELISE DE DONCKER 1180C APPL. MATH. AND PROGR. DIV.- K.U.LEUVEN 1181C***MODIFIED ELSIE M. CLARK 1182C NATIONAL BUREAU OF STANDARDS 1183C***DESCRIPTION 1184C..................................................................... 1185C 1186C THIS PROGRAM IS CALLED BY THE (DOUBLE PRECISION) QUADPACK QUICK 1187C CHECK ROUTINES, FOR PRINTING OUT THEIR MESSAGES. 1188C 1189C..................................................................... 1190C***END PROLOGUE 1191C 1192 DOUBLE PRECISION ABSERR,ERROR,EXACT,RESULT 1193 INTEGER IER,IERV,IP,KO,KPRINT,LIERV,NEVAL,NUM1 1194 DIMENSION IERV(LIERV) 1195 IER = IERV(1) 1196 ERROR = DABS(EXACT-RESULT) 1197 IF(KPRINT-2) 10,20,40 119810 IF(IP.EQ.1) GO TO 999 1199 GO TO 999 120020 IF(IP.EQ.1) GO TO 30 1201 WRITE(KO,900) NUM1 1202 IF(NUM1.EQ.0) WRITE(KO,901) 1203 IF(NUM1.GT.0) WRITE(KO,902) NUM1 1204 IF(LIERV.GT.1) WRITE(KO,903) (IERV(K),K=2,LIERV) 1205 IF(NUM1.EQ.6) WRITE(KO,904) 1206 WRITE(KO,905) 1207 WRITE(KO,906) 1208 IF(NUM1.NE.5) WRITE(KO,907) EXACT,RESULT,ERROR,ABSERR,IER,NEVAL 1209 IF(NUM1.EQ.5) WRITE(KO,910) RESULT,ABSERR,IER,NEVAL 1210 GO TO 999 121130 WRITE(KO,908) NUM1 1212 GO TO 999 121340 IF(IP.EQ.1) WRITE(KO,908) NUM1 1214 IF(NUM1.EQ.0) WRITE(KO,901) 1215 IF(NUM1.GT.0) WRITE(KO,902) NUM1 1216 IF(LIERV.GT.1) WRITE(KO,903) (IERV(K),K=2,LIERV) 1217 IF(NUM1.EQ.6) WRITE(KO,904) 1218 WRITE(KO,905) 1219 WRITE(KO,906) 1220 IF(NUM1.NE.5) WRITE(KO,907) EXACT,RESULT,ERROR,ABSERR,IER,NEVAL 1221 IF(NUM1.EQ.5) WRITE(KO,910) RESULT,ABSERR,IER,NEVAL 1222900 FORMAT(' TEST ON IER = ',I1,' FAILED.') 1223901 FORMAT(' WE MUST HAVE IER = 0, ERROR.LE.ABSERR AND, 1224 * ABSERR.LE.MAX(EPSABS,EPSREL*DABS(EXACT))') 1225902 FORMAT(' WE MUST HAVE IER = ',I1) 1226903 FORMAT(' OR IER = ',8(I1,2X)) 1227904 FORMAT(' RESULT, ABSERR, NEVAL AND EVENTUALLY LAST SHOULD, 1228 * BE ZERO') 1229905 FORMAT(' WE HAVE ') 1230906 FORMAT(' ',6X,'EXACT',11X,'RESULT',6X,'ERROR',4X,'ABSERR', 1231 * 4X,'IER NEVAL'/ 1232 * ' ',42X,'(EST.ERR.)(FLAG)(NO F-EVAL)') 1233907 FORMAT(' ',2(D15.7,1X),2(D11.2,1X),I4,4X,I6) 1234908 FORMAT(' TEST ON IER = ',I2,' PASSED') 1235910 FORMAT(5X,'INFINITY',4X,D15.7,11X,D11.2,I5,4X,I6) 1236999 RETURN 1237 END 1238 DOUBLE PRECISION FUNCTION DT0(X) 1239 DOUBLE PRECISION A,B,DF0S,X,X1,Y 1240 A = 0.0D+00 1241 B = 0.1D+01 1242 X1 = X+0.1D+01 1243 Y = (B-A)/X1+A 1244 DT0 = (B-A)*DF0S(Y)/X1/X1 1245 RETURN 1246 END 1247 DOUBLE PRECISION FUNCTION DT1(X) 1248 DOUBLE PRECISION A,B,DF1S,X,X1,Y 1249 A = 0.0D+00 1250 B = 0.1D+01 1251 X1 = X+0.1D+01 1252 Y = (B-A)/X1+A 1253 DT1 = (B-A)*DF1S(Y)/X1/X1 1254 RETURN 1255 END 1256 DOUBLE PRECISION FUNCTION DT2(X) 1257 DOUBLE PRECISION A,B,DF2S,X,X1,Y 1258 A = 0.1D+00 1259 B = 0.1D+01 1260 X1 = X+0.1D+01 1261 Y = (B-A)/X1+A 1262 DT2 = (B-A)*DF2S(Y)/X1/X1 1263 RETURN 1264 END 1265 DOUBLE PRECISION FUNCTION DT3(X) 1266 DOUBLE PRECISION A,B,DF3S,X,X1,Y 1267 A = 0.0D+00 1268 B = 0.5D+01 1269 X1 = X+0.1D+01 1270 Y = (B-A)/X1+A 1271 DT3 = (B-A)*DF3S(Y)/X1/X1 1272 RETURN 1273 END 1274 DOUBLE PRECISION FUNCTION DT4(X) 1275 DOUBLE PRECISION A,B,DF4S,X,X1,Y 1276 A = 0.0D+00 1277 B = 0.1D+01 1278 X1 = X+0.1D+01 1279 Y = (B-A)/X1+A 1280 DT4 = (B-A)*DF4S(Y)/X1/X1 1281 RETURN 1282 END 1283 DOUBLE PRECISION FUNCTION DT5(X) 1284 DOUBLE PRECISION A,B,DF5S,X,X1,Y 1285 A = 0.0D+00 1286 B = 0.1D+01 1287 X1 = X+0.1D+01 1288 Y = (B-A)/X1+A 1289 DT5 = (B-A)*DF5S(Y)/X1/X1 1290 RETURN 1291 END 1292