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