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