*DECK ZQCAI SUBROUTINE ZQCAI (LUN, KPRINT, IPASS) C***BEGIN PROLOGUE ZQCAI C***SUBSIDIARY C***PURPOSE Quick check for SLATEC subroutines C ZAIRY, ZBIRY C***LIBRARY SLATEC C***CATEGORY C10D C***TYPE COMPLEX (CQCAI-C, ZQCAI-Z) C***KEYWORDS QUICK CHECK, ZAIRY, ZBIRY 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 ZQCAI (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 ZQCAI is a quick check routine for the complex Airy functions C generated by subroutines ZAIRY and ZBIRY. C C ZQCAI generates Airy functions and their derivatives from ZAIRY C and ZBIRY and checks them against the Wronskian evaluation C in the Z plane. 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 ZAIRY, ZBIRY, ZABS, ZSQRT, ZEXP, I1MACH, D1MACH C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 890831 Revised to meet new SLATEC standards C 930122 Added ZEXP and ZSQRT to EXTERNAL statement. (RWC) C***END PROLOGUE ZQCAI 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, ZSQRT C C Declare local variables. C DOUBLE PRECISION CON1R,CON1I, CON2R,CON2I, CON3R,CON3I, CVR,CVI, * CWR,CWI, CYR,CYI, WR,WI, YR,YI, ZR,ZI, ZRR,ZRI DOUBLE PRECISION AA, AB, ACW, ACY, ALIM, ARG, ATOL, AV, AZRR, A1, * A2, CT, C23, DIG, ELIM, EPS, ER, ERTOL, FILM, FNUL, FPI, HPI, * PI, PI3, PTR, R, RL, RM, RPI, RTPI, R1M4, R1M5, SLAK, SPI, ST, * STI, STR, T, TOL, TPI, TPI3, TS INTEGER I, ICASE, ICL, IERR, IL, IR, IRB, IRSET, IT, ITL, K, KDO, * KEPS, KODE, K1, K2, LFLG, NZ1, NZ2, NZ3, NZ4 DIMENSION KDO(20), KEPS(20), T(20), WR(20), WI(20), YR(20), * YI(20) C C***FIRST EXECUTABLE STATEMENT ZQCAI IF (KPRINT.GE.2) THEN WRITE (LUN,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE AIRY FUNCTIONS FROM ', * 'ZAIRY AND ZBIRY'/) 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) SLAK = 3.0D0+4.0D0*(-LOG10(TOL)-7.0D0)/11.0D0 SLAK = MAX(SLAK,3.0D0) ERTOL = TOL*10.0D0**SLAK RL = 1.2D0*DIG + 3.0D0 RM = 0.5D0*(ALIM + ELIM) RM = MIN(RM,200.0D0) RM = MAX(RM,RL+10.0D0) FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) 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 (6D12.4/) ENDIF C----------------------------------------------------------------------- C Generate angles for construction of complex Z to be used in tests. C----------------------------------------------------------------------- FPI = ATAN(1.0D0) HPI = FPI + FPI PI = HPI + HPI TPI = PI + PI RPI = 1.0D0/PI TPI3 = TPI/3.0D0 SPI = PI/6.0D0 PI3 = SPI+SPI RTPI = 1.0D0/TPI A1 = RTPI*COS(SPI) A2 = RTPI*SIN(SPI) CON1R = COS(TPI3) CON1I = SIN(TPI3) CON2R = A1 CON2I = -A2 CON3R = RPI CON3I = 0.0D0 C23 = 2.0D0/3.0D0 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 ICL = 1 IL = 5 DO 5 I = 1,IL KDO(I) = 0 KEPS(I) = 0 5 CONTINUE ELSE ICL = 2 IL = 7 DO 6 I = 1,IL KDO(I) = 0 KEPS(I) = 0 6 CONTINUE KEPS(2) = 1 KEPS(3) = 1 KEPS(5) = 1 KEPS(6) = 1 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 ENDIF I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 C----------------------------------------------------------------------- C Test values of Z in -PI.lt.arg(Z).le.PI near formula boundaries. C----------------------------------------------------------------------- IF (KPRINT.GE.2) THEN WRITE (LUN,99996) 99996 FORMAT (' CHECKS IN THE Z PLANE'/) ENDIF LFLG = 0 DO 180 ICASE = 1,ICL C----------------------------------------------------------------------- C ICASE = 1 computes wron(AI(Z),BI(Z)) =CON3 C ICASE = 2 computes wron(AI(Z),AI(Z*CON1))=CON2 C----------------------------------------------------------------------- DO 170 KODE = 1,2 DO 160 IRSET = 1,3 IRB = MIN(IRSET,2) DO 150 IR = IRB,4 C------------ switch (irset) GO TO (40, 50, 60), IRSET 40 CONTINUE R = 2.0D0*(IR-1)/3.0D0 GO TO 70 50 CONTINUE R = (2.0D0*(4-IR)+RL*(IR-1))/3.0D0 GO TO 70 60 CONTINUE R = (RL*(4-IR)+RM*(IR-1))/3.0D0 70 CONTINUE C------------ end switch DO 140 IT = 1,ITL C---------------------------------------------------------------------- C The following values are set before the DO 30 loop: C C23 = 2/3 C CON1 = cmplx(cos(2PI/3),sin(2PI/3)) C CON2 = cmplx(cos(PI/6),-sin(PI/6)/2PI C CON3 = cmplx(1/PI,0) C---------------------------------------------------------------------- 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 ZSQRT(ZR, ZI, STR, STI) PTR = (ZR*STR-ZI*STI)*C23 ZRI = (ZR*STI+ZI*STR)*C23 ZRR = PTR AZRR = ABS(ZRR) C-------------- Check for possible underflow or overflow IF (AZRR.NE.0.0D0) THEN ARG = -AZRR - 0.5D0*LOG(AZRR) + 0.226D0 ARG = ARG + ARG C---------------- Skip test for this case? IF (ARG.LT.(-ELIM)) GO TO 140 ENDIF CALL ZAIRY(ZR, ZI, 0, KODE, YR(1), YI(1), NZ1, IERR) CALL ZAIRY(ZR, ZI, 1, KODE, YR(2), YI(2), NZ2, IERR) IF (ICASE.EQ.1) THEN C---------------- Compare 1/PI with Wronskian of ZAIRY(Z) and ZBIRY(Z). CALL ZBIRY(ZR, ZI, 0, KODE, WR(1), WI(1), IERR) CALL ZBIRY(ZR, ZI, 1, KODE, WR(2), WI(2), IERR) IF (KODE.EQ.2) THEN C----------------------------------------------------------------------- C When KODE = 2, the scaling factor exp(-zeta1-zeta2) is 1.0 for C -PI.lt.arg(Z).le.PI/3 and exp(-2.0*zeta1) for PI/3.lt.arg(Z) C .le.PI where zeta1 = zeta2 in this range. This is due to the fact C that arg(Z*CON1) is taken to be in (-PI,PI) by the principal C square root. C----------------------------------------------------------------------- C------------------ Adjust scaling factor. CVR = AZRR - ZRR CVI = -ZRI CALL ZEXP(CVR, CVI, CVR, CVI) STR = WR(1)*CVR - WI(1)*CVI WI(1) = WR(1)*CVI + WI(1)*CVR WR(1) = STR STR = WR(2)*CVR - WI(2)*CVI WI(2) = WR(2)*CVI + WI(2)*CVR WR(2) = STR ENDIF CVR = CON3R CVI = CON3I ELSE C---------------- Compare exp(-i*PI/6)/2PI with Wronskian of ZAIRY(Z) C and ZAIRY(Z*exp(2i*PI/3)). CVR = ZR*CON1R - ZI*CON1I CVI = ZR*CON1I + ZI*CON1R CALL ZAIRY(CVR, CVI, 0, KODE, WR(1), WI(1), NZ3, IERR) CALL ZAIRY(CVR, CVI, 1, KODE, WR(2), WI(2), NZ4, IERR) IF (KODE.EQ.2) THEN IF (T(IT).GE.PI3) THEN C-------------------- Adjust scaling factor. CVR = ZRR + ZRR CVI = ZRI + ZRI CALL ZEXP(-CVR, -CVI, CVR, CVI) STR = WR(1)*CVR - WI(1)*CVI WI(1) = WR(1)*CVI + WI(1)*CVR WR(1) = STR STR = WR(2)*CVR - WI(2)*CVI WI(2) = WR(2)*CVI + WI(2)*CVR WR(2) = STR ENDIF ENDIF STR = WR(2)*CON1R - WI(2)*CON1I WI(2) = WR(2)*CON1I + WI(2)*CON1R WR(2) = STR CVR = CON2R CVI = CON2I ENDIF C----------------------------------------------------------------------- C Error relative to maximum term C----------------------------------------------------------------------- AV = ZABS(CVR,CVI) CWR = YR(1)*WR(2) - YI(1)*WI(2) CWI = YR(1)*WI(2) + YI(1)*WR(2) CYR = YR(2)*WR(1) - YI(2)*WI(1) CYI = YR(2)*WI(1) + YI(2)*WR(1) CYR = CWR - CYR - CVR CYI = CWI - CYI - CVI ACY = ZABS(YR(1),YI(1))*ZABS(WR(2),WI(2)) ACW = ZABS(WR(1),WI(1))*ZABS(YR(2),YI(2)) AV = MAX(ACW,ACY,AV) ER = ZABS(CYR,CYI)/AV IF (ER.GE.ERTOL) 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 ZAIRY AND ERROR') ENDIF IF (KPRINT.GE.3) THEN WRITE (LUN,99993) 99993 FORMAT (' COMPARISON VALUE AND WRONSKIAN') WRITE (LUN,99992) 99992 FORMAT (' RESULTS FROM ZAIRY AND/OR ZBIRY') WRITE (LUN,99991) 99991 FORMAT (' TEST CASE INDICES'/) ENDIF ENDIF LFLG = 1 IF (KPRINT.GE.2) THEN WRITE (LUN,99990) ZR, ZI, ER 99990 FORMAT (12X,'INPUT: Z=',2D12.4,5X,'ERROR: ER=', * D12.4) ENDIF IF (KPRINT.GE.3) THEN WRITE (LUN,99989) CVR, CVI, CYR, CYI 99989 FORMAT (' COMPARISON VALUE: CV=',2D12.4/ * 8X,'WRONSKIAN: CY=',2D12.4) WRITE (LUN,99988) NZ1, YR(1), YI(1), * NZ2, YR(2), YI(2) 99988 FORMAT (10X,'RESULTS: NZ1=',I3,4X,'Y(1)=',2D12.4/ * 20X,'NZ2=',I3,4X,'Y(2)=',2D12.4) IF (ICASE.EQ.1) THEN WRITE (LUN,99987) WR(1), WI(1), WR(2), WI(2) 99987 FORMAT (31X,'W(1)=',2D12.4/31X,'W(2)=',2D12.4) ELSE WRITE (LUN,99986) NZ3, WR(1), WI(1), * NZ4, WR(2), WI(2) 99986 FORMAT (20X,'NZ3=',I3,4X,'W(1)=',2D12.4/ * 20X,'NZ4=',I3,4X,'W(2)=',2D12.4) ENDIF WRITE (LUN,99985) IT, IR, IRSET, ICASE 99985 FORMAT (13X,'CASE: IT=',I3,4X,'IR=',I3,4X, * 'IRSET=',I3,4X,'ICASE=',I3,4X/) ENDIF ENDIF 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE IF (KPRINT.GE.2) THEN IF (LFLG.EQ.0) THEN WRITE (LUN,99984) 99984 FORMAT (' QUICK CHECKS OK') ELSE WRITE (LUN,99983) 99983 FORMAT (' ***',5X,'FAILURE(S) FOR ZAIRY IN THE Z 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 (/' ****** ZAIRY PASSED ALL TESTS ******'/) ENDIF IF (IPASS.EQ.0.AND.KPRINT.GE.1) THEN WRITE (LUN,99981) 99981 FORMAT (/' ****** ZAIRY FAILED SOME TESTS ******'/) ENDIF RETURN END