Mercurial > octave

comparison libcruft/lapack/slazq3.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 SLAZQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
2 $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
3 $ DN2, TAU )
4 *
5 * -- LAPACK auxiliary routine (version 3.1) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 LOGICAL IEEE
11 INTEGER I0, ITER, N0, NDIV, NFAIL, PP, TTYPE
12 REAL DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, QMAX,
13 $ SIGMA, TAU
14 * ..
15 * .. Array Arguments ..
16 REAL Z( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * SLAZQ3 checks for deflation, computes a shift (TAU) and calls dqds.
23 * In case of failure it changes shifts, and tries again until output
24 * is positive.
25 *
26 * Arguments
27 * =========
28 *
29 * I0 (input) INTEGER
30 * First index.
31 *
32 * N0 (input) INTEGER
33 * Last index.
34 *
35 * Z (input) REAL array, dimension ( 4*N )
36 * Z holds the qd array.
37 *
38 * PP (input) INTEGER
39 * PP=0 for ping, PP=1 for pong.
40 *
41 * DMIN (output) REAL
42 * Minimum value of d.
43 *
44 * SIGMA (output) REAL
45 * Sum of shifts used in current segment.
46 *
47 * DESIG (input/output) REAL
48 * Lower order part of SIGMA
49 *
50 * QMAX (input) REAL
51 * Maximum value of q.
52 *
53 * NFAIL (output) INTEGER
54 * Number of times shift was too big.
55 *
56 * ITER (output) INTEGER
57 * Number of iterations.
58 *
59 * NDIV (output) INTEGER
60 * Number of divisions.
61 *
62 * IEEE (input) LOGICAL
63 * Flag for IEEE or non IEEE arithmetic (passed to SLASQ5).
64 *
65 * TTYPE (input/output) INTEGER
66 * Shift type. TTYPE is passed as an argument in order to save
67 * its value between calls to SLAZQ3
68 *
69 * DMIN1 (input/output) REAL
70 * DMIN2 (input/output) REAL
71 * DN (input/output) REAL
72 * DN1 (input/output) REAL
73 * DN2 (input/output) REAL
74 * TAU (input/output) REAL
75 * These are passed as arguments in order to save their values
76 * between calls to SLAZQ3
77 *
78 * This is a thread safe version of SLASQ3, which passes TTYPE, DMIN1,
79 * DMIN2, DN, DN1. DN2 and TAU through the argument list in place of
80 * declaring them in a SAVE statment.
81 *
82 * =====================================================================
83 *
84 * .. Parameters ..
85 REAL CBIAS
86 PARAMETER ( CBIAS = 1.50E0 )
87 REAL ZERO, QURTR, HALF, ONE, TWO, HUNDRD
88 PARAMETER ( ZERO = 0.0E0, QURTR = 0.250E0, HALF = 0.5E0,
89 $ ONE = 1.0E0, TWO = 2.0E0, HUNDRD = 100.0E0 )
90 * ..
91 * .. Local Scalars ..
92 INTEGER IPN4, J4, N0IN, NN
93 REAL EPS, G, S, SAFMIN, T, TEMP, TOL, TOL2
94 * ..
95 * .. External Subroutines ..
96 EXTERNAL SLASQ5, SLASQ6, SLAZQ4
97 * ..
98 * .. External Function ..
99 REAL SLAMCH
100 EXTERNAL SLAMCH
101 * ..
102 * .. Intrinsic Functions ..
103 INTRINSIC ABS, MIN, SQRT
104 * ..
105 * .. Executable Statements ..
106 *
107 N0IN = N0
108 EPS = SLAMCH( 'Precision' )
109 SAFMIN = SLAMCH( 'Safe minimum' )
110 TOL = EPS*HUNDRD
111 TOL2 = TOL**2
112 G = ZERO
113 *
114 * Check for deflation.
115 *
116 10 CONTINUE
117 *
118 IF( N0.LT.I0 )
119 $ RETURN
120 IF( N0.EQ.I0 )
121 $ GO TO 20
122 NN = 4*N0 + PP
123 IF( N0.EQ.( I0+1 ) )
124 $ GO TO 40
125 *
126 * Check whether E(N0-1) is negligible, 1 eigenvalue.
127 *
128 IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
129 $ Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
130 $ GO TO 30
131 *
132 20 CONTINUE
133 *
134 Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
135 N0 = N0 - 1
136 GO TO 10
137 *
138 * Check whether E(N0-2) is negligible, 2 eigenvalues.
139 *
140 30 CONTINUE
141 *
142 IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
143 $ Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
144 $ GO TO 50
145 *
146 40 CONTINUE
147 *
148 IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
149 S = Z( NN-3 )
150 Z( NN-3 ) = Z( NN-7 )
151 Z( NN-7 ) = S
152 END IF
153 IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN
154 T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
155 S = Z( NN-3 )*( Z( NN-5 ) / T )
156 IF( S.LE.T ) THEN
157 S = Z( NN-3 )*( Z( NN-5 ) /
158 $ ( T*( ONE+SQRT( ONE+S / T ) ) ) )
159 ELSE
160 S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
161 END IF
162 T = Z( NN-7 ) + ( S+Z( NN-5 ) )
163 Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
164 Z( NN-7 ) = T
165 END IF
166 Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
167 Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
168 N0 = N0 - 2
169 GO TO 10
170 *
171 50 CONTINUE
172 *
173 * Reverse the qd-array, if warranted.
174 *
175 IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
176 IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
177 IPN4 = 4*( I0+N0 )
178 DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
179 TEMP = Z( J4-3 )
180 Z( J4-3 ) = Z( IPN4-J4-3 )
181 Z( IPN4-J4-3 ) = TEMP
182 TEMP = Z( J4-2 )
183 Z( J4-2 ) = Z( IPN4-J4-2 )
184 Z( IPN4-J4-2 ) = TEMP
185 TEMP = Z( J4-1 )
186 Z( J4-1 ) = Z( IPN4-J4-5 )
187 Z( IPN4-J4-5 ) = TEMP
188 TEMP = Z( J4 )
189 Z( J4 ) = Z( IPN4-J4-4 )
190 Z( IPN4-J4-4 ) = TEMP
191 60 CONTINUE
192 IF( N0-I0.LE.4 ) THEN
193 Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
194 Z( 4*N0-PP ) = Z( 4*I0-PP )
195 END IF
196 DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
197 Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
198 $ Z( 4*I0+PP+3 ) )
199 Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
200 $ Z( 4*I0-PP+4 ) )
201 QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
202 DMIN = -ZERO
203 END IF
204 END IF
205 *
206 IF( DMIN.LT.ZERO .OR. SAFMIN*QMAX.LT.MIN( Z( 4*N0+PP-1 ),
207 $ Z( 4*N0+PP-9 ), DMIN2+Z( 4*N0-PP ) ) ) THEN
208 *
209 * Choose a shift.
210 *
211 CALL SLAZQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
212 $ DN2, TAU, TTYPE, G )
213 *
214 * Call dqds until DMIN > 0.
215 *
216 80 CONTINUE
217 *
218 CALL SLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
219 $ DN1, DN2, IEEE )
220 *
221 NDIV = NDIV + ( N0-I0+2 )
222 ITER = ITER + 1
223 *
224 * Check status.
225 *
226 IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN
227 *
228 * Success.
229 *
230 GO TO 100
231 *
232 ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND.
233 $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
234 $ ABS( DN ).LT.TOL*SIGMA ) THEN
235 *
236 * Convergence hidden by negative DN.
237 *
238 Z( 4*( N0-1 )-PP+2 ) = ZERO
239 DMIN = ZERO
240 GO TO 100
241 ELSE IF( DMIN.LT.ZERO ) THEN
242 *
243 * TAU too big. Select new TAU and try again.
244 *
245 NFAIL = NFAIL + 1
246 IF( TTYPE.LT.-22 ) THEN
247 *
248 * Failed twice. Play it safe.
249 *
250 TAU = ZERO
251 ELSE IF( DMIN1.GT.ZERO ) THEN
252 *
253 * Late failure. Gives excellent shift.
254 *
255 TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
256 TTYPE = TTYPE - 11
257 ELSE
258 *
259 * Early failure. Divide by 4.
260 *
261 TAU = QURTR*TAU
262 TTYPE = TTYPE - 12
263 END IF
264 GO TO 80
265 ELSE IF( DMIN.NE.DMIN ) THEN
266 *
267 * NaN.
268 *
269 TAU = ZERO
270 GO TO 80
271 ELSE
272 *
273 * Possible underflow. Play it safe.
274 *
275 GO TO 90
276 END IF
277 END IF
278 *
279 * Risk of underflow.
280 *
281 90 CONTINUE
282 CALL SLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
283 NDIV = NDIV + ( N0-I0+2 )
284 ITER = ITER + 1
285 TAU = ZERO
286 *
287 100 CONTINUE
288 IF( TAU.LT.SIGMA ) THEN
289 DESIG = DESIG + TAU
290 T = SIGMA + DESIG
291 DESIG = DESIG - ( T-SIGMA )
292 ELSE
293 T = SIGMA + TAU
294 DESIG = SIGMA - ( T-TAU ) + DESIG
295 END IF
296 SIGMA = T
297 *
298 RETURN
299 *
300 * End of SLAZQ3
301 *
302 END