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