Mercurial > octave-nkf

comparison libcruft/amos/cbesy.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 CBESY(Z, FNU, KODE, N, CY, NZ, CWRK, IERR)
2 C***BEGIN PROLOGUE CBESY
3 C***DATE WRITTEN 830501 (YYMMDD)
4 C***REVISION DATE 890801 (YYMMDD)
5 C***CATEGORY NO. B5K
6 C***KEYWORDS Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT,
7 C BESSEL FUNCTION OF SECOND KIND
8 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
9 C***PURPOSE TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT
10 C***DESCRIPTION
11 C
12 C ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
13 C BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE
14 C ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
15 C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESY RETURNS THE SCALED
16 C FUNCTIONS
17 C
18 C CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z) I = 1,...,N , Y=AIMAG(Z)
19 C
20 C WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND
21 C LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
22 C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
23 C (REF. 1).
24 C
25 C INPUT
26 C Z - Z=CMPLX(X,Y), Z.NE.CMPLX(0.,0.),-PI.LT.ARG(Z).LE.PI
27 C FNU - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0E0
28 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
29 C KODE= 1 RETURNS
30 C CY(I)=Y(FNU+I-1,Z), I=1,...,N
31 C = 2 RETURNS
32 C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N
33 C WHERE Y=AIMAG(Z)
34 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
35 C CWRK - A COMPLEX WORK VECTOR OF DIMENSION AT LEAST N
36 C
37 C OUTPUT
38 C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
39 C VALUES FOR THE SEQUENCE
40 C CY(I)=Y(FNU+I-1,Z) OR
41 C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)) I=1,...,N
42 C DEPENDING ON KODE.
43 C NZ - NZ=0 , A NORMAL RETURN
44 C NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO
45 C UNDERFLOW (GENERALLY ON KODE=2)
46 C IERR - ERROR FLAG
47 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
48 C IERR=1, INPUT ERROR - NO COMPUTATION
49 C IERR=2, OVERFLOW - NO COMPUTATION, FNU+N-1 IS
50 C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
51 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
52 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
53 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE
54 C ACCURACY
55 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
56 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
57 C CANCE BY ARGUMENT REDUCTION
58 C IERR=5, ERROR - NO COMPUTATION,
59 C ALGORITHM TERMINATION CONDITION NOT MET
60 C
61 C***LONG DESCRIPTION
62 C
63 C THE COMPUTATION IS CARRIED OUT BY THE FORMULA
64 C
65 C Y(FNU,Z)=0.5*(H(1,FNU,Z)-H(2,FNU,Z))/I
66 C
67 C WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z)
68 C AND H(2,FNU,Z) ARE CALCULATED IN CBESH.
69 C
70 C FOR NEGATIVE ORDERS,THE FORMULA
71 C
72 C Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU)
73 C
74 C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD
75 C INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE
76 C POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)*
77 C SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS
78 C NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A
79 C LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM
80 C CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS,
81 C WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF
82 C ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z).
83 C
84 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
85 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
86 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
87 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
88 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
89 C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO
90 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
91 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
92 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
93 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
94 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
95 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
96 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
97 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
98 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
99 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
100 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
101 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
102 C
103 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
104 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
105 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
106 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
107 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
108 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
109 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
110 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
111 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
112 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
113 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
114 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
115 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
116 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
117 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
118 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
119 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
120 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
121 C OR -PI/2+P.
122 C
123 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
124 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
125 C COMMERCE, 1955.
126 C
127 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
128 C BY D. E. AMOS, SAND83-0083, MAY, 1983.
129 C
130 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
131 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
132 C
133 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
134 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
135 C 1018, MAY, 1985
136 C
137 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
138 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
139 C MATH. SOFTWARE, 1986
140 C
141 C***ROUTINES CALLED CBESH,I1MACH,R1MACH
142 C***END PROLOGUE CBESY
143 C
144 COMPLEX CWRK, CY, C1, C2, EX, HCI, Z, ZU, ZV
145 REAL ELIM, EY, FNU, R1, R2, TAY, XX, YY, R1MACH, ASCLE, RTOL,
146 * ATOL, AA, BB
147 INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH
148 DIMENSION CY(N), CWRK(N)
149 C***FIRST EXECUTABLE STATEMENT CBESY
150 XX = REAL(Z)
151 YY = AIMAG(Z)
152 IERR = 0
153 NZ=0
154 IF (XX.EQ.0.0E0 .AND. YY.EQ.0.0E0) IERR=1
155 IF (FNU.LT.0.0E0) IERR=1
156 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
157 IF (N.LT.1) IERR=1
158 IF (IERR.NE.0) RETURN
159 HCI = CMPLX(0.0E0,0.5E0)
160 CALL CBESH(Z, FNU, KODE, 1, N, CY, NZ1, IERR)
161 IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
162 CALL CBESH(Z, FNU, KODE, 2, N, CWRK, NZ2, IERR)
163 IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
164 NZ = MIN0(NZ1,NZ2)
165 IF (KODE.EQ.2) GO TO 60
166 DO 50 I=1,N
167 CY(I) = HCI*(CWRK(I)-CY(I))
168 50 CONTINUE
169 RETURN
170 60 CONTINUE
171 TOL = AMAX1(R1MACH(4),1.0E-18)
172 K1 = I1MACH(12)
173 K2 = I1MACH(13)
174 K = MIN0(IABS(K1),IABS(K2))
175 R1M5 = R1MACH(5)
176 C-----------------------------------------------------------------------
177 C ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT
178 C-----------------------------------------------------------------------
179 ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
180 R1 = COS(XX)
181 R2 = SIN(XX)
182 EX = CMPLX(R1,R2)
183 EY = 0.0E0
184 TAY = ABS(YY+YY)
185 IF (TAY.LT.ELIM) EY = EXP(-TAY)
186 IF (YY.LT.0.0E0) GO TO 90
187 C1 = EX*CMPLX(EY,0.0E0)
188 C2 = CONJG(EX)
189 70 CONTINUE
190 NZ = 0
191 RTOL = 1.0E0/TOL
192 ASCLE = R1MACH(1)*RTOL*1.0E+3
193 DO 80 I=1,N
194 C CY(I) = HCI*(C2*CWRK(I)-C1*CY(I))
195 ZV = CWRK(I)
196 AA=REAL(ZV)
197 BB=AIMAG(ZV)
198 ATOL=1.0E0
199 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 75
200 ZV = ZV*CMPLX(RTOL,0.0E0)
201 ATOL = TOL
202 75 CONTINUE
203 ZV = ZV*C2*HCI
204 ZV = ZV*CMPLX(ATOL,0.0E0)
205 ZU=CY(I)
206 AA=REAL(ZU)
207 BB=AIMAG(ZU)
208 ATOL=1.0E0
209 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 85
210 ZU = ZU*CMPLX(RTOL,0.0E0)
211 ATOL = TOL
212 85 CONTINUE
213 ZU = ZU*C1*HCI
214 ZU = ZU*CMPLX(ATOL,0.0E0)
215 CY(I) = ZV - ZU
216 IF (CY(I).EQ.CMPLX(0.0E0,0.0E0) .AND. EY.EQ.0.0E0) NZ = NZ + 1
217 80 CONTINUE
218 RETURN
219 90 CONTINUE
220 C1 = EX
221 C2 = CONJG(EX)*CMPLX(EY,0.0E0)
222 GO TO 70
223 170 CONTINUE
224 NZ = 0
225 RETURN
226 END