Mercurial > octave-nkf

comparison libcruft/amos/crati.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 SUBROUTINE CRATI(Z, FNU, N, CY, TOL)
2 C***BEGIN PROLOGUE CRATI
3 C***REFER TO CBESI,CBESK,CBESH
4 C
5 C CRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD
6 C RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD
7 C RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B,
8 C MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973,
9 C BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER,
10 C BY D. J. SOOKNE.
11 C
12 C***ROUTINES CALLED (NONE)
13 C***END PROLOGUE CRATI
14 COMPLEX CDFNU, CONE, CY, CZERO, PT, P1, P2, RZ, T1, Z
15 REAL AK, AMAGZ, AP1, AP2, ARG, AZ, DFNU, FDNU, FLAM, FNU, FNUP,
16 * RAP1, RHO, TEST, TEST1, TOL
17 INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N
18 DIMENSION CY(N)
19 DATA CZERO, CONE / (0.0E0,0.0E0), (1.0E0,0.0E0) /
20 AZ = CABS(Z)
21 INU = INT(FNU)
22 IDNU = INU + N - 1
23 FDNU = FLOAT(IDNU)
24 MAGZ = INT(AZ)
25 AMAGZ = FLOAT(MAGZ+1)
26 FNUP = AMAX1(AMAGZ,FDNU)
27 ID = IDNU - MAGZ - 1
28 ITIME = 1
29 K = 1
30 RZ = (CONE+CONE)/Z
31 T1 = CMPLX(FNUP,0.0E0)*RZ
32 P2 = -T1
33 P1 = CONE
34 T1 = T1 + RZ
35 IF (ID.GT.0) ID = 0
36 AP2 = CABS(P2)
37 AP1 = CABS(P1)
38 C-----------------------------------------------------------------------
39 C THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNX
40 C GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT
41 C P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR
42 C PREMATURELY.
43 C-----------------------------------------------------------------------
44 ARG = (AP2+AP2)/(AP1*TOL)
45 TEST1 = SQRT(ARG)
46 TEST = TEST1
47 RAP1 = 1.0E0/AP1
48 P1 = P1*CMPLX(RAP1,0.0E0)
49 P2 = P2*CMPLX(RAP1,0.0E0)
50 AP2 = AP2*RAP1
51 10 CONTINUE
52 K = K + 1
53 AP1 = AP2
54 PT = P2
55 P2 = P1 - T1*P2
56 P1 = PT
57 T1 = T1 + RZ
58 AP2 = CABS(P2)
59 IF (AP1.LE.TEST) GO TO 10
60 IF (ITIME.EQ.2) GO TO 20
61 AK = CABS(T1)*0.5E0
62 FLAM = AK + SQRT(AK*AK-1.0E0)
63 RHO = AMIN1(AP2/AP1,FLAM)
64 TEST = TEST1*SQRT(RHO/(RHO*RHO-1.0E0))
65 ITIME = 2
66 GO TO 10
67 20 CONTINUE
68 KK = K + 1 - ID
69 AK = FLOAT(KK)
70 DFNU = FNU + FLOAT(N-1)
71 CDFNU = CMPLX(DFNU,0.0E0)
72 T1 = CMPLX(AK,0.0E0)
73 P1 = CMPLX(1.0E0/AP2,0.0E0)
74 P2 = CZERO
75 DO 30 I=1,KK
76 PT = P1
77 P1 = RZ*(CDFNU+T1)*P1 + P2
78 P2 = PT
79 T1 = T1 - CONE
80 30 CONTINUE
81 IF (REAL(P1).NE.0.0E0 .OR. AIMAG(P1).NE.0.0E0) GO TO 40
82 P1 = CMPLX(TOL,TOL)
83 40 CONTINUE
84 CY(N) = P2/P1
85 IF (N.EQ.1) RETURN
86 K = N - 1
87 AK = FLOAT(K)
88 T1 = CMPLX(AK,0.0E0)
89 CDFNU = CMPLX(FNU,0.0E0)*RZ
90 DO 60 I=2,N
91 PT = CDFNU + T1*RZ + CY(K+1)
92 IF (REAL(PT).NE.0.0E0 .OR. AIMAG(PT).NE.0.0E0) GO TO 50
93 PT = CMPLX(TOL,TOL)
94 50 CONTINUE
95 CY(K) = CONE/PT
96 T1 = T1 - CONE
97 K = K - 1
98 60 CONTINUE
99 RETURN
100 END