![[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|
Return to dlasq2.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 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