Mercurial > octave-nkf

comparison libcruft/amos/cuni1.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 CUNI1(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM,
2 * ALIM)
3 C***BEGIN PROLOGUE CUNI1
4 C***REFER TO CBESI,CBESK
5 C
6 C CUNI1 COMPUTES I(FNU,Z) BY MEANS OF THE UNIFORM ASYMPTOTIC
7 C EXPANSION FOR I(FNU,Z) IN -PI/3.LE.ARG Z.LE.PI/3.
8 C
9 C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
10 C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
11 C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
12 C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
13 C Y(I)=CZERO FOR I=NLAST+1,N
14 C
15 C***ROUTINES CALLED CUCHK,CUNIK,CUOIK,R1MACH
16 C***END PROLOGUE CUNI1
17 COMPLEX CFN, CONE, CRSC, CSCL, CSR, CSS, CWRK, CZERO, C1, C2,
18 * PHI, RZ, SUM, S1, S2, Y, Z, ZETA1, ZETA2, CY
19 REAL ALIM, APHI, ASCLE, BRY, C2I, C2M, C2R, ELIM, FN, FNU, FNUL,
20 * RS1, TOL, YY, R1MACH
21 INTEGER I, IFLAG, INIT, K, KODE, M, N, ND, NLAST, NN, NUF, NW, NZ
22 DIMENSION BRY(3), Y(N), CWRK(16), CSS(3), CSR(3), CY(2)
23 DATA CZERO, CONE / (0.0E0,0.0E0), (1.0E0,0.0E0) /
24 C
25 NZ = 0
26 ND = N
27 NLAST = 0
28 C-----------------------------------------------------------------------
29 C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
30 C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
31 C EXP(ALIM)=EXP(ELIM)*TOL
32 C-----------------------------------------------------------------------
33 CSCL = CMPLX(1.0E0/TOL,0.0E0)
34 CRSC = CMPLX(TOL,0.0E0)
35 CSS(1) = CSCL
36 CSS(2) = CONE
37 CSS(3) = CRSC
38 CSR(1) = CRSC
39 CSR(2) = CONE
40 CSR(3) = CSCL
41 BRY(1) = 1.0E+3*R1MACH(1)/TOL
42 C-----------------------------------------------------------------------
43 C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
44 C-----------------------------------------------------------------------
45 FN = AMAX1(FNU,1.0E0)
46 INIT = 0
47 CALL CUNIK(Z, FN, 1, 1, TOL, INIT, PHI, ZETA1, ZETA2, SUM, CWRK)
48 IF (KODE.EQ.1) GO TO 10
49 CFN = CMPLX(FN,0.0E0)
50 S1 = -ZETA1 + CFN*(CFN/(Z+ZETA2))
51 GO TO 20
52 10 CONTINUE
53 S1 = -ZETA1 + ZETA2
54 20 CONTINUE
55 RS1 = REAL(S1)
56 IF (ABS(RS1).GT.ELIM) GO TO 130
57 30 CONTINUE
58 NN = MIN0(2,ND)
59 DO 80 I=1,NN
60 FN = FNU + FLOAT(ND-I)
61 INIT = 0
62 CALL CUNIK(Z, FN, 1, 0, TOL, INIT, PHI, ZETA1, ZETA2, SUM, CWRK)
63 IF (KODE.EQ.1) GO TO 40
64 CFN = CMPLX(FN,0.0E0)
65 YY = AIMAG(Z)
66 S1 = -ZETA1 + CFN*(CFN/(Z+ZETA2)) + CMPLX(0.0E0,YY)
67 GO TO 50
68 40 CONTINUE
69 S1 = -ZETA1 + ZETA2
70 50 CONTINUE
71 C-----------------------------------------------------------------------
72 C TEST FOR UNDERFLOW AND OVERFLOW
73 C-----------------------------------------------------------------------
74 RS1 = REAL(S1)
75 IF (ABS(RS1).GT.ELIM) GO TO 110
76 IF (I.EQ.1) IFLAG = 2
77 IF (ABS(RS1).LT.ALIM) GO TO 60
78 C-----------------------------------------------------------------------
79 C REFINE TEST AND SCALE
80 C-----------------------------------------------------------------------
81 APHI = CABS(PHI)
82 RS1 = RS1 + ALOG(APHI)
83 IF (ABS(RS1).GT.ELIM) GO TO 110
84 IF (I.EQ.1) IFLAG = 1
85 IF (RS1.LT.0.0E0) GO TO 60
86 IF (I.EQ.1) IFLAG = 3
87 60 CONTINUE
88 C-----------------------------------------------------------------------
89 C SCALE S1 IF CABS(S1).LT.ASCLE
90 C-----------------------------------------------------------------------
91 S2 = PHI*SUM
92 C2R = REAL(S1)
93 C2I = AIMAG(S1)
94 C2M = EXP(C2R)*REAL(CSS(IFLAG))
95 S1 = CMPLX(C2M,0.0E0)*CMPLX(COS(C2I),SIN(C2I))
96 S2 = S2*S1
97 IF (IFLAG.NE.1) GO TO 70
98 CALL CUCHK(S2, NW, BRY(1), TOL)
99 IF (NW.NE.0) GO TO 110
100 70 CONTINUE
101 M = ND - I + 1
102 CY(I) = S2
103 Y(M) = S2*CSR(IFLAG)
104 80 CONTINUE
105 IF (ND.LE.2) GO TO 100
106 RZ = CMPLX(2.0E0,0.0E0)/Z
107 BRY(2) = 1.0E0/BRY(1)
108 BRY(3) = R1MACH(2)
109 S1 = CY(1)
110 S2 = CY(2)
111 C1 = CSR(IFLAG)
112 ASCLE = BRY(IFLAG)
113 K = ND - 2
114 FN = FLOAT(K)
115 DO 90 I=3,ND
116 C2 = S2
117 S2 = S1 + CMPLX(FNU+FN,0.0E0)*RZ*S2
118 S1 = C2
119 C2 = S2*C1
120 Y(K) = C2
121 K = K - 1
122 FN = FN - 1.0E0
123 IF (IFLAG.GE.3) GO TO 90
124 C2R = REAL(C2)
125 C2I = AIMAG(C2)
126 C2R = ABS(C2R)
127 C2I = ABS(C2I)
128 C2M = AMAX1(C2R,C2I)
129 IF (C2M.LE.ASCLE) GO TO 90
130 IFLAG = IFLAG + 1
131 ASCLE = BRY(IFLAG)
132 S1 = S1*C1
133 S2 = C2
134 S1 = S1*CSS(IFLAG)
135 S2 = S2*CSS(IFLAG)
136 C1 = CSR(IFLAG)
137 90 CONTINUE
138 100 CONTINUE
139 RETURN
140 C-----------------------------------------------------------------------
141 C SET UNDERFLOW AND UPDATE PARAMETERS
142 C-----------------------------------------------------------------------
143 110 CONTINUE
144 IF (RS1.GT.0.0E0) GO TO 120
145 Y(ND) = CZERO
146 NZ = NZ + 1
147 ND = ND - 1
148 IF (ND.EQ.0) GO TO 100
149 CALL CUOIK(Z, FNU, KODE, 1, ND, Y, NUF, TOL, ELIM, ALIM)
150 IF (NUF.LT.0) GO TO 120
151 ND = ND - NUF
152 NZ = NZ + NUF
153 IF (ND.EQ.0) GO TO 100
154 FN = FNU + FLOAT(ND-1)
155 IF (FN.GE.FNUL) GO TO 30
156 NLAST = ND
157 RETURN
158 120 CONTINUE
159 NZ = -1
160 RETURN
161 130 CONTINUE
162 IF (RS1.GT.0.0E0) GO TO 120
163 NZ = N
164 DO 140 I=1,N
165 Y(I) = CZERO
166 140 CONTINUE
167 RETURN
168 END