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