1*DECK ZQCBJ
2      SUBROUTINE ZQCBJ (LUN, KPRINT, IPASS)
3C***BEGIN PROLOGUE  ZQCBJ
4C***SUBSIDIARY
5C***PURPOSE  Quick check for SLATEC subroutine
6C            ZBESJ
7C***LIBRARY   SLATEC
8C***CATEGORY  C10A4
9C***TYPE      COMPLEX (CQCBJ-C, ZQCBJ-Z)
10C***KEYWORDS  QUICK CHECK, ZBESJ
11C***AUTHOR  Amos, Don, (SNL)
12C           Goudy, Sue, (SNL)
13C           Walton, Lee, (SNL)
14C***DESCRIPTION
15C
16C *Usage:
17C
18C        INTEGER  LUN, KPRINT, IPASS
19C
20C        CALL ZQCBJ (LUN, KPRINT, IPASS)
21C
22C *Arguments:
23C
24C     LUN    :IN  is the unit number to which output is to be written.
25C
26C     KPRINT :IN  controls the amount of output, as specified in the
27C                 SLATEC Guidelines.
28C
29C     IPASS  :OUT indicates whether the test passed or failed.
30C                 A value of one is good, indicating no failures.
31C
32C *Description:
33C
34C                 *** A DOUBLE PRECISION ROUTINE ***
35C
36C   ZQCBJ is a quick check routine for the complex J Bessel function
37C    generated by subroutine ZBESJ.
38C
39C   ZQCBJ generates sequences of J Bessel functions from ZBESJ
40C    and checks them against the evaluation from the formula
41C
42C            J(FNU,Z) = 0.5*( H(1,FNU,Z) + H(2,FNU,Z) )
43C
44C    where -PI.lt.arg(Z).le.PI for abs(Z).ge.FNU.
45C
46C   For abs(Z).lt.FNU, the first N members of a sequence of length
47C    N+16 are checked against a corresponding N member sequence where
48C    both sequences are generated by ZBESJ beginning at order FNU.
49C
50C***REFERENCES  Abramowitz, M. and Stegun, I. A., Handbook
51C                 of Mathematical Functions, Dover Publications,
52C                 New York, 1964.
53C               Amos, D. E., A Subroutine Package for Bessel
54C                 Functions of a Complex Argument and Nonnegative
55C                 Order, SAND85-1018, May, 1985.
56C***ROUTINES CALLED  ZBESH, ZBESJ, ZABS, ZEXP, I1MACH, D1MACH
57C***REVISION HISTORY  (YYMMDD)
58C   830501  DATE WRITTEN
59C   890831  Revised to meet new SLATEC standards
60C   930122  Added ZEXP to EXTERNAL statement.  (RWC)
61C***END PROLOGUE  ZQCBJ
62C
63C*Internal Notes:
64C   Machine constants are defined by functions I1MACH and D1MACH.
65C
66C   The parameter MQC can have values 1 (the default) for a faster,
67C   less definitive test or 2 for a slower, more definitive test.
68C
69C**End
70C
71C  Set test complexity parameter.
72C
73      INTEGER  MQC
74      PARAMETER (MQC=1)
75C
76C  Declare arguments.
77C
78      INTEGER  LUN, KPRINT, IPASS
79C
80C  Declare external functions.
81C
82      INTEGER  I1MACH
83      DOUBLE PRECISION  D1MACH, ZABS
84      EXTERNAL  I1MACH, D1MACH, ZABS, ZEXP
85C
86C  Declare local variables.
87C
88      DOUBLE PRECISION  COE1R,COE1I, COE2R,COE2I, CWR,CWI, HALFR,HALFI,
89     *   VR,VI, WR,WI, YR,YI, ZR,ZI
90      DOUBLE PRECISION  AA, AB, AER, AI, ALIM, AR, ATOL, AV, CC, CT, DD,
91     *   DIG, ELIM, EPS, ER, ERTOL, FILM, FNU, FNUL, GNU, HPI, PI, R,
92     *   RL, RM, R1M4, R1M5, R2, SLAK, ST, STR, T, TOL, TS, XNU
93      INTEGER  I, ICASE, IERR, IL, IR, IRB, IT, ITL, K, KDO, KEPS, KK,
94     *   KODE, K1, K2, LFLG, M, MFLG, N, NL, NU, NUL, NZ, NZ1, NZ2
95      DIMENSION  AER(20), KDO(20), KEPS(20), T(20), VR(20), VI(20),
96     *   WR(20), WI(20), XNU(20), YR(20), YI(20)
97C
98C***FIRST EXECUTABLE STATEMENT  ZQCBJ
99      IF (KPRINT.GE.2) THEN
100        WRITE (LUN,99999)
10199999   FORMAT (' QUICK CHECK ROUTINE FOR THE J BESSEL FUNCTION FROM ',
102     *     'ZBESJ'/)
103      ENDIF
104C-----------------------------------------------------------------------
105C     Set parameters related to machine constants.
106C     TOL is the approximate unit roundoff limited to 1.0D-18.
107C     ELIM is the approximate exponential over- and underflow limit.
108C     exp(-ELIM).lt.exp(-ALIM)=exp(-ELIM)/TOL    and
109C     exp(ELIM).gt.exp(ALIM)=exp(ELIM)*TOL       are intervals near
110C     underflow and overflow limits where scaled arithmetic is done.
111C     RL is the lower boundary of the asymptotic expansion for large Z.
112C     DIG = number of base 10 digits in TOL = 10**(-DIG).
113C     FNUL is the lower boundary of the asymptotic series for large FNU.
114C-----------------------------------------------------------------------
115      R1M4 = D1MACH(4)
116      TOL = MAX(R1M4,1.0D-18)
117      ATOL = 100.0D0*TOL
118      AA = -LOG10(R1M4)
119      K1 = I1MACH(12)
120      K2 = I1MACH(13)
121      R1M5 = D1MACH(5)
122      K = MIN(ABS(K1),ABS(K2))
123      ELIM = 2.303D0*(K*R1M5-3.0D0)
124      AB = AA*2.303D0
125      ALIM = ELIM + MAX(-AB,-41.45D0)
126      DIG = MIN(AA,18.0D0)
127      FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
128      RL = 1.2D0*DIG + 3.0D0
129      SLAK = 3.0D0+4.0D0*(-LOG10(TOL)-7.0D0)/11.0D0
130      SLAK = MAX(SLAK,3.0D0)
131      ERTOL = TOL*10.0D0**SLAK
132      RM = 0.5D0*(ALIM + ELIM)
133      RM = MIN(RM,200.0D0)
134      RM = MAX(RM,RL+10.0D0)
135      R2 = MIN(RM,FNUL)
136      IF (KPRINT.GE.2) THEN
137        WRITE (LUN,99998)
13899998   FORMAT (' PARAMETERS'/
139     *     5X,'TOL ',8X,'ELIM',8X,'ALIM',8X,'RL  ',8X,'FNUL',8X,'DIG')
140        WRITE (LUN,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
14199997   FORMAT (1X,6D12.4/)
142      ENDIF
143C-----------------------------------------------------------------------
144C     Set other constants needed in the tests.
145C-----------------------------------------------------------------------
146      HALFR = 0.5D0
147      HALFI = 0.0D0
148      HPI = 2.0D0*ATAN(1.0D0)
149      PI = HPI + HPI
150C-----------------------------------------------------------------------
151C     Generate angles for construction of complex Z to be used in tests.
152C-----------------------------------------------------------------------
153C     KDO(K), K = 1,IL  determines which of the IL angles in -PI to PI
154C     are used to compute values of Z.
155C       KDO(K) = 0  means that the index K will be used for one or two
156C                   values of Z, depending on the choice of KEPS(K)
157C              = 1  means that the index K and the corresponding angle
158C                   will be skipped
159C     KEPS(K), K = 1,IL determines which of the angles get incremented
160C     up and down to put values of Z in regions where different
161C     formulae are used.
162C       KEPS(K)  = 0  means that the angle will be used without change
163C                = 1  means that the angle will be incremented up and
164C                   down by EPS
165C     The angles to be used are stored in the T(I) array, I = 1,ITL.
166C-----------------------------------------------------------------------
167      IF (MQC.NE.2) THEN
168        NL = 2
169        IL = 5
170        DO 5 I = 1,IL
171          KEPS(I) = 0
172          KDO(I) = 0
173    5   CONTINUE
174        NUL = 5
175        XNU(1) = 0.0D0
176        XNU(2) = 1.0D0
177        XNU(3) = 2.0D0
178        XNU(4) = 0.5D0*FNUL
179        XNU(5) = FNUL + 1.1D0
180      ELSE
181        NL = 4
182        IL = 13
183        DO 6 I = 1,IL
184          KDO(I) = 0
185          KEPS(I) = 0
186    6   CONTINUE
187        KDO(2) = 1
188        KDO(6) = 1
189        KDO(8) = 1
190        KDO(12) = 1
191        KEPS(3) = 1
192        KEPS(4) = 1
193        KEPS(5) = 1
194        KEPS(9) = 1
195        KEPS(10) = 1
196        KEPS(11) = 1
197        NUL = 6
198        XNU(1) = 0.0D0
199        XNU(2) = 0.6D0
200        XNU(3) = 1.3D0
201        XNU(4) = 2.0D0
202        XNU(5) = 0.5D0*FNUL
203        XNU(6) = FNUL + 1.1D0
204      ENDIF
205      I = 2
206      EPS = 0.01D0
207      FILM = IL - 1
208      T(1) = -PI + EPS
209      DO 30 K = 2,IL
210        IF (KDO(K).EQ.0) THEN
211          T(I) = PI*(-IL+2*K-1)/FILM
212          IF (KEPS(K).EQ.0) THEN
213            TS = T(I)
214            T(I) = TS - EPS
215            I = I + 1
216            T(I) = TS + EPS
217          ELSE
218            I = I + 1
219          ENDIF
220        ENDIF
221   30 CONTINUE
222      ITL = I - 1
223C-----------------------------------------------------------------------
224C     Test values of Z in -PI.lt.arg(Z).le.PI.
225C-----------------------------------------------------------------------
226      IF (KPRINT.GE.2) THEN
227        WRITE (LUN,99996)
22899996   FORMAT (' CHECKS IN THE (Z,FNU) SPACE'/)
229      ENDIF
230      LFLG = 0
231      DO 260 KODE = 1,2
232        DO 250 N = 1,NL
233          DO 240 NU = 1,NUL
234            FNU = XNU(NU)
235            DO 230 ICASE = 1,3
236              IRB = MIN(2,ICASE)
237              DO 220 IR = IRB,4
238C-------------- switch (icase)
239                GO TO (50, 60, 70), ICASE
240   50           CONTINUE
241                  R = (EPS*(4-IR)+2.0D0*(IR-1))/3.0D0
242                  GO TO 80
243   60           CONTINUE
244                  R = (2.0D0*(4-IR)+R2*(IR-1))/3.0D0
245                  GO TO 80
246   70           CONTINUE
247                  IF (R2.GE.RM) GO TO 230
248                  R = (R2*(4-IR)+RM*(IR-1))/3.0D0
249   80           CONTINUE
250C-------------- end switch
251                GNU = FNU + (N-1)
252                DO 210 IT = 1,ITL
253                  CT = COS(T(IT))
254                  ST = SIN(T(IT))
255                  IF (ABS(CT).LT.ATOL) CT = 0.0D0
256                  IF (ABS(ST).LT.ATOL) ST = 0.0D0
257                  ZR = R*CT
258                  ZI = R*ST
259                  IF (R.GE.GNU) THEN
260C------------------ Cases for abs(Z).ge.FNU+N-1
261                    CALL ZBESJ(ZR, ZI, FNU, KODE, N, VR, VI, NZ, IERR)
262C------------------ Underflow - skip test for this case.
263                    IF (NZ.NE.0) GO TO 210
264                    CALL ZBESH(ZR, ZI, FNU, KODE, 1, N, WR, WI, NZ1,
265     *                 IERR)
266                    CALL ZBESH(ZR, ZI, FNU, KODE, 2, N, YR, YI, NZ2,
267     *                 IERR)
268                    IF (KODE.EQ.2) THEN
269C-------------------- Adjust scaling of H functions.
270                      CC = -ZI - ABS(ZI)
271                      IF (CC.GT.(-ALIM)) THEN
272                        CWR = CC
273                        CWI = ZR
274                        CALL ZEXP(CWR, CWI, COE1R, COE1I)
275                      ELSE
276                        COE1R = 0.0D0
277                        COE1I = 0.0D0
278                      ENDIF
279                      DD = ZI - ABS(ZI)
280                      IF (DD.GT.(-ALIM)) THEN
281                        CWR = DD
282                        CWI = -ZR
283                        CALL ZEXP(CWR, CWI, COE2R, COE2I)
284                      ELSE
285                        COE2R = 0.0D0
286                        COE2I = 0.0D0
287                      ENDIF
288                      DO 130 KK = 1,N
289                        STR = YR(KK)*COE2R - YI(KK)*COE2I
290                        YI(KK) = YR(KK)*COE2I + YI(KK)*COE2R
291                        YR(KK) = STR
292                        STR = WR(KK)*COE1R - WI(KK)*COE1I
293                        WI(KK) = WR(KK)*COE1I + WI(KK)*COE1R
294                        WR(KK) = STR
295  130                 CONTINUE
296                    ENDIF
297                  ELSE
298C------------------ Cases for abs(Z).lt.FNU+N-1
299                    M = N + 16
300                    CALL ZBESJ(ZR, ZI, FNU, KODE, M, VR, VI, NZ, IERR)
301C------------------ Underflow at end of sequence - skip test
302                    IF (NZ.GT.10) GO TO 210
303                    CALL ZBESJ(ZR, ZI, FNU, KODE, N, WR, WI, NZ, IERR)
304                    DO 150 KK = 1,N
305                      YR(KK) = WR(KK)
306                      YI(KK) = WI(KK)
307  150               CONTINUE
308                  ENDIF
309C-----------------------------------------------------------------------
310C     If abs(Z).ge.FNU+N-1 then the error test compares J(Z<FNU) with
311C     0.5*(H1(Z,FNU) + H2(Z,FNU)).
312C     If abs(Z).lt.FNU+N-1 then the error test ensures that calculations
313C     begun in one region of the (Z,FNU) plane and carried into another
314C     region do not diverge from calculations carried out entirely in
315C     one region.  This is an internal consistency check.
316C-----------------------------------------------------------------------
317                  MFLG = 0
318                  DO 190 I = 1,N
319                    AB = FNU + I - 1
320                    AA = MAX(2.0D0,AB)
321                    CWR = (WR(I)+YR(I))*HALFR - (WI(I)+YI(I))*HALFI
322                    CWI = (WR(I)+YR(I))*HALFI + (WI(I)+YI(I))*HALFR
323                    AV = ZABS(VR(I),VI(I))
324                    AR = CWR - VR(I)
325                    AI = CWI - VI(I)
326                    ER = ZABS(AR,AI)
327                    IF (AV.NE.0.0D0) THEN
328                      IF (ZI.EQ.0.0D0) THEN
329                        IF (DABS(ZR).LT.AA) ER = ER/AV
330                      ELSE
331                        ER = ER/AV
332                      ENDIF
333                    ENDIF
334                    AER(I) = ER
335                    IF (ER.GT.ERTOL) MFLG = 1
336  190             CONTINUE
337                  IF (MFLG.NE.0) THEN
338                    IF (LFLG.EQ.0) THEN
339                      IF (KPRINT.GE.2) THEN
340                        WRITE (LUN,99995) ERTOL
34199995                   FORMAT (' CASES WHICH VIOLATE THE RELATIVE ',
342     *                     'ERROR TEST WITH ERTOL=', D12.4/)
343                        WRITE (LUN,99994)
34499994                   FORMAT (' INPUT TO ZBESJ   Z, FNU, KODE, N')
345                      ENDIF
346                      IF (KPRINT.GE.3) THEN
347                        IF (R.GE.GNU) THEN
348                          WRITE (LUN,99993)
34999993                     FORMAT (' COMPARE WITH AVERAGE OF H1 AND H2 ',
350     *                       'FUNCTIONS FOR THE SAME INPUT')
351                          WRITE (LUN,99992)
35299992                     FORMAT (' RESULTS FROM ZBESJ    NZ, V(KK)')
353                          WRITE (LUN,99991)
35499991                     FORMAT (' RESULTS FROM ZBESH   NZ1, W(KK)')
355                          WRITE (LUN,99990)
35699990                     FORMAT (' RESULTS FROM ZBESH   NZ2, Y(KK)')
357                        ELSE
358                          WRITE (LUN,99989)
35999989                     FORMAT (' RESULTS FROM ZBESJ    NZ, W(KK)')
360                        ENDIF
361                        WRITE (LUN,99988)
36299988                   FORMAT (' TEST CASE INDICES   IR, IT, ICASE'/)
363                      ENDIF
364                    ENDIF
365                    LFLG = LFLG + 1
366                    IF (KPRINT.GE.2) THEN
367                      WRITE (LUN,99987) ZR, ZI, FNU, KODE, N
36899987                 FORMAT ('   INPUT:   Z=',2D12.4,3X,'FNU=',D12.4,
369     *                   3X,'KODE=',I3,3X,'N=',I3)
370                    ENDIF
371                    IF (KPRINT.GE.3) THEN
372                      WRITE (LUN,99986) (AER(K),K=1,N)
37399986                 FORMAT ('   ERROR:   AER(K)=',4D12.4)
374                      IF (R.GE.GNU) THEN
375                        KK = MAX(NZ1,NZ2) + 1
376                        KK = MIN(N,KK)
377                        WRITE (LUN,99985) NZ, VR(KK), VI(KK)
37899985                   FORMAT (' RESULTS:   NZ=',I3,3X,'V(KK)=',2D12.4)
379                        WRITE (LUN,99984) NZ1, WR(KK), WI(KK)
38099984                   FORMAT (' RESULTS:  NZ1=',I3,3X,'W(KK)=',2D12.4)
381                        WRITE (LUN,99983) NZ2, YR(KK), YI(KK)
38299983                   FORMAT (' RESULTS:  NZ2=',I3,3X,'Y(KK)=',2D12.4)
383                      ELSE
384                        KK = N - NZ
385                        WRITE (LUN,99982) NZ, WR(KK), WI(KK)
38699982                   FORMAT (' RESULTS:   NZ=',I3,3X,'W(KK)=',2D12.4)
387                      ENDIF
388                      WRITE (LUN,99981) IR, IT, ICASE
38999981                 FORMAT ('    CASE:  IR=',I3,3X,'IT=',I3,3X,
390     *                   'ICASE=',I3/)
391                    ENDIF
392                  ENDIF
393  210           CONTINUE
394  220         CONTINUE
395  230       CONTINUE
396  240     CONTINUE
397  250   CONTINUE
398  260 CONTINUE
399      IF (KPRINT.GE.2) THEN
400        IF (LFLG.EQ.0) THEN
401          WRITE (LUN,99980)
40299980     FORMAT (' QUICK CHECKS OK')
403        ELSE
404          WRITE (LUN,99979) LFLG
40599979     FORMAT (' ***',I5,' FAILURE(S) FOR ZBESJ IN THE (Z,FNU)',
406     *       ' PLANE')
407        ENDIF
408      ENDIF
409      IPASS=0
410      IF (LFLG.EQ.0) THEN
411        IPASS=1
412      ENDIF
413      IF (IPASS.EQ.1.AND.KPRINT.GE.2) THEN
414        WRITE (LUN,99978)
41599978   FORMAT (/' ****** ZBESJ  PASSED ALL TESTS  ******'/)
416      ENDIF
417      IF (IPASS.EQ.0.AND.KPRINT.GE.1) THEN
418        WRITE (LUN,99977)
41999977   FORMAT (/' ****** ZBESJ  FAILED SOME TESTS ******'/)
420      ENDIF
421      RETURN
422      END
423