|
3333
|
1 SUBROUTINE DLASQ2( N, Z, INFO ) |
|
2329
|
2 * |
|
3596
|
3 * -- LAPACK routine (version 3.0) -- |
|
2329
|
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
5 * Courant Institute, Argonne National Lab, and Rice University |
|
3596
|
6 * October 31, 1999 |
|
2329
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
3333
|
9 INTEGER INFO, N |
|
2329
|
10 * .. |
|
|
11 * .. Array Arguments .. |
|
3333
|
12 DOUBLE PRECISION Z( * ) |
|
2329
|
13 * .. |
|
|
14 * |
|
3333
|
15 * Purpose |
|
|
16 * ======= |
|
2329
|
17 * |
|
3596
|
18 * DLASQ2 computes all the eigenvalues of the symmetric positive |
|
3333
|
19 * definite tridiagonal matrix associated with the qd array Z to high |
|
|
20 * relative accuracy are computed to high relative accuracy, in the |
|
|
21 * absence of denormalization, underflow and overflow. |
|
2329
|
22 * |
|
3333
|
23 * To see the relation of Z to the tridiagonal matrix, let L be a |
|
|
24 * unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and |
|
|
25 * let U be an upper bidiagonal matrix with 1's above and diagonal |
|
|
26 * Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the |
|
|
27 * symmetric tridiagonal to which it is similar. |
|
2329
|
28 * |
|
3596
|
29 * Note : DLASQ2 defines a logical variable, IEEE, which is true |
|
|
30 * on machines which follow ieee-754 floating-point standard in their |
|
|
31 * handling of infinities and NaNs, and false otherwise. This variable |
|
|
32 * is passed to DLASQ3. |
|
2329
|
33 * |
|
3333
|
34 * Arguments |
|
|
35 * ========= |
|
|
36 * |
|
|
37 * N (input) INTEGER |
|
|
38 * The number of rows and columns in the matrix. N >= 0. |
|
2329
|
39 * |
|
3333
|
40 * Z (workspace) DOUBLE PRECISION array, dimension ( 4*N ) |
|
|
41 * On entry Z holds the qd array. On exit, entries 1 to N hold |
|
|
42 * the eigenvalues in decreasing order, Z( 2*N+1 ) holds the |
|
3596
|
43 * trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If |
|
|
44 * N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) |
|
|
45 * holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of |
|
|
46 * shifts that failed. |
|
2329
|
47 * |
|
3333
|
48 * INFO (output) INTEGER |
|
|
49 * = 0: successful exit |
|
|
50 * < 0: if the i-th argument is a scalar and had an illegal |
|
|
51 * value, then INFO = -i, if the i-th argument is an |
|
|
52 * array and the j-entry had an illegal value, then |
|
|
53 * INFO = -(i*100+j) |
|
|
54 * > 0: the algorithm failed |
|
|
55 * = 1, a split was marked by a positive value in E |
|
|
56 * = 2, current block of Z not diagonalized after 30*N |
|
|
57 * iterations (in inner while loop) |
|
3596
|
58 * = 3, termination criterion of outer while loop not met |
|
3333
|
59 * (program created more than N unreduced blocks) |
|
2329
|
60 * |
|
3333
|
61 * Further Details |
|
|
62 * =============== |
|
|
63 * Local Variables: I0:N0 defines a current unreduced segment of Z. |
|
|
64 * The shifts are accumulated in SIGMA. Iteration count is in ITER. |
|
|
65 * Ping-pong is controlled by PP (alternates between 0 and 1). |
|
2329
|
66 * |
|
|
67 * ===================================================================== |
|
|
68 * |
|
|
69 * .. Parameters .. |
|
3333
|
70 DOUBLE PRECISION CBIAS |
|
|
71 PARAMETER ( CBIAS = 1.50D0 ) |
|
3596
|
72 DOUBLE PRECISION ZERO, HALF, ONE, TWO, FOUR, HUNDRD |
|
3333
|
73 PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, |
|
3596
|
74 $ TWO = 2.0D0, FOUR = 4.0D0, HUNDRD = 100.0D0 ) |
|
2329
|
75 * .. |
|
|
76 * .. Local Scalars .. |
|
3596
|
77 LOGICAL IEEE |
|
|
78 INTEGER I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K, |
|
3333
|
79 $ N0, NBIG, NDIV, NFAIL, PP, SPLT |
|
3596
|
80 DOUBLE PRECISION D, DESIG, DMIN, E, EMAX, EMIN, EPS, OLDEMN, |
|
|
81 $ QMAX, QMIN, S, SAFMIN, SIGMA, T, TEMP, TOL, |
|
|
82 $ TOL2, TRACE, ZMAX |
|
2329
|
83 * .. |
|
|
84 * .. External Subroutines .. |
|
3596
|
85 EXTERNAL DLASQ3, DLASRT, XERBLA |
|
3333
|
86 * .. |
|
|
87 * .. External Functions .. |
|
3596
|
88 INTEGER ILAENV |
|
3333
|
89 DOUBLE PRECISION DLAMCH |
|
3596
|
90 EXTERNAL DLAMCH, ILAENV |
|
2329
|
91 * .. |
|
|
92 * .. Intrinsic Functions .. |
|
3596
|
93 INTRINSIC DBLE, MAX, MIN, SQRT |
|
2329
|
94 * .. |
|
|
95 * .. Executable Statements .. |
|
3596
|
96 * |
|
3333
|
97 * Test the input arguments. |
|
|
98 * (in case DLASQ2 is not called by DLASQ1) |
|
2329
|
99 * |
|
3333
|
100 INFO = 0 |
|
3596
|
101 EPS = DLAMCH( 'Precision' ) |
|
|
102 SAFMIN = DLAMCH( 'Safe minimum' ) |
|
|
103 TOL = EPS*HUNDRD |
|
|
104 TOL2 = TOL**2 |
|
2329
|
105 * |
|
3333
|
106 IF( N.LT.0 ) THEN |
|
|
107 INFO = -1 |
|
|
108 CALL XERBLA( 'DLASQ2', 1 ) |
|
|
109 RETURN |
|
|
110 ELSE IF( N.EQ.0 ) THEN |
|
|
111 RETURN |
|
|
112 ELSE IF( N.EQ.1 ) THEN |
|
2329
|
113 * |
|
3333
|
114 * 1-by-1 case. |
|
|
115 * |
|
|
116 IF( Z( 1 ).LT.ZERO ) THEN |
|
|
117 INFO = -201 |
|
|
118 CALL XERBLA( 'DLASQ2', 2 ) |
|
|
119 END IF |
|
|
120 RETURN |
|
|
121 ELSE IF( N.EQ.2 ) THEN |
|
|
122 * |
|
|
123 * 2-by-2 case. |
|
2329
|
124 * |
|
3333
|
125 IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN |
|
|
126 INFO = -2 |
|
|
127 CALL XERBLA( 'DLASQ2', 2 ) |
|
|
128 RETURN |
|
|
129 ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN |
|
|
130 D = Z( 3 ) |
|
|
131 Z( 3 ) = Z( 1 ) |
|
|
132 Z( 1 ) = D |
|
|
133 END IF |
|
|
134 Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 ) |
|
3596
|
135 IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN |
|
|
136 T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) |
|
3333
|
137 S = Z( 3 )*( Z( 2 ) / T ) |
|
|
138 IF( S.LE.T ) THEN |
|
|
139 S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) ) |
|
|
140 ELSE |
|
|
141 S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) ) |
|
|
142 END IF |
|
|
143 T = Z( 1 ) + ( S+Z( 2 ) ) |
|
|
144 Z( 3 ) = Z( 3 )*( Z( 1 ) / T ) |
|
|
145 Z( 1 ) = T |
|
|
146 END IF |
|
|
147 Z( 2 ) = Z( 3 ) |
|
|
148 Z( 6 ) = Z( 2 ) + Z( 1 ) |
|
|
149 RETURN |
|
2329
|
150 END IF |
|
|
151 * |
|
3333
|
152 * Check for negative data and compute sums of q's and e's. |
|
|
153 * |
|
|
154 Z( 2*N ) = ZERO |
|
|
155 EMIN = Z( 2 ) |
|
|
156 QMAX = ZERO |
|
3596
|
157 ZMAX = ZERO |
|
3333
|
158 D = ZERO |
|
|
159 E = ZERO |
|
|
160 * |
|
3596
|
161 DO 10 K = 1, 2*( N-1 ), 2 |
|
3333
|
162 IF( Z( K ).LT.ZERO ) THEN |
|
|
163 INFO = -( 200+K ) |
|
|
164 CALL XERBLA( 'DLASQ2', 2 ) |
|
|
165 RETURN |
|
3596
|
166 ELSE IF( Z( K+1 ).LT.ZERO ) THEN |
|
|
167 INFO = -( 200+K+1 ) |
|
3333
|
168 CALL XERBLA( 'DLASQ2', 2 ) |
|
|
169 RETURN |
|
|
170 END IF |
|
|
171 D = D + Z( K ) |
|
3596
|
172 E = E + Z( K+1 ) |
|
3333
|
173 QMAX = MAX( QMAX, Z( K ) ) |
|
3596
|
174 EMIN = MIN( EMIN, Z( K+1 ) ) |
|
|
175 ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) ) |
|
3333
|
176 10 CONTINUE |
|
3596
|
177 IF( Z( 2*N-1 ).LT.ZERO ) THEN |
|
|
178 INFO = -( 200+2*N-1 ) |
|
|
179 CALL XERBLA( 'DLASQ2', 2 ) |
|
|
180 RETURN |
|
|
181 END IF |
|
|
182 D = D + Z( 2*N-1 ) |
|
|
183 QMAX = MAX( QMAX, Z( 2*N-1 ) ) |
|
|
184 ZMAX = MAX( QMAX, ZMAX ) |
|
3333
|
185 * |
|
|
186 * Check for diagonality. |
|
2329
|
187 * |
|
3333
|
188 IF( E.EQ.ZERO ) THEN |
|
3596
|
189 DO 20 K = 2, N |
|
|
190 Z( K ) = Z( 2*K-1 ) |
|
|
191 20 CONTINUE |
|
3333
|
192 CALL DLASRT( 'D', N, Z, IINFO ) |
|
|
193 Z( 2*N-1 ) = D |
|
|
194 RETURN |
|
|
195 END IF |
|
|
196 * |
|
|
197 TRACE = D + E |
|
|
198 * |
|
|
199 * Check for zero data. |
|
|
200 * |
|
|
201 IF( TRACE.EQ.ZERO ) THEN |
|
|
202 Z( 2*N-1 ) = ZERO |
|
|
203 RETURN |
|
|
204 END IF |
|
3596
|
205 * |
|
|
206 * Check whether the machine is IEEE conformable. |
|
|
207 * |
|
|
208 IEEE = ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND. |
|
|
209 $ ILAENV( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 |
|
|
210 * |
|
3333
|
211 * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...). |
|
|
212 * |
|
|
213 DO 30 K = 2*N, 2, -2 |
|
3596
|
214 Z( 2*K ) = ZERO |
|
|
215 Z( 2*K-1 ) = Z( K ) |
|
|
216 Z( 2*K-2 ) = ZERO |
|
|
217 Z( 2*K-3 ) = Z( K-1 ) |
|
3333
|
218 30 CONTINUE |
|
|
219 * |
|
3596
|
220 I0 = 1 |
|
|
221 N0 = N |
|
|
222 * |
|
3333
|
223 * Reverse the qd-array, if warranted. |
|
|
224 * |
|
|
225 IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN |
|
|
226 IPN4 = 4*( I0+N0 ) |
|
|
227 DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4 |
|
|
228 TEMP = Z( I4-3 ) |
|
|
229 Z( I4-3 ) = Z( IPN4-I4-3 ) |
|
|
230 Z( IPN4-I4-3 ) = TEMP |
|
|
231 TEMP = Z( I4-1 ) |
|
|
232 Z( I4-1 ) = Z( IPN4-I4-5 ) |
|
|
233 Z( IPN4-I4-5 ) = TEMP |
|
|
234 40 CONTINUE |
|
2329
|
235 END IF |
|
|
236 * |
|
3333
|
237 * Initial split checking via dqd and Li's test. |
|
|
238 * |
|
|
239 PP = 0 |
|
2329
|
240 * |
|
3333
|
241 DO 80 K = 1, 2 |
|
2329
|
242 * |
|
3596
|
243 D = Z( 4*N0+PP-3 ) |
|
|
244 DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4 |
|
|
245 IF( Z( I4-1 ).LE.TOL2*D ) THEN |
|
|
246 Z( I4-1 ) = -ZERO |
|
|
247 D = Z( I4-3 ) |
|
|
248 ELSE |
|
|
249 D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) ) |
|
|
250 END IF |
|
|
251 50 CONTINUE |
|
2329
|
252 * |
|
3596
|
253 * dqd maps Z to ZZ plus Li's test. |
|
2329
|
254 * |
|
3596
|
255 EMIN = Z( 4*I0+PP+1 ) |
|
|
256 D = Z( 4*I0+PP-3 ) |
|
|
257 DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4 |
|
|
258 Z( I4-2*PP-2 ) = D + Z( I4-1 ) |
|
|
259 IF( Z( I4-1 ).LE.TOL2*D ) THEN |
|
|
260 Z( I4-1 ) = -ZERO |
|
|
261 Z( I4-2*PP-2 ) = D |
|
|
262 Z( I4-2*PP ) = ZERO |
|
|
263 D = Z( I4+1 ) |
|
|
264 ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND. |
|
|
265 $ SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN |
|
|
266 TEMP = Z( I4+1 ) / Z( I4-2*PP-2 ) |
|
|
267 Z( I4-2*PP ) = Z( I4-1 )*TEMP |
|
|
268 D = D*TEMP |
|
|
269 ELSE |
|
|
270 Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) ) |
|
|
271 D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) ) |
|
|
272 END IF |
|
|
273 EMIN = MIN( EMIN, Z( I4-2*PP ) ) |
|
|
274 60 CONTINUE |
|
|
275 Z( 4*N0-PP-2 ) = D |
|
3333
|
276 * |
|
|
277 * Now find qmax. |
|
|
278 * |
|
|
279 QMAX = Z( 4*I0-PP-2 ) |
|
|
280 DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4 |
|
|
281 QMAX = MAX( QMAX, Z( I4 ) ) |
|
|
282 70 CONTINUE |
|
|
283 * |
|
|
284 * Prepare for the next iteration on K. |
|
|
285 * |
|
|
286 PP = 1 - PP |
|
|
287 80 CONTINUE |
|
2329
|
288 * |
|
3333
|
289 ITER = 2 |
|
|
290 NFAIL = 0 |
|
|
291 NDIV = 2*( N0-I0 ) |
|
|
292 * |
|
|
293 DO 140 IWHILA = 1, N + 1 |
|
3596
|
294 IF( N0.LT.1 ) |
|
3333
|
295 $ GO TO 150 |
|
2329
|
296 * |
|
3596
|
297 * While array unfinished do |
|
3333
|
298 * |
|
|
299 * E(N0) holds the value of SIGMA when submatrix in I0:N0 |
|
|
300 * splits from the rest of the array, but is negated. |
|
3596
|
301 * |
|
3333
|
302 DESIG = ZERO |
|
|
303 IF( N0.EQ.N ) THEN |
|
|
304 SIGMA = ZERO |
|
2329
|
305 ELSE |
|
3333
|
306 SIGMA = -Z( 4*N0-1 ) |
|
2329
|
307 END IF |
|
3333
|
308 IF( SIGMA.LT.ZERO ) THEN |
|
|
309 INFO = 1 |
|
|
310 RETURN |
|
|
311 END IF |
|
2329
|
312 * |
|
3333
|
313 * Find last unreduced submatrix's top index I0, find QMAX and |
|
|
314 * EMIN. Find Gershgorin-type bound if Q's much greater than E's. |
|
2329
|
315 * |
|
3596
|
316 EMAX = ZERO |
|
|
317 IF( N0.GT.I0 ) THEN |
|
|
318 EMIN = ABS( Z( 4*N0-5 ) ) |
|
|
319 ELSE |
|
|
320 EMIN = ZERO |
|
|
321 END IF |
|
3333
|
322 QMIN = Z( 4*N0-3 ) |
|
|
323 QMAX = QMIN |
|
|
324 DO 90 I4 = 4*N0, 8, -4 |
|
|
325 IF( Z( I4-5 ).LE.ZERO ) |
|
|
326 $ GO TO 100 |
|
|
327 IF( QMIN.GE.FOUR*EMAX ) THEN |
|
|
328 QMIN = MIN( QMIN, Z( I4-3 ) ) |
|
|
329 EMAX = MAX( EMAX, Z( I4-5 ) ) |
|
2329
|
330 END IF |
|
3333
|
331 QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) ) |
|
|
332 EMIN = MIN( EMIN, Z( I4-5 ) ) |
|
|
333 90 CONTINUE |
|
3596
|
334 I4 = 4 |
|
3333
|
335 * |
|
|
336 100 CONTINUE |
|
|
337 I0 = I4 / 4 |
|
|
338 * |
|
|
339 * Store EMIN for passing to DLASQ3. |
|
|
340 * |
|
|
341 Z( 4*N0-1 ) = EMIN |
|
|
342 * |
|
|
343 * Put -(initial shift) into DMIN. |
|
|
344 * |
|
|
345 DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) ) |
|
|
346 * |
|
|
347 * Now I0:N0 is unreduced. PP = 0 for ping, PP = 1 for pong. |
|
|
348 * |
|
3596
|
349 PP = 0 |
|
3333
|
350 * |
|
|
351 NBIG = 30*( N0-I0+1 ) |
|
|
352 DO 120 IWHILB = 1, NBIG |
|
3596
|
353 IF( I0.GT.N0 ) |
|
3333
|
354 $ GO TO 130 |
|
|
355 * |
|
|
356 * While submatrix unfinished take a good dqds step. |
|
|
357 * |
|
|
358 CALL DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, |
|
3596
|
359 $ ITER, NDIV, IEEE ) |
|
3333
|
360 * |
|
|
361 PP = 1 - PP |
|
|
362 * |
|
|
363 * When EMIN is very small check for splits. |
|
|
364 * |
|
|
365 IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN |
|
3596
|
366 IF( Z( 4*N0 ).LE.TOL2*QMAX .OR. |
|
|
367 $ Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN |
|
3333
|
368 SPLT = I0 - 1 |
|
|
369 QMAX = Z( 4*I0-3 ) |
|
|
370 EMIN = Z( 4*I0-1 ) |
|
|
371 OLDEMN = Z( 4*I0 ) |
|
|
372 DO 110 I4 = 4*I0, 4*( N0-3 ), 4 |
|
3596
|
373 IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR. |
|
|
374 $ Z( I4-1 ).LE.TOL2*SIGMA ) THEN |
|
3333
|
375 Z( I4-1 ) = -SIGMA |
|
|
376 SPLT = I4 / 4 |
|
|
377 QMAX = ZERO |
|
|
378 EMIN = Z( I4+3 ) |
|
|
379 OLDEMN = Z( I4+4 ) |
|
|
380 ELSE |
|
|
381 QMAX = MAX( QMAX, Z( I4+1 ) ) |
|
|
382 EMIN = MIN( EMIN, Z( I4-1 ) ) |
|
|
383 OLDEMN = MIN( OLDEMN, Z( I4 ) ) |
|
|
384 END IF |
|
|
385 110 CONTINUE |
|
|
386 Z( 4*N0-1 ) = EMIN |
|
|
387 Z( 4*N0 ) = OLDEMN |
|
|
388 I0 = SPLT + 1 |
|
2329
|
389 END IF |
|
|
390 END IF |
|
3333
|
391 * |
|
|
392 120 CONTINUE |
|
|
393 * |
|
|
394 INFO = 2 |
|
2329
|
395 RETURN |
|
3333
|
396 * |
|
|
397 * end IWHILB |
|
|
398 * |
|
|
399 130 CONTINUE |
|
|
400 * |
|
|
401 140 CONTINUE |
|
|
402 * |
|
|
403 INFO = 3 |
|
|
404 RETURN |
|
|
405 * |
|
3596
|
406 * end IWHILA |
|
3333
|
407 * |
|
|
408 150 CONTINUE |
|
3596
|
409 * |
|
3333
|
410 * Move q's to the front. |
|
3596
|
411 * |
|
3333
|
412 DO 160 K = 2, N |
|
|
413 Z( K ) = Z( 4*K-3 ) |
|
|
414 160 CONTINUE |
|
3596
|
415 * |
|
3333
|
416 * Sort and compute sum of eigenvalues. |
|
|
417 * |
|
|
418 CALL DLASRT( 'D', N, Z, IINFO ) |
|
|
419 * |
|
|
420 E = ZERO |
|
|
421 DO 170 K = N, 1, -1 |
|
|
422 E = E + Z( K ) |
|
|
423 170 CONTINUE |
|
|
424 * |
|
|
425 * Store trace, sum(eigenvalues) and information on performance. |
|
|
426 * |
|
3596
|
427 Z( 2*N+1 ) = TRACE |
|
3333
|
428 Z( 2*N+2 ) = E |
|
|
429 Z( 2*N+3 ) = DBLE( ITER ) |
|
|
430 Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 ) |
|
3596
|
431 Z( 2*N+5 ) = HUNDRD*NFAIL / DBLE( ITER ) |
|
3333
|
432 RETURN |
|
2329
|
433 * |
|
|
434 * End of DLASQ2 |
|
|
435 * |
|
|
436 END |
