1: SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
2: $ DN1, DN2, TAU, TTYPE, G )
3: *
4: * -- LAPACK routine (version 3.3.1) --
5: *
6: * -- Contributed by Osni Marques of the Lawrence Berkeley National --
7: * -- Laboratory and Beresford Parlett of the Univ. of California at --
8: * -- Berkeley --
9: * -- November 2008 --
10: *
11: * -- LAPACK is a software package provided by Univ. of Tennessee, --
12: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
13: *
14: * .. Scalar Arguments ..
15: INTEGER I0, N0, N0IN, PP, TTYPE
16: DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
17: * ..
18: * .. Array Arguments ..
19: DOUBLE PRECISION Z( * )
20: * ..
21: *
22: * Purpose
23: * =======
24: *
25: * DLASQ4 computes an approximation TAU to the smallest eigenvalue
26: * using values of d from the previous transform.
27: *
28: * Arguments
29: * =========
30: *
31: * I0 (input) INTEGER
32: * First index.
33: *
34: * N0 (input) INTEGER
35: * Last index.
36: *
37: * Z (input) DOUBLE PRECISION array, dimension ( 4*N )
38: * Z holds the qd array.
39: *
40: * PP (input) INTEGER
41: * PP=0 for ping, PP=1 for pong.
42: *
43: * NOIN (input) INTEGER
44: * The value of N0 at start of EIGTEST.
45: *
46: * DMIN (input) DOUBLE PRECISION
47: * Minimum value of d.
48: *
49: * DMIN1 (input) DOUBLE PRECISION
50: * Minimum value of d, excluding D( N0 ).
51: *
52: * DMIN2 (input) DOUBLE PRECISION
53: * Minimum value of d, excluding D( N0 ) and D( N0-1 ).
54: *
55: * DN (input) DOUBLE PRECISION
56: * d(N)
57: *
58: * DN1 (input) DOUBLE PRECISION
59: * d(N-1)
60: *
61: * DN2 (input) DOUBLE PRECISION
62: * d(N-2)
63: *
64: * TAU (output) DOUBLE PRECISION
65: * This is the shift.
66: *
67: * TTYPE (output) INTEGER
68: * Shift type.
69: *
70: * G (input/output) REAL
71: * G is passed as an argument in order to save its value between
72: * calls to DLASQ4.
73: *
74: * Further Details
75: * ===============
76: * CNST1 = 9/16
77: *
78: * =====================================================================
79: *
80: * .. Parameters ..
81: DOUBLE PRECISION CNST1, CNST2, CNST3
82: PARAMETER ( CNST1 = 0.5630D0, CNST2 = 1.010D0,
83: $ CNST3 = 1.050D0 )
84: DOUBLE PRECISION QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
85: PARAMETER ( QURTR = 0.250D0, THIRD = 0.3330D0,
86: $ HALF = 0.50D0, ZERO = 0.0D0, ONE = 1.0D0,
87: $ TWO = 2.0D0, HUNDRD = 100.0D0 )
88: * ..
89: * .. Local Scalars ..
90: INTEGER I4, NN, NP
91: DOUBLE PRECISION A2, B1, B2, GAM, GAP1, GAP2, S
92: * ..
93: * .. Intrinsic Functions ..
94: INTRINSIC MAX, MIN, SQRT
95: * ..
96: * .. Executable Statements ..
97: *
98: * A negative DMIN forces the shift to take that absolute value
99: * TTYPE records the type of shift.
100: *
101: IF( DMIN.LE.ZERO ) THEN
102: TAU = -DMIN
103: TTYPE = -1
104: RETURN
105: END IF
106: *
107: NN = 4*N0 + PP
108: IF( N0IN.EQ.N0 ) THEN
109: *
110: * No eigenvalues deflated.
111: *
112: IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN
113: *
114: B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) )
115: B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) )
116: A2 = Z( NN-7 ) + Z( NN-5 )
117: *
118: * Cases 2 and 3.
119: *
120: IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN
121: GAP2 = DMIN2 - A2 - DMIN2*QURTR
122: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN
123: GAP1 = A2 - DN - ( B2 / GAP2 )*B2
124: ELSE
125: GAP1 = A2 - DN - ( B1+B2 )
126: END IF
127: IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN
128: S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN )
129: TTYPE = -2
130: ELSE
131: S = ZERO
132: IF( DN.GT.B1 )
133: $ S = DN - B1
134: IF( A2.GT.( B1+B2 ) )
135: $ S = MIN( S, A2-( B1+B2 ) )
136: S = MAX( S, THIRD*DMIN )
137: TTYPE = -3
138: END IF
139: ELSE
140: *
141: * Case 4.
142: *
143: TTYPE = -4
144: S = QURTR*DMIN
145: IF( DMIN.EQ.DN ) THEN
146: GAM = DN
147: A2 = ZERO
148: IF( Z( NN-5 ) .GT. Z( NN-7 ) )
149: $ RETURN
150: B2 = Z( NN-5 ) / Z( NN-7 )
151: NP = NN - 9
152: ELSE
153: NP = NN - 2*PP
154: B2 = Z( NP-2 )
155: GAM = DN1
156: IF( Z( NP-4 ) .GT. Z( NP-2 ) )
157: $ RETURN
158: A2 = Z( NP-4 ) / Z( NP-2 )
159: IF( Z( NN-9 ) .GT. Z( NN-11 ) )
160: $ RETURN
161: B2 = Z( NN-9 ) / Z( NN-11 )
162: NP = NN - 13
163: END IF
164: *
165: * Approximate contribution to norm squared from I < NN-1.
166: *
167: A2 = A2 + B2
168: DO 10 I4 = NP, 4*I0 - 1 + PP, -4
169: IF( B2.EQ.ZERO )
170: $ GO TO 20
171: B1 = B2
172: IF( Z( I4 ) .GT. Z( I4-2 ) )
173: $ RETURN
174: B2 = B2*( Z( I4 ) / Z( I4-2 ) )
175: A2 = A2 + B2
176: IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
177: $ GO TO 20
178: 10 CONTINUE
179: 20 CONTINUE
180: A2 = CNST3*A2
181: *
182: * Rayleigh quotient residual bound.
183: *
184: IF( A2.LT.CNST1 )
185: $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
186: END IF
187: ELSE IF( DMIN.EQ.DN2 ) THEN
188: *
189: * Case 5.
190: *
191: TTYPE = -5
192: S = QURTR*DMIN
193: *
194: * Compute contribution to norm squared from I > NN-2.
195: *
196: NP = NN - 2*PP
197: B1 = Z( NP-2 )
198: B2 = Z( NP-6 )
199: GAM = DN2
200: IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 )
201: $ RETURN
202: A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 )
203: *
204: * Approximate contribution to norm squared from I < NN-2.
205: *
206: IF( N0-I0.GT.2 ) THEN
207: B2 = Z( NN-13 ) / Z( NN-15 )
208: A2 = A2 + B2
209: DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4
210: IF( B2.EQ.ZERO )
211: $ GO TO 40
212: B1 = B2
213: IF( Z( I4 ) .GT. Z( I4-2 ) )
214: $ RETURN
215: B2 = B2*( Z( I4 ) / Z( I4-2 ) )
216: A2 = A2 + B2
217: IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
218: $ GO TO 40
219: 30 CONTINUE
220: 40 CONTINUE
221: A2 = CNST3*A2
222: END IF
223: *
224: IF( A2.LT.CNST1 )
225: $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
226: ELSE
227: *
228: * Case 6, no information to guide us.
229: *
230: IF( TTYPE.EQ.-6 ) THEN
231: G = G + THIRD*( ONE-G )
232: ELSE IF( TTYPE.EQ.-18 ) THEN
233: G = QURTR*THIRD
234: ELSE
235: G = QURTR
236: END IF
237: S = G*DMIN
238: TTYPE = -6
239: END IF
240: *
241: ELSE IF( N0IN.EQ.( N0+1 ) ) THEN
242: *
243: * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
244: *
245: IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN
246: *
247: * Cases 7 and 8.
248: *
249: TTYPE = -7
250: S = THIRD*DMIN1
251: IF( Z( NN-5 ).GT.Z( NN-7 ) )
252: $ RETURN
253: B1 = Z( NN-5 ) / Z( NN-7 )
254: B2 = B1
255: IF( B2.EQ.ZERO )
256: $ GO TO 60
257: DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
258: A2 = B1
259: IF( Z( I4 ).GT.Z( I4-2 ) )
260: $ RETURN
261: B1 = B1*( Z( I4 ) / Z( I4-2 ) )
262: B2 = B2 + B1
263: IF( HUNDRD*MAX( B1, A2 ).LT.B2 )
264: $ GO TO 60
265: 50 CONTINUE
266: 60 CONTINUE
267: B2 = SQRT( CNST3*B2 )
268: A2 = DMIN1 / ( ONE+B2**2 )
269: GAP2 = HALF*DMIN2 - A2
270: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
271: S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
272: ELSE
273: S = MAX( S, A2*( ONE-CNST2*B2 ) )
274: TTYPE = -8
275: END IF
276: ELSE
277: *
278: * Case 9.
279: *
280: S = QURTR*DMIN1
281: IF( DMIN1.EQ.DN1 )
282: $ S = HALF*DMIN1
283: TTYPE = -9
284: END IF
285: *
286: ELSE IF( N0IN.EQ.( N0+2 ) ) THEN
287: *
288: * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
289: *
290: * Cases 10 and 11.
291: *
292: IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN
293: TTYPE = -10
294: S = THIRD*DMIN2
295: IF( Z( NN-5 ).GT.Z( NN-7 ) )
296: $ RETURN
297: B1 = Z( NN-5 ) / Z( NN-7 )
298: B2 = B1
299: IF( B2.EQ.ZERO )
300: $ GO TO 80
301: DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
302: IF( Z( I4 ).GT.Z( I4-2 ) )
303: $ RETURN
304: B1 = B1*( Z( I4 ) / Z( I4-2 ) )
305: B2 = B2 + B1
306: IF( HUNDRD*B1.LT.B2 )
307: $ GO TO 80
308: 70 CONTINUE
309: 80 CONTINUE
310: B2 = SQRT( CNST3*B2 )
311: A2 = DMIN2 / ( ONE+B2**2 )
312: GAP2 = Z( NN-7 ) + Z( NN-9 ) -
313: $ SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2
314: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
315: S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
316: ELSE
317: S = MAX( S, A2*( ONE-CNST2*B2 ) )
318: END IF
319: ELSE
320: S = QURTR*DMIN2
321: TTYPE = -11
322: END IF
323: ELSE IF( N0IN.GT.( N0+2 ) ) THEN
324: *
325: * Case 12, more than two eigenvalues deflated. No information.
326: *
327: S = ZERO
328: TTYPE = -12
329: END IF
330: *
331: TAU = S
332: RETURN
333: *
334: * End of DLASQ4
335: *
336: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq4.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack