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