1: *> \brief \b DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLASQ2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASQ2( N, Z, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, N
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION Z( * )
28: * ..
29: *
30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> DLASQ2 computes all the eigenvalues of the symmetric positive
37: *> definite tridiagonal matrix associated with the qd array Z to high
38: *> relative accuracy are computed to high relative accuracy, in the
39: *> absence of denormalization, underflow and overflow.
40: *>
41: *> To see the relation of Z to the tridiagonal matrix, let L be a
42: *> unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
43: *> let U be an upper bidiagonal matrix with 1's above and diagonal
44: *> Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
45: *> symmetric tridiagonal to which it is similar.
46: *>
47: *> Note : DLASQ2 defines a logical variable, IEEE, which is true
48: *> on machines which follow ieee-754 floating-point standard in their
49: *> handling of infinities and NaNs, and false otherwise. This variable
50: *> is passed to DLASQ3.
51: *> \endverbatim
52: *
53: * Arguments:
54: * ==========
55: *
56: *> \param[in] N
57: *> \verbatim
58: *> N is INTEGER
59: *> The number of rows and columns in the matrix. N >= 0.
60: *> \endverbatim
61: *>
62: *> \param[in,out] Z
63: *> \verbatim
64: *> Z is DOUBLE PRECISION array, dimension ( 4*N )
65: *> On entry Z holds the qd array. On exit, entries 1 to N hold
66: *> the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
67: *> trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
68: *> N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
69: *> holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
70: *> shifts that failed.
71: *> \endverbatim
72: *>
73: *> \param[out] INFO
74: *> \verbatim
75: *> INFO is INTEGER
76: *> = 0: successful exit
77: *> < 0: if the i-th argument is a scalar and had an illegal
78: *> value, then INFO = -i, if the i-th argument is an
79: *> array and the j-entry had an illegal value, then
80: *> INFO = -(i*100+j)
81: *> > 0: the algorithm failed
82: *> = 1, a split was marked by a positive value in E
83: *> = 2, current block of Z not diagonalized after 100*N
84: *> iterations (in inner while loop). On exit Z holds
85: *> a qd array with the same eigenvalues as the given Z.
86: *> = 3, termination criterion of outer while loop not met
87: *> (program created more than N unreduced blocks)
88: *> \endverbatim
89: *
90: * Authors:
91: * ========
92: *
93: *> \author Univ. of Tennessee
94: *> \author Univ. of California Berkeley
95: *> \author Univ. of Colorado Denver
96: *> \author NAG Ltd.
97: *
98: *> \ingroup auxOTHERcomputational
99: *
100: *> \par Further Details:
101: * =====================
102: *>
103: *> \verbatim
104: *>
105: *> Local Variables: I0:N0 defines a current unreduced segment of Z.
106: *> The shifts are accumulated in SIGMA. Iteration count is in ITER.
107: *> Ping-pong is controlled by PP (alternates between 0 and 1).
108: *> \endverbatim
109: *>
110: * =====================================================================
111: SUBROUTINE DLASQ2( N, Z, INFO )
112: *
113: * -- LAPACK computational routine --
114: * -- LAPACK is a software package provided by Univ. of Tennessee, --
115: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
116: *
117: * .. Scalar Arguments ..
118: INTEGER INFO, N
119: * ..
120: * .. Array Arguments ..
121: DOUBLE PRECISION Z( * )
122: * ..
123: *
124: * =====================================================================
125: *
126: * .. Parameters ..
127: DOUBLE PRECISION CBIAS
128: PARAMETER ( CBIAS = 1.50D0 )
129: DOUBLE PRECISION ZERO, HALF, ONE, TWO, FOUR, HUNDRD
130: PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
131: $ TWO = 2.0D0, FOUR = 4.0D0, HUNDRD = 100.0D0 )
132: * ..
133: * .. Local Scalars ..
134: LOGICAL IEEE
135: INTEGER I0, I1, I4, IINFO, IPN4, ITER, IWHILA, IWHILB,
136: $ K, KMIN, N0, N1, NBIG, NDIV, NFAIL, PP, SPLT,
137: $ TTYPE
138: DOUBLE PRECISION D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN,
139: $ DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX,
140: $ QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL,
141: $ TOL2, TRACE, ZMAX, TEMPE, TEMPQ
142: * ..
143: * .. External Subroutines ..
144: EXTERNAL DLASQ3, DLASRT, XERBLA
145: * ..
146: * .. External Functions ..
147: INTEGER ILAENV
148: DOUBLE PRECISION DLAMCH
149: EXTERNAL DLAMCH, ILAENV
150: * ..
151: * .. Intrinsic Functions ..
152: INTRINSIC ABS, DBLE, MAX, MIN, SQRT
153: * ..
154: * .. Executable Statements ..
155: *
156: * Test the input arguments.
157: * (in case DLASQ2 is not called by DLASQ1)
158: *
159: INFO = 0
160: EPS = DLAMCH( 'Precision' )
161: SAFMIN = DLAMCH( 'Safe minimum' )
162: TOL = EPS*HUNDRD
163: TOL2 = TOL**2
164: *
165: IF( N.LT.0 ) THEN
166: INFO = -1
167: CALL XERBLA( 'DLASQ2', 1 )
168: RETURN
169: ELSE IF( N.EQ.0 ) THEN
170: RETURN
171: ELSE IF( N.EQ.1 ) THEN
172: *
173: * 1-by-1 case.
174: *
175: IF( Z( 1 ).LT.ZERO ) THEN
176: INFO = -201
177: CALL XERBLA( 'DLASQ2', 2 )
178: END IF
179: RETURN
180: ELSE IF( N.EQ.2 ) THEN
181: *
182: * 2-by-2 case.
183: *
184: IF( Z( 1 ).LT.ZERO ) THEN
185: INFO = -201
186: CALL XERBLA( 'DLASQ2', 2 )
187: RETURN
188: ELSE IF( Z( 2 ).LT.ZERO ) THEN
189: INFO = -202
190: CALL XERBLA( 'DLASQ2', 2 )
191: RETURN
192: ELSE IF( Z( 3 ).LT.ZERO ) THEN
193: INFO = -203
194: CALL XERBLA( 'DLASQ2', 2 )
195: RETURN
196: ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN
197: D = Z( 3 )
198: Z( 3 ) = Z( 1 )
199: Z( 1 ) = D
200: END IF
201: Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )
202: IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN
203: T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) )
204: S = Z( 3 )*( Z( 2 ) / T )
205: IF( S.LE.T ) THEN
206: S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )
207: ELSE
208: S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
209: END IF
210: T = Z( 1 ) + ( S+Z( 2 ) )
211: Z( 3 ) = Z( 3 )*( Z( 1 ) / T )
212: Z( 1 ) = T
213: END IF
214: Z( 2 ) = Z( 3 )
215: Z( 6 ) = Z( 2 ) + Z( 1 )
216: RETURN
217: END IF
218: *
219: * Check for negative data and compute sums of q's and e's.
220: *
221: Z( 2*N ) = ZERO
222: EMIN = Z( 2 )
223: QMAX = ZERO
224: ZMAX = ZERO
225: D = ZERO
226: E = ZERO
227: *
228: DO 10 K = 1, 2*( N-1 ), 2
229: IF( Z( K ).LT.ZERO ) THEN
230: INFO = -( 200+K )
231: CALL XERBLA( 'DLASQ2', 2 )
232: RETURN
233: ELSE IF( Z( K+1 ).LT.ZERO ) THEN
234: INFO = -( 200+K+1 )
235: CALL XERBLA( 'DLASQ2', 2 )
236: RETURN
237: END IF
238: D = D + Z( K )
239: E = E + Z( K+1 )
240: QMAX = MAX( QMAX, Z( K ) )
241: EMIN = MIN( EMIN, Z( K+1 ) )
242: ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) )
243: 10 CONTINUE
244: IF( Z( 2*N-1 ).LT.ZERO ) THEN
245: INFO = -( 200+2*N-1 )
246: CALL XERBLA( 'DLASQ2', 2 )
247: RETURN
248: END IF
249: D = D + Z( 2*N-1 )
250: QMAX = MAX( QMAX, Z( 2*N-1 ) )
251: ZMAX = MAX( QMAX, ZMAX )
252: *
253: * Check for diagonality.
254: *
255: IF( E.EQ.ZERO ) THEN
256: DO 20 K = 2, N
257: Z( K ) = Z( 2*K-1 )
258: 20 CONTINUE
259: CALL DLASRT( 'D', N, Z, IINFO )
260: Z( 2*N-1 ) = D
261: RETURN
262: END IF
263: *
264: TRACE = D + E
265: *
266: * Check for zero data.
267: *
268: IF( TRACE.EQ.ZERO ) THEN
269: Z( 2*N-1 ) = ZERO
270: RETURN
271: END IF
272: *
273: * Check whether the machine is IEEE conformable.
274: *
275: IEEE = ( ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 )
276: *
277: * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
278: *
279: DO 30 K = 2*N, 2, -2
280: Z( 2*K ) = ZERO
281: Z( 2*K-1 ) = Z( K )
282: Z( 2*K-2 ) = ZERO
283: Z( 2*K-3 ) = Z( K-1 )
284: 30 CONTINUE
285: *
286: I0 = 1
287: N0 = N
288: *
289: * Reverse the qd-array, if warranted.
290: *
291: IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN
292: IPN4 = 4*( I0+N0 )
293: DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4
294: TEMP = Z( I4-3 )
295: Z( I4-3 ) = Z( IPN4-I4-3 )
296: Z( IPN4-I4-3 ) = TEMP
297: TEMP = Z( I4-1 )
298: Z( I4-1 ) = Z( IPN4-I4-5 )
299: Z( IPN4-I4-5 ) = TEMP
300: 40 CONTINUE
301: END IF
302: *
303: * Initial split checking via dqd and Li's test.
304: *
305: PP = 0
306: *
307: DO 80 K = 1, 2
308: *
309: D = Z( 4*N0+PP-3 )
310: DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4
311: IF( Z( I4-1 ).LE.TOL2*D ) THEN
312: Z( I4-1 ) = -ZERO
313: D = Z( I4-3 )
314: ELSE
315: D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) )
316: END IF
317: 50 CONTINUE
318: *
319: * dqd maps Z to ZZ plus Li's test.
320: *
321: EMIN = Z( 4*I0+PP+1 )
322: D = Z( 4*I0+PP-3 )
323: DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4
324: Z( I4-2*PP-2 ) = D + Z( I4-1 )
325: IF( Z( I4-1 ).LE.TOL2*D ) THEN
326: Z( I4-1 ) = -ZERO
327: Z( I4-2*PP-2 ) = D
328: Z( I4-2*PP ) = ZERO
329: D = Z( I4+1 )
330: ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND.
331: $ SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN
332: TEMP = Z( I4+1 ) / Z( I4-2*PP-2 )
333: Z( I4-2*PP ) = Z( I4-1 )*TEMP
334: D = D*TEMP
335: ELSE
336: Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) )
337: D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) )
338: END IF
339: EMIN = MIN( EMIN, Z( I4-2*PP ) )
340: 60 CONTINUE
341: Z( 4*N0-PP-2 ) = D
342: *
343: * Now find qmax.
344: *
345: QMAX = Z( 4*I0-PP-2 )
346: DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4
347: QMAX = MAX( QMAX, Z( I4 ) )
348: 70 CONTINUE
349: *
350: * Prepare for the next iteration on K.
351: *
352: PP = 1 - PP
353: 80 CONTINUE
354: *
355: * Initialise variables to pass to DLASQ3.
356: *
357: TTYPE = 0
358: DMIN1 = ZERO
359: DMIN2 = ZERO
360: DN = ZERO
361: DN1 = ZERO
362: DN2 = ZERO
363: G = ZERO
364: TAU = ZERO
365: *
366: ITER = 2
367: NFAIL = 0
368: NDIV = 2*( N0-I0 )
369: *
370: DO 160 IWHILA = 1, N + 1
371: IF( N0.LT.1 )
372: $ GO TO 170
373: *
374: * While array unfinished do
375: *
376: * E(N0) holds the value of SIGMA when submatrix in I0:N0
377: * splits from the rest of the array, but is negated.
378: *
379: DESIG = ZERO
380: IF( N0.EQ.N ) THEN
381: SIGMA = ZERO
382: ELSE
383: SIGMA = -Z( 4*N0-1 )
384: END IF
385: IF( SIGMA.LT.ZERO ) THEN
386: INFO = 1
387: RETURN
388: END IF
389: *
390: * Find last unreduced submatrix's top index I0, find QMAX and
391: * EMIN. Find Gershgorin-type bound if Q's much greater than E's.
392: *
393: EMAX = ZERO
394: IF( N0.GT.I0 ) THEN
395: EMIN = ABS( Z( 4*N0-5 ) )
396: ELSE
397: EMIN = ZERO
398: END IF
399: QMIN = Z( 4*N0-3 )
400: QMAX = QMIN
401: DO 90 I4 = 4*N0, 8, -4
402: IF( Z( I4-5 ).LE.ZERO )
403: $ GO TO 100
404: IF( QMIN.GE.FOUR*EMAX ) THEN
405: QMIN = MIN( QMIN, Z( I4-3 ) )
406: EMAX = MAX( EMAX, Z( I4-5 ) )
407: END IF
408: QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) )
409: EMIN = MIN( EMIN, Z( I4-5 ) )
410: 90 CONTINUE
411: I4 = 4
412: *
413: 100 CONTINUE
414: I0 = I4 / 4
415: PP = 0
416: *
417: IF( N0-I0.GT.1 ) THEN
418: DEE = Z( 4*I0-3 )
419: DEEMIN = DEE
420: KMIN = I0
421: DO 110 I4 = 4*I0+1, 4*N0-3, 4
422: DEE = Z( I4 )*( DEE /( DEE+Z( I4-2 ) ) )
423: IF( DEE.LE.DEEMIN ) THEN
424: DEEMIN = DEE
425: KMIN = ( I4+3 )/4
426: END IF
427: 110 CONTINUE
428: IF( (KMIN-I0)*2.LT.N0-KMIN .AND.
429: $ DEEMIN.LE.HALF*Z(4*N0-3) ) THEN
430: IPN4 = 4*( I0+N0 )
431: PP = 2
432: DO 120 I4 = 4*I0, 2*( I0+N0-1 ), 4
433: TEMP = Z( I4-3 )
434: Z( I4-3 ) = Z( IPN4-I4-3 )
435: Z( IPN4-I4-3 ) = TEMP
436: TEMP = Z( I4-2 )
437: Z( I4-2 ) = Z( IPN4-I4-2 )
438: Z( IPN4-I4-2 ) = TEMP
439: TEMP = Z( I4-1 )
440: Z( I4-1 ) = Z( IPN4-I4-5 )
441: Z( IPN4-I4-5 ) = TEMP
442: TEMP = Z( I4 )
443: Z( I4 ) = Z( IPN4-I4-4 )
444: Z( IPN4-I4-4 ) = TEMP
445: 120 CONTINUE
446: END IF
447: END IF
448: *
449: * Put -(initial shift) into DMIN.
450: *
451: DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) )
452: *
453: * Now I0:N0 is unreduced.
454: * PP = 0 for ping, PP = 1 for pong.
455: * PP = 2 indicates that flipping was applied to the Z array and
456: * and that the tests for deflation upon entry in DLASQ3
457: * should not be performed.
458: *
459: NBIG = 100*( N0-I0+1 )
460: DO 140 IWHILB = 1, NBIG
461: IF( I0.GT.N0 )
462: $ GO TO 150
463: *
464: * While submatrix unfinished take a good dqds step.
465: *
466: CALL DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
467: $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
468: $ DN2, G, TAU )
469: *
470: PP = 1 - PP
471: *
472: * When EMIN is very small check for splits.
473: *
474: IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN
475: IF( Z( 4*N0 ).LE.TOL2*QMAX .OR.
476: $ Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN
477: SPLT = I0 - 1
478: QMAX = Z( 4*I0-3 )
479: EMIN = Z( 4*I0-1 )
480: OLDEMN = Z( 4*I0 )
481: DO 130 I4 = 4*I0, 4*( N0-3 ), 4
482: IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR.
483: $ Z( I4-1 ).LE.TOL2*SIGMA ) THEN
484: Z( I4-1 ) = -SIGMA
485: SPLT = I4 / 4
486: QMAX = ZERO
487: EMIN = Z( I4+3 )
488: OLDEMN = Z( I4+4 )
489: ELSE
490: QMAX = MAX( QMAX, Z( I4+1 ) )
491: EMIN = MIN( EMIN, Z( I4-1 ) )
492: OLDEMN = MIN( OLDEMN, Z( I4 ) )
493: END IF
494: 130 CONTINUE
495: Z( 4*N0-1 ) = EMIN
496: Z( 4*N0 ) = OLDEMN
497: I0 = SPLT + 1
498: END IF
499: END IF
500: *
501: 140 CONTINUE
502: *
503: INFO = 2
504: *
505: * Maximum number of iterations exceeded, restore the shift
506: * SIGMA and place the new d's and e's in a qd array.
507: * This might need to be done for several blocks
508: *
509: I1 = I0
510: N1 = N0
511: 145 CONTINUE
512: TEMPQ = Z( 4*I0-3 )
513: Z( 4*I0-3 ) = Z( 4*I0-3 ) + SIGMA
514: DO K = I0+1, N0
515: TEMPE = Z( 4*K-5 )
516: Z( 4*K-5 ) = Z( 4*K-5 ) * (TEMPQ / Z( 4*K-7 ))
517: TEMPQ = Z( 4*K-3 )
518: Z( 4*K-3 ) = Z( 4*K-3 ) + SIGMA + TEMPE - Z( 4*K-5 )
519: END DO
520: *
521: * Prepare to do this on the previous block if there is one
522: *
523: IF( I1.GT.1 ) THEN
524: N1 = I1-1
525: DO WHILE( ( I1.GE.2 ) .AND. ( Z(4*I1-5).GE.ZERO ) )
526: I1 = I1 - 1
527: END DO
528: SIGMA = -Z(4*N1-1)
529: GO TO 145
530: END IF
531:
532: DO K = 1, N
533: Z( 2*K-1 ) = Z( 4*K-3 )
534: *
535: * Only the block 1..N0 is unfinished. The rest of the e's
536: * must be essentially zero, although sometimes other data
537: * has been stored in them.
538: *
539: IF( K.LT.N0 ) THEN
540: Z( 2*K ) = Z( 4*K-1 )
541: ELSE
542: Z( 2*K ) = 0
543: END IF
544: END DO
545: RETURN
546: *
547: * end IWHILB
548: *
549: 150 CONTINUE
550: *
551: 160 CONTINUE
552: *
553: INFO = 3
554: RETURN
555: *
556: * end IWHILA
557: *
558: 170 CONTINUE
559: *
560: * Move q's to the front.
561: *
562: DO 180 K = 2, N
563: Z( K ) = Z( 4*K-3 )
564: 180 CONTINUE
565: *
566: * Sort and compute sum of eigenvalues.
567: *
568: CALL DLASRT( 'D', N, Z, IINFO )
569: *
570: E = ZERO
571: DO 190 K = N, 1, -1
572: E = E + Z( K )
573: 190 CONTINUE
574: *
575: * Store trace, sum(eigenvalues) and information on performance.
576: *
577: Z( 2*N+1 ) = TRACE
578: Z( 2*N+2 ) = E
579: Z( 2*N+3 ) = DBLE( ITER )
580: Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 )
581: Z( 2*N+5 ) = HUNDRD*NFAIL / DBLE( ITER )
582: RETURN
583: *
584: * End of DLASQ2
585: *
586: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack