*DECK ZQCBJ SUBROUTINE ZQCBJ (LUN, KPRINT, IPASS) C***BEGIN PROLOGUE ZQCBJ C***SUBSIDIARY C***PURPOSE Quick check for SLATEC subroutine C ZBESJ C***LIBRARY SLATEC C***CATEGORY C10A4 C***TYPE COMPLEX (CQCBJ-C, ZQCBJ-Z) C***KEYWORDS QUICK CHECK, ZBESJ C***AUTHOR Amos, Don, (SNL) C Goudy, Sue, (SNL) C Walton, Lee, (SNL) C***DESCRIPTION C C *Usage: C C INTEGER LUN, KPRINT, IPASS C C CALL ZQCBJ (LUN, KPRINT, IPASS) C C *Arguments: C C LUN :IN is the unit number to which output is to be written. C C KPRINT :IN controls the amount of output, as specified in the C SLATEC Guidelines. C C IPASS :OUT indicates whether the test passed or failed. C A value of one is good, indicating no failures. C C *Description: C C *** A DOUBLE PRECISION ROUTINE *** C C ZQCBJ is a quick check routine for the complex J Bessel function C generated by subroutine ZBESJ. C C ZQCBJ generates sequences of J Bessel functions from ZBESJ C and checks them against the evaluation from the formula C C J(FNU,Z) = 0.5*( H(1,FNU,Z) + H(2,FNU,Z) ) C C where -PI.lt.arg(Z).le.PI for abs(Z).ge.FNU. C C For abs(Z).lt.FNU, the first N members of a sequence of length C N+16 are checked against a corresponding N member sequence where C both sequences are generated by ZBESJ beginning at order FNU. C C***REFERENCES Abramowitz, M. and Stegun, I. A., Handbook C of Mathematical Functions, Dover Publications, C New York, 1964. C Amos, D. E., A Subroutine Package for Bessel C Functions of a Complex Argument and Nonnegative C Order, SAND85-1018, May, 1985. C***ROUTINES CALLED ZBESH, ZBESJ, ZABS, ZEXP, I1MACH, D1MACH C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 890831 Revised to meet new SLATEC standards C 930122 Added ZEXP to EXTERNAL statement. (RWC) C***END PROLOGUE ZQCBJ C C*Internal Notes: C Machine constants are defined by functions I1MACH and D1MACH. C C The parameter MQC can have values 1 (the default) for a faster, C less definitive test or 2 for a slower, more definitive test. C C**End C C Set test complexity parameter. C INTEGER MQC PARAMETER (MQC=1) C C Declare arguments. C INTEGER LUN, KPRINT, IPASS C C Declare external functions. C INTEGER I1MACH DOUBLE PRECISION D1MACH, ZABS EXTERNAL I1MACH, D1MACH, ZABS, ZEXP C C Declare local variables. C DOUBLE PRECISION COE1R,COE1I, COE2R,COE2I, CWR,CWI, HALFR,HALFI, * VR,VI, WR,WI, YR,YI, ZR,ZI DOUBLE PRECISION AA, AB, AER, AI, ALIM, AR, ATOL, AV, CC, CT, DD, * DIG, ELIM, EPS, ER, ERTOL, FILM, FNU, FNUL, GNU, HPI, PI, R, * RL, RM, R1M4, R1M5, R2, SLAK, ST, STR, T, TOL, TS, XNU INTEGER I, ICASE, IERR, IL, IR, IRB, IT, ITL, K, KDO, KEPS, KK, * KODE, K1, K2, LFLG, M, MFLG, N, NL, NU, NUL, NZ, NZ1, NZ2 DIMENSION AER(20), KDO(20), KEPS(20), T(20), VR(20), VI(20), * WR(20), WI(20), XNU(20), YR(20), YI(20) C C***FIRST EXECUTABLE STATEMENT ZQCBJ IF (KPRINT.GE.2) THEN WRITE (LUN,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE J BESSEL FUNCTION FROM ', * 'ZBESJ'/) ENDIF C----------------------------------------------------------------------- C Set parameters related to machine constants. C TOL is the approximate unit roundoff limited to 1.0D-18. C ELIM is the approximate exponential over- and underflow limit. C exp(-ELIM).lt.exp(-ALIM)=exp(-ELIM)/TOL and C exp(ELIM).gt.exp(ALIM)=exp(ELIM)*TOL are intervals near C underflow and overflow limits where scaled arithmetic is done. C RL is the lower boundary of the asymptotic expansion for large Z. C DIG = number of base 10 digits in TOL = 10**(-DIG). C FNUL is the lower boundary of the asymptotic series for large FNU. C----------------------------------------------------------------------- R1M4 = D1MACH(4) TOL = MAX(R1M4,1.0D-18) ATOL = 100.0D0*TOL AA = -LOG10(R1M4) K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = D1MACH(5) K = MIN(ABS(K1),ABS(K2)) ELIM = 2.303D0*(K*R1M5-3.0D0) AB = AA*2.303D0 ALIM = ELIM + MAX(-AB,-41.45D0) DIG = MIN(AA,18.0D0) FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) RL = 1.2D0*DIG + 3.0D0 SLAK = 3.0D0+4.0D0*(-LOG10(TOL)-7.0D0)/11.0D0 SLAK = MAX(SLAK,3.0D0) ERTOL = TOL*10.0D0**SLAK RM = 0.5D0*(ALIM + ELIM) RM = MIN(RM,200.0D0) RM = MAX(RM,RL+10.0D0) R2 = MIN(RM,FNUL) IF (KPRINT.GE.2) THEN WRITE (LUN,99998) 99998 FORMAT (' PARAMETERS'/ * 5X,'TOL ',8X,'ELIM',8X,'ALIM',8X,'RL ',8X,'FNUL',8X,'DIG') WRITE (LUN,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (1X,6D12.4/) ENDIF C----------------------------------------------------------------------- C Set other constants needed in the tests. C----------------------------------------------------------------------- HALFR = 0.5D0 HALFI = 0.0D0 HPI = 2.0D0*ATAN(1.0D0) PI = HPI + HPI C----------------------------------------------------------------------- C Generate angles for construction of complex Z to be used in tests. C----------------------------------------------------------------------- C KDO(K), K = 1,IL determines which of the IL angles in -PI to PI C are used to compute values of Z. C KDO(K) = 0 means that the index K will be used for one or two C values of Z, depending on the choice of KEPS(K) C = 1 means that the index K and the corresponding angle C will be skipped C KEPS(K), K = 1,IL determines which of the angles get incremented C up and down to put values of Z in regions where different C formulae are used. C KEPS(K) = 0 means that the angle will be used without change C = 1 means that the angle will be incremented up and C down by EPS C The angles to be used are stored in the T(I) array, I = 1,ITL. C----------------------------------------------------------------------- IF (MQC.NE.2) THEN NL = 2 IL = 5 DO 5 I = 1,IL KEPS(I) = 0 KDO(I) = 0 5 CONTINUE NUL = 5 XNU(1) = 0.0D0 XNU(2) = 1.0D0 XNU(3) = 2.0D0 XNU(4) = 0.5D0*FNUL XNU(5) = FNUL + 1.1D0 ELSE NL = 4 IL = 13 DO 6 I = 1,IL KDO(I) = 0 KEPS(I) = 0 6 CONTINUE KDO(2) = 1 KDO(6) = 1 KDO(8) = 1 KDO(12) = 1 KEPS(3) = 1 KEPS(4) = 1 KEPS(5) = 1 KEPS(9) = 1 KEPS(10) = 1 KEPS(11) = 1 NUL = 6 XNU(1) = 0.0D0 XNU(2) = 0.6D0 XNU(3) = 1.3D0 XNU(4) = 2.0D0 XNU(5) = 0.5D0*FNUL XNU(6) = FNUL + 1.1D0 ENDIF I = 2 EPS = 0.01D0 FILM = IL - 1 T(1) = -PI + EPS DO 30 K = 2,IL IF (KDO(K).EQ.0) THEN T(I) = PI*(-IL+2*K-1)/FILM IF (KEPS(K).EQ.0) THEN TS = T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS ELSE I = I + 1 ENDIF ENDIF 30 CONTINUE ITL = I - 1 C----------------------------------------------------------------------- C Test values of Z in -PI.lt.arg(Z).le.PI. C----------------------------------------------------------------------- IF (KPRINT.GE.2) THEN WRITE (LUN,99996) 99996 FORMAT (' CHECKS IN THE (Z,FNU) SPACE'/) ENDIF LFLG = 0 DO 260 KODE = 1,2 DO 250 N = 1,NL DO 240 NU = 1,NUL FNU = XNU(NU) DO 230 ICASE = 1,3 IRB = MIN(2,ICASE) DO 220 IR = IRB,4 C-------------- switch (icase) GO TO (50, 60, 70), ICASE 50 CONTINUE R = (EPS*(4-IR)+2.0D0*(IR-1))/3.0D0 GO TO 80 60 CONTINUE R = (2.0D0*(4-IR)+R2*(IR-1))/3.0D0 GO TO 80 70 CONTINUE IF (R2.GE.RM) GO TO 230 R = (R2*(4-IR)+RM*(IR-1))/3.0D0 80 CONTINUE C-------------- end switch GNU = FNU + (N-1) DO 210 IT = 1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (ABS(CT).LT.ATOL) CT = 0.0D0 IF (ABS(ST).LT.ATOL) ST = 0.0D0 ZR = R*CT ZI = R*ST IF (R.GE.GNU) THEN C------------------ Cases for abs(Z).ge.FNU+N-1 CALL ZBESJ(ZR, ZI, FNU, KODE, N, VR, VI, NZ, IERR) C------------------ Underflow - skip test for this case. IF (NZ.NE.0) GO TO 210 CALL ZBESH(ZR, ZI, FNU, KODE, 1, N, WR, WI, NZ1, * IERR) CALL ZBESH(ZR, ZI, FNU, KODE, 2, N, YR, YI, NZ2, * IERR) IF (KODE.EQ.2) THEN C-------------------- Adjust scaling of H functions. CC = -ZI - ABS(ZI) IF (CC.GT.(-ALIM)) THEN CWR = CC CWI = ZR CALL ZEXP(CWR, CWI, COE1R, COE1I) ELSE COE1R = 0.0D0 COE1I = 0.0D0 ENDIF DD = ZI - ABS(ZI) IF (DD.GT.(-ALIM)) THEN CWR = DD CWI = -ZR CALL ZEXP(CWR, CWI, COE2R, COE2I) ELSE COE2R = 0.0D0 COE2I = 0.0D0 ENDIF DO 130 KK = 1,N STR = YR(KK)*COE2R - YI(KK)*COE2I YI(KK) = YR(KK)*COE2I + YI(KK)*COE2R YR(KK) = STR STR = WR(KK)*COE1R - WI(KK)*COE1I WI(KK) = WR(KK)*COE1I + WI(KK)*COE1R WR(KK) = STR 130 CONTINUE ENDIF ELSE C------------------ Cases for abs(Z).lt.FNU+N-1 M = N + 16 CALL ZBESJ(ZR, ZI, FNU, KODE, M, VR, VI, NZ, IERR) C------------------ Underflow at end of sequence - skip test IF (NZ.GT.10) GO TO 210 CALL ZBESJ(ZR, ZI, FNU, KODE, N, WR, WI, NZ, IERR) DO 150 KK = 1,N YR(KK) = WR(KK) YI(KK) = WI(KK) 150 CONTINUE ENDIF C----------------------------------------------------------------------- C If abs(Z).ge.FNU+N-1 then the error test compares J(Z