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