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