*DECK ZQCBY SUBROUTINE ZQCBY (LUN, KPRINT, IPASS) C***BEGIN PROLOGUE ZQCBY C***SUBSIDIARY C***PURPOSE Quick check for SLATEC subroutine C ZBESY C***LIBRARY SLATEC C***CATEGORY C10A4 C***TYPE COMPLEX (CQCBY-C, ZQCBY-Z) C***KEYWORDS QUICK CHECK, ZBESY 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 ZQCBY (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 ZQCBY is a quick check routine for the complex Y Bessel function C generated by subroutine ZBESY. C C ZQCBY generates sequences of Y Bessel functions from ZBESY C and checks them against the evaluation from the formula C C Y(FNU,Z*ROT) = C(FNU+1)*I(FNU,Z)-(2/PI)*CONJG(C(FNU))*K(FNU,Z) C C where ROT = EXP(PI*I/2) , C(FNU)=EXP(PI*FNU*I/2) , I**2=-1 C C and -PI.LT.ARG(Z).LE.PI/2, in the (Z,FNU) space. 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 ZBESI, ZBESK, ZBESY, 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 ZQCBY 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 CIPR,CIPI, COE1R,COE1I, COE2R,COE2I, * CSGNR,CSGNI, CSPNR,CSPNI, CWR,CWI, CWRKR,CWRKI, VR,VI, WR,WI, * YR,YI, ZR,ZI, ZNR,ZNI DOUBLE PRECISION AA, AB, AER, AI, ALIM, AR, ARG, ATOL, AV, CC, * CT, DIG, ELIM, EPS, ER, ERTOL, FFNU, FILM, FNU, FNUL, HPI, PI, * PTR, R, RHPI, RL, RM, R1M4, R1M5, R2, SLAK, ST, STI, STR, T, * TOL, TS, XNU INTEGER I, ICASE, IERR, IFNU, IL, IR, IRB, IT, ITL, I4, K, KDO, * KEPS, KK, KODE, K1, K2, LFLG, MFLG, N, NL, NU, NUL, NZ, NZ1, * NZ2 DIMENSION AER(20), CIPR(4), CIPI(4), CWRKR(20), CWRKI(20), * KDO(20), KEPS(20), T(20), VR(20), VI(20), WR(20), WI(20), * XNU(20), YR(20), YI(20) DATA CIPR(1), CIPI(1), CIPR(2), CIPI(2), CIPR(3), CIPI(3), * CIPR(4), CIPI(4) / 1.0D0,0.0D0, 0.0D0,1.0D0, -1.0D0,0.0D0, * 0.0D0,-1.0D0 / C C***FIRST EXECUTABLE STATEMENT ZQCBY IF (KPRINT.GE.2) THEN WRITE (LUN,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE Y BESSEL FUNCTION FROM ', * 'ZBESY'/) 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----------------------------------------------------------------------- HPI = 2.0D0*ATAN(1.0D0) RHPI = 1.0D0/HPI 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 KDO(5) = 1 NUL = 5 XNU(1) = 0.0D0 XNU(2) = 1.0D0 XNU(3) = 2.0D0 XNU(4) = 0.5D0*FNUL XNU(5) = FNUL + 1.2D0 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(11) = 1 KDO(12) = 1 KDO(13) = 1 KEPS(3) = 1 KEPS(4) = 1 KEPS(5) = 1 KEPS(9) = 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.2D0 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).NE.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/2.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 190 KODE = 1,2 DO 180 N = 1,NL DO 170 NU = 1,NUL C----------------------------------------------------------------------- C Construct values which will be used to set C COE1 = exp(i*(FNU+1)*PI/2) and C COE2 = (2/pi)*exp(-i*FNU*PI/2). C----------------------------------------------------------------------- FNU = XNU(NU) IFNU = INT(FNU) FFNU = FNU - IFNU ARG = HPI*FFNU CSGNR = COS(ARG) CSGNI = SIN(ARG) I4 = MOD(IFNU,4) + 1 STR = CSGNR*CIPR(I4) - CSGNI*CIPI(I4) CSGNI = CSGNR*CIPI(I4) + CSGNI*CIPR(I4) CSGNR = STR CSPNR = CSGNR*RHPI CSPNI = -CSGNI*RHPI C---------- CSGN=CSGN*CI in CQCBY STR = -CSGNI CSGNI = CSGNR CSGNR = STR DO 160 ICASE = 1,3 IRB = MIN(2,ICASE) DO 150 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 160 R = (R2*(4-IR)+RM*(IR-1))/3.0D0 80 CONTINUE C-------------- end switch DO 140 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 CALL ZBESI(ZR, ZI, FNU, KODE, N, WR, WI, NZ2, IERR) C---------------- Underflow in ZBESI - skip test for this case. IF (NZ2.NE.0) GO TO 140 CALL ZBESK(ZR, ZI, FNU, KODE, N, YR, YI, NZ1, IERR) C---------------- Underflow in ZBESK - skip test for this case. IF (NZ1.NE.0) GO TO 140 ZNR = -ZI ZNI = ZR CALL ZBESY(ZNR, ZNI, FNU, KODE, N, VR, VI, NZ, CWRKR, * CWRKI, IERR) C---------------- Underflow in ZBESY - skip test for this case. IF (NZ.NE.0) GO TO 140 COE1R = CSGNR COE1I = CSGNI COE2R = CSPNR COE2I = CSPNI IF (KODE.EQ.2) THEN C------------------ Adjust scale for I and K functions. CC = -ZR - ABS(ZR) IF (CC.GT.(-ALIM)) THEN ZNR = CC ZNI = -ZI CALL ZEXP(ZNR, ZNI, STR, STI) PTR = STR*COE2R - STI*COE2I COE2I = STR*COE2I + STI*COE2R COE2R = PTR ELSE C-------------------- Scaling problem - skip test for this case COE2R = 0.0D0 COE2I = 0.0D0 GO TO 140 ENDIF ENDIF DO 110 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 STR = -COE1I COE1I = COE1R COE1R = STR STR = COE2I COE2I = -COE2R COE2R = STR 110 CONTINUE C----------------------------------------------------------------------- C Compare Y(ZN,FNU) with COE1*I(Z,FNU)-COE2*K(Z,FNU). C----------------------------------------------------------------------- MFLG = 0 DO 120 I = 1,N AB = FNU + I - 1 AA = MAX(0.5D0,AB) CWR = WR(I) - YR(I) CWI = WI(I) - YI(I) AV = ZABS(VR(I),VI(I)) AR = CWR - VR(I) AI = CWI - VI(I) ER = ZABS(AR,AI) IF (AV.NE.0.0D0) THEN IF (ZNI.EQ.0.0D0) THEN IF (ZNR.GT.0.0D0) THEN IF (DABS(ZNR).LT.AA) ER = ER/AV ELSE IF (DABS(FFNU-0.5D0).LT.0.125D0) THEN IF (DABS(ZNR).LT.AA) ER = ER/AV ELSE ER = ER/AV ENDIF ENDIF ELSE ER = ER/AV ENDIF ENDIF AER(I) = ER IF (ER.GT.ERTOL) MFLG = 1 120 CONTINUE IF (MFLG.NE.0) THEN IF (LFLG.EQ.0) THEN IF (KPRINT.GE.2) THEN WRITE (LUN,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ', * 'ERROR TEST WITH ERTOL = ',D12.4/) WRITE (LUN,99994) 99994 FORMAT (' INPUT TO ZBESY ZN, FNU, KODE, N') ENDIF IF (KPRINT.GE.3) THEN WRITE (LUN,99993) 99993 FORMAT (' COMPARE Y(ZN,FNU) WITH COE1*I(Z,FNU)', * '-COE2*K(Z,FNU)') WRITE (LUN,99992) 99992 FORMAT (' Z = ZN*EXP(-i*PI/2)'/ * ' COE1 = EXP(i*(FNU+1)*PI/2) ', * ' COE2 = (2/PI)*EXP(-i*FNU*PI/2)') WRITE (LUN,99991) 99991 FORMAT (' RESULTS FROM ZBESY V(KK)'/ * 9X,'FROM ZBESI W(KK)'/ * 9X,'FROM ZBESK Y(KK)') WRITE (LUN,99990) 99990 FORMAT (' TEST CASE INDICES IR, IT, ICASE'/) ENDIF ENDIF LFLG = LFLG + 1 IF (KPRINT.GE.2) THEN WRITE (LUN,99989) ZNR, ZNI, FNU, KODE, N 99989 FORMAT (' INPUT: ZN=',2D12.4,3X,'FNU=',D12.4, * 3X,'KODE=',I3,3X,'N=',I3) ENDIF IF (KPRINT.GE.3) THEN WRITE (LUN,99988) (AER(K),K=1,N) 99988 FORMAT (' ERROR: AER(K)=',4D12.4) WRITE (LUN,99987) ZR, ZI, COE1R, COE1I, COE2R, * COE2R 99987 FORMAT (12X,'Z=',2D12.4/12X,'COE1=',2D12.4,3X, * 'COE2=',2D12.4) KK = MAX(NZ1,NZ2) + 1 KK = MIN(N,KK) WRITE (LUN,99986) VR(KK), VI(KK), WR(KK), WI(KK), * YR(KK), YI(KK) 99986 FORMAT (' RESULTS: V(KK)=',2D12.4/ * 12X,'W(KK)=',2D12.4/12X,'Y(KK)=',2D12.4) WRITE (LUN,99985) IR, IT, ICASE 99985 FORMAT (' CASE: IR=',I3,3X,'IT=',I3,3X, * 'ICASE=',I3/) ENDIF ENDIF 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE IF (KPRINT.GE.2) THEN IF (LFLG.EQ.0) THEN WRITE (LUN,99984) 99984 FORMAT (' QUICK CHECKS OK') ELSE WRITE (LUN,99983) LFLG 99983 FORMAT (' ***',I5,' FAILURE(S) FOR ZBESY IN THE (Z,FNU) ', * 'PLANE') ENDIF ENDIF IPASS = 0 IF (LFLG.EQ.0) THEN IPASS = 1 ENDIF IF (IPASS.EQ.1.AND.KPRINT.GE.2) THEN WRITE (LUN,99982) 99982 FORMAT (/' ****** ZBESY PASSED ALL TESTS ******'/) ENDIF IF (IPASS.EQ.0.AND.KPRINT.GE.1) THEN WRITE (LUN,99981) 99981 FORMAT (/' ****** ZBESY FAILED SOME TESTS ******'/) ENDIF RETURN END