1: *> \brief \b DLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLASQ4 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq4.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq4.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq4.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
22: * DN1, DN2, TAU, TTYPE, G )
23: *
24: * .. Scalar Arguments ..
25: * INTEGER I0, N0, N0IN, PP, TTYPE
26: * DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION Z( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DLASQ4 computes an approximation TAU to the smallest eigenvalue
39: *> using values of d from the previous transform.
40: *> \endverbatim
41: *
42: * Arguments:
43: * ==========
44: *
45: *> \param[in] I0
46: *> \verbatim
47: *> I0 is INTEGER
48: *> First index.
49: *> \endverbatim
50: *>
51: *> \param[in] N0
52: *> \verbatim
53: *> N0 is INTEGER
54: *> Last index.
55: *> \endverbatim
56: *>
57: *> \param[in] Z
58: *> \verbatim
59: *> Z is DOUBLE PRECISION array, dimension ( 4*N0 )
60: *> Z holds the qd array.
61: *> \endverbatim
62: *>
63: *> \param[in] PP
64: *> \verbatim
65: *> PP is INTEGER
66: *> PP=0 for ping, PP=1 for pong.
67: *> \endverbatim
68: *>
69: *> \param[in] N0IN
70: *> \verbatim
71: *> N0IN is INTEGER
72: *> The value of N0 at start of EIGTEST.
73: *> \endverbatim
74: *>
75: *> \param[in] DMIN
76: *> \verbatim
77: *> DMIN is DOUBLE PRECISION
78: *> Minimum value of d.
79: *> \endverbatim
80: *>
81: *> \param[in] DMIN1
82: *> \verbatim
83: *> DMIN1 is DOUBLE PRECISION
84: *> Minimum value of d, excluding D( N0 ).
85: *> \endverbatim
86: *>
87: *> \param[in] DMIN2
88: *> \verbatim
89: *> DMIN2 is DOUBLE PRECISION
90: *> Minimum value of d, excluding D( N0 ) and D( N0-1 ).
91: *> \endverbatim
92: *>
93: *> \param[in] DN
94: *> \verbatim
95: *> DN is DOUBLE PRECISION
96: *> d(N)
97: *> \endverbatim
98: *>
99: *> \param[in] DN1
100: *> \verbatim
101: *> DN1 is DOUBLE PRECISION
102: *> d(N-1)
103: *> \endverbatim
104: *>
105: *> \param[in] DN2
106: *> \verbatim
107: *> DN2 is DOUBLE PRECISION
108: *> d(N-2)
109: *> \endverbatim
110: *>
111: *> \param[out] TAU
112: *> \verbatim
113: *> TAU is DOUBLE PRECISION
114: *> This is the shift.
115: *> \endverbatim
116: *>
117: *> \param[out] TTYPE
118: *> \verbatim
119: *> TTYPE is INTEGER
120: *> Shift type.
121: *> \endverbatim
122: *>
123: *> \param[in,out] G
124: *> \verbatim
125: *> G is DOUBLE PRECISION
126: *> G is passed as an argument in order to save its value between
127: *> calls to DLASQ4.
128: *> \endverbatim
129: *
130: * Authors:
131: * ========
132: *
133: *> \author Univ. of Tennessee
134: *> \author Univ. of California Berkeley
135: *> \author Univ. of Colorado Denver
136: *> \author NAG Ltd.
137: *
138: *> \ingroup auxOTHERcomputational
139: *
140: *> \par Further Details:
141: * =====================
142: *>
143: *> \verbatim
144: *>
145: *> CNST1 = 9/16
146: *> \endverbatim
147: *>
148: * =====================================================================
149: SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
150: $ DN1, DN2, TAU, TTYPE, G )
151: *
152: * -- LAPACK computational routine --
153: * -- LAPACK is a software package provided by Univ. of Tennessee, --
154: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155: *
156: * .. Scalar Arguments ..
157: INTEGER I0, N0, N0IN, PP, TTYPE
158: DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
159: * ..
160: * .. Array Arguments ..
161: DOUBLE PRECISION Z( * )
162: * ..
163: *
164: * =====================================================================
165: *
166: * .. Parameters ..
167: DOUBLE PRECISION CNST1, CNST2, CNST3
168: PARAMETER ( CNST1 = 0.5630D0, CNST2 = 1.010D0,
169: $ CNST3 = 1.050D0 )
170: DOUBLE PRECISION QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
171: PARAMETER ( QURTR = 0.250D0, THIRD = 0.3330D0,
172: $ HALF = 0.50D0, ZERO = 0.0D0, ONE = 1.0D0,
173: $ TWO = 2.0D0, HUNDRD = 100.0D0 )
174: * ..
175: * .. Local Scalars ..
176: INTEGER I4, NN, NP
177: DOUBLE PRECISION A2, B1, B2, GAM, GAP1, GAP2, S
178: * ..
179: * .. Intrinsic Functions ..
180: INTRINSIC MAX, MIN, SQRT
181: * ..
182: * .. Executable Statements ..
183: *
184: * A negative DMIN forces the shift to take that absolute value
185: * TTYPE records the type of shift.
186: *
187: IF( DMIN.LE.ZERO ) THEN
188: TAU = -DMIN
189: TTYPE = -1
190: RETURN
191: END IF
192: *
193: NN = 4*N0 + PP
194: IF( N0IN.EQ.N0 ) THEN
195: *
196: * No eigenvalues deflated.
197: *
198: IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN
199: *
200: B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) )
201: B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) )
202: A2 = Z( NN-7 ) + Z( NN-5 )
203: *
204: * Cases 2 and 3.
205: *
206: IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN
207: GAP2 = DMIN2 - A2 - DMIN2*QURTR
208: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN
209: GAP1 = A2 - DN - ( B2 / GAP2 )*B2
210: ELSE
211: GAP1 = A2 - DN - ( B1+B2 )
212: END IF
213: IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN
214: S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN )
215: TTYPE = -2
216: ELSE
217: S = ZERO
218: IF( DN.GT.B1 )
219: $ S = DN - B1
220: IF( A2.GT.( B1+B2 ) )
221: $ S = MIN( S, A2-( B1+B2 ) )
222: S = MAX( S, THIRD*DMIN )
223: TTYPE = -3
224: END IF
225: ELSE
226: *
227: * Case 4.
228: *
229: TTYPE = -4
230: S = QURTR*DMIN
231: IF( DMIN.EQ.DN ) THEN
232: GAM = DN
233: A2 = ZERO
234: IF( Z( NN-5 ) .GT. Z( NN-7 ) )
235: $ RETURN
236: B2 = Z( NN-5 ) / Z( NN-7 )
237: NP = NN - 9
238: ELSE
239: NP = NN - 2*PP
240: GAM = DN1
241: IF( Z( NP-4 ) .GT. Z( NP-2 ) )
242: $ RETURN
243: A2 = Z( NP-4 ) / Z( NP-2 )
244: IF( Z( NN-9 ) .GT. Z( NN-11 ) )
245: $ RETURN
246: B2 = Z( NN-9 ) / Z( NN-11 )
247: NP = NN - 13
248: END IF
249: *
250: * Approximate contribution to norm squared from I < NN-1.
251: *
252: A2 = A2 + B2
253: DO 10 I4 = NP, 4*I0 - 1 + PP, -4
254: IF( B2.EQ.ZERO )
255: $ GO TO 20
256: B1 = B2
257: IF( Z( I4 ) .GT. Z( I4-2 ) )
258: $ RETURN
259: B2 = B2*( Z( I4 ) / Z( I4-2 ) )
260: A2 = A2 + B2
261: IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
262: $ GO TO 20
263: 10 CONTINUE
264: 20 CONTINUE
265: A2 = CNST3*A2
266: *
267: * Rayleigh quotient residual bound.
268: *
269: IF( A2.LT.CNST1 )
270: $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
271: END IF
272: ELSE IF( DMIN.EQ.DN2 ) THEN
273: *
274: * Case 5.
275: *
276: TTYPE = -5
277: S = QURTR*DMIN
278: *
279: * Compute contribution to norm squared from I > NN-2.
280: *
281: NP = NN - 2*PP
282: B1 = Z( NP-2 )
283: B2 = Z( NP-6 )
284: GAM = DN2
285: IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 )
286: $ RETURN
287: A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 )
288: *
289: * Approximate contribution to norm squared from I < NN-2.
290: *
291: IF( N0-I0.GT.2 ) THEN
292: B2 = Z( NN-13 ) / Z( NN-15 )
293: A2 = A2 + B2
294: DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4
295: IF( B2.EQ.ZERO )
296: $ GO TO 40
297: B1 = B2
298: IF( Z( I4 ) .GT. Z( I4-2 ) )
299: $ RETURN
300: B2 = B2*( Z( I4 ) / Z( I4-2 ) )
301: A2 = A2 + B2
302: IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
303: $ GO TO 40
304: 30 CONTINUE
305: 40 CONTINUE
306: A2 = CNST3*A2
307: END IF
308: *
309: IF( A2.LT.CNST1 )
310: $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
311: ELSE
312: *
313: * Case 6, no information to guide us.
314: *
315: IF( TTYPE.EQ.-6 ) THEN
316: G = G + THIRD*( ONE-G )
317: ELSE IF( TTYPE.EQ.-18 ) THEN
318: G = QURTR*THIRD
319: ELSE
320: G = QURTR
321: END IF
322: S = G*DMIN
323: TTYPE = -6
324: END IF
325: *
326: ELSE IF( N0IN.EQ.( N0+1 ) ) THEN
327: *
328: * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
329: *
330: IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN
331: *
332: * Cases 7 and 8.
333: *
334: TTYPE = -7
335: S = THIRD*DMIN1
336: IF( Z( NN-5 ).GT.Z( NN-7 ) )
337: $ RETURN
338: B1 = Z( NN-5 ) / Z( NN-7 )
339: B2 = B1
340: IF( B2.EQ.ZERO )
341: $ GO TO 60
342: DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
343: A2 = B1
344: IF( Z( I4 ).GT.Z( I4-2 ) )
345: $ RETURN
346: B1 = B1*( Z( I4 ) / Z( I4-2 ) )
347: B2 = B2 + B1
348: IF( HUNDRD*MAX( B1, A2 ).LT.B2 )
349: $ GO TO 60
350: 50 CONTINUE
351: 60 CONTINUE
352: B2 = SQRT( CNST3*B2 )
353: A2 = DMIN1 / ( ONE+B2**2 )
354: GAP2 = HALF*DMIN2 - A2
355: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
356: S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
357: ELSE
358: S = MAX( S, A2*( ONE-CNST2*B2 ) )
359: TTYPE = -8
360: END IF
361: ELSE
362: *
363: * Case 9.
364: *
365: S = QURTR*DMIN1
366: IF( DMIN1.EQ.DN1 )
367: $ S = HALF*DMIN1
368: TTYPE = -9
369: END IF
370: *
371: ELSE IF( N0IN.EQ.( N0+2 ) ) THEN
372: *
373: * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
374: *
375: * Cases 10 and 11.
376: *
377: IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN
378: TTYPE = -10
379: S = THIRD*DMIN2
380: IF( Z( NN-5 ).GT.Z( NN-7 ) )
381: $ RETURN
382: B1 = Z( NN-5 ) / Z( NN-7 )
383: B2 = B1
384: IF( B2.EQ.ZERO )
385: $ GO TO 80
386: DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
387: IF( Z( I4 ).GT.Z( I4-2 ) )
388: $ RETURN
389: B1 = B1*( Z( I4 ) / Z( I4-2 ) )
390: B2 = B2 + B1
391: IF( HUNDRD*B1.LT.B2 )
392: $ GO TO 80
393: 70 CONTINUE
394: 80 CONTINUE
395: B2 = SQRT( CNST3*B2 )
396: A2 = DMIN2 / ( ONE+B2**2 )
397: GAP2 = Z( NN-7 ) + Z( NN-9 ) -
398: $ SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2
399: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
400: S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
401: ELSE
402: S = MAX( S, A2*( ONE-CNST2*B2 ) )
403: END IF
404: ELSE
405: S = QURTR*DMIN2
406: TTYPE = -11
407: END IF
408: ELSE IF( N0IN.GT.( N0+2 ) ) THEN
409: *
410: * Case 12, more than two eigenvalues deflated. No information.
411: *
412: S = ZERO
413: TTYPE = -12
414: END IF
415: *
416: TAU = S
417: RETURN
418: *
419: * End of DLASQ4
420: *
421: 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