*DECK CQCBI SUBROUTINE CQCBI (LUN, KPRINT, IPASS) C***BEGIN PROLOGUE CQCBI C***SUBSIDIARY C***PURPOSE Quick check for SLATEC subroutine C CBESI C***LIBRARY SLATEC C***CATEGORY C10B4 C***TYPE COMPLEX (CQCBI-C, ZQCBI-Z) C***KEYWORDS QUICK CHECK, CBESI 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 CQCBI (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 CQCBI is a quick check routine for the complex I Bessel function C generated by subroutine CBESI. C C CQCBI generates sequences crossing formula boundaries which C are started by one formula and terminated in a region where C another formula applies. The terminated value is checked by C the formula appropriate to that region. 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 CBESI, CBESK, CWRSK, I1MACH, R1MACH C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 890831 Revised to meet new SLATEC standards C C***END PROLOGUE CQCBI C C*Internal Notes: C Machine constants are defined by functions I1MACH and R1MACH. 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 REAL R1MACH EXTERNAL I1MACH, R1MACH C C Declare local variables. C COMPLEX CK, CONE, CSGN, CW, CY, W, Y, Z, ZN, ZSC, ZT, ZZ REAL AA, AB, AER, ALIM, ARG, ATOL, AW, CARG, CT, DIG, ELIM, EPS, * ER, ERTOL, FILM, FNU, FNUL, GNU, HPI, PI, R, RL, RLT, RM, R1, * R1M4, R1M5, R2, SARG, SLAK, ST, T, TOL, TS, XX, YY INTEGER I, ICASE, IERR, IL, INU, IPRNT, IR, IT, ITL, K, KDO, * KEPS, KK, KODE, K1, K2, LFLG, MFLG, N, NL, NZI, NZK, NZ1, NZ2, * N1 DIMENSION AER(20), CK(2), KDO(20), KEPS(20), T(20), W(20), Y(20) C C***FIRST EXECUTABLE STATEMENT CQCBI IF (KPRINT.GE.2) THEN WRITE (LUN,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE I BESSEL FUNCTION FROM ', * 'CBESI'/) ENDIF C----------------------------------------------------------------------- C Set parameters related to machine constants. C TOL is the approximate unit roundoff limited to 1.0E-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 = R1MACH(4) TOL = MAX(R1M4,1.0E-18) ATOL = 100.0E0*TOL AA = -LOG10(R1M4) K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN(ABS(K1),ABS(K2)) ELIM = 2.303E0*(K*R1M5-3.0E0) AB = AA*2.303E0 ALIM = ELIM + MAX(-AB,-41.45E0) DIG = MIN(AA,18.0E0) FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) RL = 1.2E0*DIG + 3.0E0 SLAK = 3.0E0+4.0E0*(-LOG10(TOL)-7.0E0)/11.0E0 SLAK = MAX(SLAK,3.0E0) ERTOL = TOL*10.0E0**SLAK RM = 0.5E0*(ALIM + ELIM) RM = MIN(RM,200.0E0) RM = MAX(RM,RL+10.0E0) R2 = MIN(RM,FNUL) R1 = 2.0E0*SQRT(FNUL+1.0E0) 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,6E12.4/) ENDIF C----------------------------------------------------------------------- C Set other constants needed in the tests. C----------------------------------------------------------------------- ZZ = CMPLX(1.0E0/TOL,0.0E0) CONE = CMPLX(1.0E0,0.0E0) HPI = 2.0E0*ATAN(1.0E0) 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 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 ENDIF I = 2 EPS = 0.01E0 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.lt.arg(Z).le.PI near formula boundaries. C----------------------------------------------------------------------- IF (KPRINT.GE.2) THEN WRITE (LUN,99996) 99996 FORMAT (' CHECKS ACROSS FORMULA BOUNDARIES') ENDIF LFLG = 0 DO 220 ICASE = 1,6 DO 210 KODE = 1,2 DO 200 N = 1,NL N1 = N + 2 C----------------------------------------------------------------------- C Set values for R = magnitude of Z and for FNU to test computation C methods for the various regions of the (Z,FNU) plane. C----------------------------------------------------------------------- DO 190 IR = 1,3 C------------ switch (icase) GO TO (50, 60, 70, 80, 90, 100), ICASE 50 CONTINUE R = (2.0E0*(3-IR)+RL*(IR-1))/2.0E0 GNU = R*R/4.0E0 - 0.2E0 - (N-1) FNU = MAX(0.0E0,GNU) GO TO 110 60 CONTINUE R = (RL*(3-IR)+R2*(IR-1))/2.0E0 GNU = SQRT(R+R) - 0.2E0 - (N-1) FNU = MAX(0.0E0,GNU) GO TO 110 70 CONTINUE IF (R2.GE.RM) GO TO 220 R = (R2*(3-IR)+RM*(IR-1))/2.0E0 GNU = SQRT(R+R) - 0.2E0 - (N-1) FNU = MAX(0.0E0,GNU) GO TO 110 80 CONTINUE IF (R1.GE.RL) GO TO 220 R = (R1*(3-IR)+RL*(IR-1))/2.0E0 FNU = FNUL - 0.2E0 - (N-1) GO TO 110 90 CONTINUE R = (RL*(3-IR)+R2*(IR-1))/2.0E0 FNU = FNUL - 0.2E0 - (N-1) GO TO 110 100 CONTINUE IF (R2.GE.RM) GO TO 220 R = (R2*(3-IR)+RM*(IR-1))/2.0E0 FNU = FNUL - 0.2E0 - (N-1) 110 CONTINUE C------------ end switch DO 180 IT = 1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (ABS(CT).LT.ATOL) CT = 0.0E0 IF (ABS(ST).LT.ATOL) ST = 0.0E0 Z = CMPLX(R*CT,R*ST) XX = REAL(Z) YY = AIMAG(Z) CALL CBESI(Z, FNU, KODE, N1, Y, NZ1, IERR) IF (NZ1.NE.0) GO TO 180 C----------------------------------------------------------------------- C Compare values from CBESI with values from CWRSK, an alternative C method for calculating the complex Bessel I function. C----------------------------------------------------------------------- ZN = Z IF (XX.GE.0.0E0) THEN CALL CWRSK(ZN, FNU, KODE, N, W, NZ2, CK, TOL, ELIM, * ALIM) IF (NZ2.NE.0) GO TO 180 ELSE ZN = -Z INU = INT(FNU) ARG = (FNU-INU)*PI IF (YY.LT.0.0E0) ARG = -ARG CARG = COS(ARG) SARG = SIN(ARG) CSGN = CMPLX(CARG,SARG) IF (MOD(INU,2).EQ.1) CSGN = -CSGN CALL CWRSK(ZN, FNU, KODE, N, W, NZ2, CK, TOL, ELIM, * ALIM) IF (NZ2.NE.0) GO TO 180 DO 130 I = 1,N W(I) = W(I)*CSGN CSGN = -CSGN 130 CONTINUE ENDIF MFLG = 0 DO 160 I = 1,N AB = FNU + I - 1 AA = MAX(2.0E0,AB) ZT = W(I) IF (CABS(ZT).GT.1.0E0) THEN ZSC = CMPLX(TOL,0.0E0) ELSE ZSC = ZZ C------------------ ZZ = CMPLX(1.0/TOL,0.0) ENDIF CW = W(I)*ZSC CY = Y(I)*ZSC ER = CABS(CY-CW) AW = CABS(CW) IF (AW.NE.0.0E0) THEN IF (XX.EQ.0.0E0) THEN IF (ABS(YY).LT.AA) THEN ER=ER/AW ELSE ER=CABS(W(I)-Y(I)) ENDIF ELSE ER=ER/AW ENDIF ELSE ER=CABS(Y(I)) ENDIF AER(I) = ER IF (ER.GE.ERTOL) MFLG = 1 160 CONTINUE C----------------------------------------------------------------------- C Write failure reports for KPRINT.ge.2 and KPRINT.ge.3 C----------------------------------------------------------------------- IF (MFLG.NE.0) THEN IF (LFLG.EQ.0) THEN IF (KPRINT.GE.2) THEN WRITE (LUN,99995) ERTOL 99995 FORMAT (/' CASES WHICH UNDERFLOW OR VIOLATE THE ', * 'RELATIVE ERROR TEST'/' WITH ERTOL = ', E12.4/) WRITE (LUN,99994) 99994 FORMAT (' INPUT TO CBESI Z, FNU, KODE, N') ENDIF IF (KPRINT.GE.3) THEN WRITE (LUN,99993) 99993 FORMAT (' ERROR TEST ON RESULTS FROM CBESI AND ', * 'CWRSK AER(K)') WRITE (LUN,99992) 99992 FORMAT (' RESULTS FROM CBESI NZ1, Y(KK)'/, * ' RESULTS FROM CWRSK NZ2, W(KK)') WRITE (LUN,99991) 99991 FORMAT (' TEST CASE INDICES IT, IR, ICASE'/) ENDIF ENDIF LFLG = LFLG + 1 IF (KPRINT.GE.2) THEN WRITE (LUN,99990) Z, FNU, KODE, N 99990 FORMAT (' INPUT: Z=',2E12.4,4X,'FNU=',E12.4,4X, * 'KODE=',I3,4X,'N=',I3) ENDIF IF (KPRINT.GE.3) THEN WRITE (LUN,99989) (AER(K),K=1,N) 99989 FORMAT (' ERROR: AER(K)=',4E12.4) KK = MAX(NZ1,NZ2) + 1 KK = MIN(N,KK) WRITE (LUN,99988) NZ1, Y(KK), NZ2, W(KK) 99988 FORMAT (' RESULTS: NZ1=',I3,4X,'Y(KK)=',2E12.4, * /11X,'NZ2=',I3,4X,'W(KK)=',2E12.4) WRITE (LUN,99987) IT, IR, ICASE 99987 FORMAT (' CASE: IT=',I3,4X,'IR=',I3,4X, * 'ICASE=',I3/) ENDIF ENDIF 180 CONTINUE 190 CONTINUE 200 CONTINUE 210 CONTINUE 220 CONTINUE IF (KPRINT.GE.2) THEN IF (LFLG.EQ.0) THEN WRITE (LUN,99986) 99986 FORMAT (' QUICK CHECKS OK') ELSE WRITE (LUN,99985) LFLG 99985 FORMAT(' ***',I5,' FAILURE(S) FOR CBESI CHECKS NEAR FORMULA ', * 'BOUNDARIES') ENDIF ENDIF C C IPRNT = 0 IF (MQC.EQ.1) GO TO 900 C----------------------------------------------------------------------- C Checks near underflow limits on series(I=1) and uniform C asymptotic expansion(I=2) C Compare 1/Z with I(Z,FNU)*K(Z,FNU+1) + I(Z,FNU+1)*K(Z,FNU) and C report cases for which the relative error is greater than ERTOL. C----------------------------------------------------------------------- IF (KPRINT.GE.2) THEN WRITE (LUN,99984) 99984 FORMAT (/' CHECKS NEAR UNDERFLOW AND OVERFLOW LIMITS'/) ENDIF Z = (1.4E0,1.4E0) KODE = 1 N = 20 DO 280 I = 1,2 FNU = 10.2E0 C----------------------------------------------------------------------- C Adjust FNU by repeating until 0.lt.NZI.lt.10 C----------------------------------------------------------------------- 230 CONTINUE CALL CBESI(Z, FNU, KODE, N, Y, NZI, IERR) IF (NZI.NE.0) THEN IF (NZI.GE.10) THEN FNU = FNU - 1.0E0 GO TO 230 ENDIF ELSE FNU = FNU + 5.0E0 GO TO 230 ENDIF C------ End repeat CALL CBESK(Z, FNU, KODE, 2, W, NZK, IERR) ZT = CONE/Z CY = W(1)*Y(2) CW = W(2)*Y(1) CW = CW + CY - ZT ER = ABS(CW)/ABS(ZT) C----------------------------------------------------------------------- C Write failure reports for KPRINT.ge.2 and KPRINT.ge.3 C----------------------------------------------------------------------- IF (ER.GE.ERTOL) THEN IF (IPRNT.EQ.0) THEN IF (KPRINT.GE.2) THEN WRITE (LUN,99983) 99983 FORMAT (' INPUT TO CBESI Z, FNU, KODE, N') ENDIF IF (KPRINT.GE.3) THEN WRITE (LUN,99982) 99982 FORMAT (' COMPARE 1/Z WITH WRONSKIAN(CBESI(Z,FNU),', * 'CBESK(Z,FNU))'/) ENDIF ENDIF IPRNT = 1 IF (KPRINT.GE.2) THEN WRITE (LUN,99981) Z, FNU, KODE, N 99981 FORMAT (' INPUT: Z=',2E12.4,3X,'FNU=',E12.4,3X,'KODE=',I3, * 3X,'N=',I3) ENDIF IF (KPRINT.GE.3) THEN WRITE (LUN,99980) ZT, CW+CY 99980 FORMAT (' RESULTS:',15X,'1/Z=',2E12.4/ * 10X,'WRON(CBESI,CBESK)=',2E12.4) WRITE (LUN,99979) ER 99979 FORMAT (' RELATIVE ERROR:',9X,'ER=',E12.4/) ENDIF ENDIF RLT = RL + RL Z = CMPLX(RLT,0.0E0) 280 CONTINUE C----------------------------------------------------------------------- C Check near overflow limits C Compare 1/Z with I(Z,FNU)*K(Z,FNU+1) + I(Z,FNU+1)*K(Z,FNU) and C report cases for which the relative error is greater than ERTOL. C----------------------------------------------------------------------- Z = CMPLX(ELIM,0.0E0) FNU = 0.0E0 C----------------------------------------------------------------------- C Adjust FNU by repeating until NZK.lt.10 C N = 20 set before DO 280 loop C----------------------------------------------------------------------- 290 CONTINUE CALL CBESK(Z, FNU, KODE, N, Y, NZK, IERR) IF (NZK.GE.10) THEN IF (NZK.EQ.N) THEN FNU = FNU + 3.0E0 ELSE FNU = FNU + 2.0E0 ENDIF GO TO 290 ENDIF C---- End repeat GNU = FNU + (N-2) CALL CBESI(Z, GNU, KODE, 2, W, NZI, IERR) ZT = CONE/Z CY = Y(N-1)*W(2) CW = Y(N)*W(1) CW = CW + CY - ZT ER = ABS(CW)/ABS(ZT) IF (ER.GE.ERTOL) THEN IF (IPRNT.EQ.0) THEN IF (KPRINT.GE.2) THEN WRITE (LUN,99983) ENDIF IF (KPRINT.GE.3) THEN WRITE (LUN,99982) ENDIF ENDIF IPRNT = 1 IF (KPRINT.GE.2) THEN WRITE (LUN,99981) Z, FNU, KODE, N ENDIF IF (KPRINT.GE.3) THEN WRITE (LUN,99980) ZT, CW+CY WRITE (LUN,99979) ER ENDIF ENDIF IF (KPRINT.GE.2) THEN IF (IPRNT.EQ.0) THEN WRITE (LUN,99986) C 99986 FORMAT (' QUICK CHECKS OK') ELSE WRITE (LUN,99978) 99978 FORMAT (' ***',5X,'FAILURE(S) FOR CBESI NEAR UNDERFLOW AND ', * 'OVERFLOW LIMITS') ENDIF ENDIF 900 CONTINUE IPASS = 0 IF (IPRNT.EQ.0.AND.LFLG.EQ.0) THEN IPASS = 1 ENDIF IF (IPASS.EQ.1.AND.KPRINT.GE.2) THEN WRITE (LUN,99977) 99977 FORMAT (/' ****** CBESI PASSED ALL TESTS ******'/) ENDIF IF (IPASS.EQ.0.AND.KPRINT.GE.1) THEN WRITE (LUN,99976) 99976 FORMAT (/' ****** CBESI FAILED SOME TESTS ******'/) ENDIF RETURN END