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