![[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|
Return to dlasq4.f CVS log | Up to [local] / rpl / lapack / lapack |
1.12 bertrand 1: *> \brief \b DLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr.
1.9 bertrand 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*N )
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 REAL
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: *
1.12 bertrand 138: *> \date September 2012
1.9 bertrand 139: *
140: *> \ingroup auxOTHERcomputational
141: *
142: *> \par Further Details:
143: * =====================
144: *>
145: *> \verbatim
146: *>
147: *> CNST1 = 9/16
148: *> \endverbatim
149: *>
150: * =====================================================================
1.1 bertrand 151: SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
152: $ DN1, DN2, TAU, TTYPE, G )
153: *
1.12 bertrand 154: * -- LAPACK computational routine (version 3.4.2) --
1.1 bertrand 155: * -- LAPACK is a software package provided by Univ. of Tennessee, --
156: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.12 bertrand 157: * September 2012
1.1 bertrand 158: *
159: * .. Scalar Arguments ..
160: INTEGER I0, N0, N0IN, PP, TTYPE
161: DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
162: * ..
163: * .. Array Arguments ..
164: DOUBLE PRECISION Z( * )
165: * ..
166: *
167: * =====================================================================
168: *
169: * .. Parameters ..
170: DOUBLE PRECISION CNST1, CNST2, CNST3
171: PARAMETER ( CNST1 = 0.5630D0, CNST2 = 1.010D0,
172: $ CNST3 = 1.050D0 )
173: DOUBLE PRECISION QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
174: PARAMETER ( QURTR = 0.250D0, THIRD = 0.3330D0,
175: $ HALF = 0.50D0, ZERO = 0.0D0, ONE = 1.0D0,
176: $ TWO = 2.0D0, HUNDRD = 100.0D0 )
177: * ..
178: * .. Local Scalars ..
179: INTEGER I4, NN, NP
180: DOUBLE PRECISION A2, B1, B2, GAM, GAP1, GAP2, S
181: * ..
182: * .. Intrinsic Functions ..
183: INTRINSIC MAX, MIN, SQRT
184: * ..
185: * .. Executable Statements ..
186: *
187: * A negative DMIN forces the shift to take that absolute value
188: * TTYPE records the type of shift.
189: *
190: IF( DMIN.LE.ZERO ) THEN
191: TAU = -DMIN
192: TTYPE = -1
193: RETURN
194: END IF
195: *
196: NN = 4*N0 + PP
197: IF( N0IN.EQ.N0 ) THEN
198: *
199: * No eigenvalues deflated.
200: *
201: IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN
202: *
203: B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) )
204: B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) )
205: A2 = Z( NN-7 ) + Z( NN-5 )
206: *
207: * Cases 2 and 3.
208: *
209: IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN
210: GAP2 = DMIN2 - A2 - DMIN2*QURTR
211: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN
212: GAP1 = A2 - DN - ( B2 / GAP2 )*B2
213: ELSE
214: GAP1 = A2 - DN - ( B1+B2 )
215: END IF
216: IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN
217: S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN )
218: TTYPE = -2
219: ELSE
220: S = ZERO
221: IF( DN.GT.B1 )
222: $ S = DN - B1
223: IF( A2.GT.( B1+B2 ) )
224: $ S = MIN( S, A2-( B1+B2 ) )
225: S = MAX( S, THIRD*DMIN )
226: TTYPE = -3
227: END IF
228: ELSE
229: *
230: * Case 4.
231: *
232: TTYPE = -4
233: S = QURTR*DMIN
234: IF( DMIN.EQ.DN ) THEN
235: GAM = DN
236: A2 = ZERO
237: IF( Z( NN-5 ) .GT. Z( NN-7 ) )
238: $ RETURN
239: B2 = Z( NN-5 ) / Z( NN-7 )
240: NP = NN - 9
241: ELSE
242: NP = NN - 2*PP
243: B2 = Z( NP-2 )
244: GAM = DN1
245: IF( Z( NP-4 ) .GT. Z( NP-2 ) )
246: $ RETURN
247: A2 = Z( NP-4 ) / Z( NP-2 )
248: IF( Z( NN-9 ) .GT. Z( NN-11 ) )
249: $ RETURN
250: B2 = Z( NN-9 ) / Z( NN-11 )
251: NP = NN - 13
252: END IF
253: *
254: * Approximate contribution to norm squared from I < NN-1.
255: *
256: A2 = A2 + B2
257: DO 10 I4 = NP, 4*I0 - 1 + PP, -4
258: IF( B2.EQ.ZERO )
259: $ GO TO 20
260: B1 = B2
261: IF( Z( I4 ) .GT. Z( I4-2 ) )
262: $ RETURN
263: B2 = B2*( Z( I4 ) / Z( I4-2 ) )
264: A2 = A2 + B2
265: IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
266: $ GO TO 20
267: 10 CONTINUE
268: 20 CONTINUE
269: A2 = CNST3*A2
270: *
271: * Rayleigh quotient residual bound.
272: *
273: IF( A2.LT.CNST1 )
274: $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
275: END IF
276: ELSE IF( DMIN.EQ.DN2 ) THEN
277: *
278: * Case 5.
279: *
280: TTYPE = -5
281: S = QURTR*DMIN
282: *
283: * Compute contribution to norm squared from I > NN-2.
284: *
285: NP = NN - 2*PP
286: B1 = Z( NP-2 )
287: B2 = Z( NP-6 )
288: GAM = DN2
289: IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 )
290: $ RETURN
291: A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 )
292: *
293: * Approximate contribution to norm squared from I < NN-2.
294: *
295: IF( N0-I0.GT.2 ) THEN
296: B2 = Z( NN-13 ) / Z( NN-15 )
297: A2 = A2 + B2
298: DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4
299: IF( B2.EQ.ZERO )
300: $ GO TO 40
301: B1 = B2
302: IF( Z( I4 ) .GT. Z( I4-2 ) )
303: $ RETURN
304: B2 = B2*( Z( I4 ) / Z( I4-2 ) )
305: A2 = A2 + B2
306: IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
307: $ GO TO 40
308: 30 CONTINUE
309: 40 CONTINUE
310: A2 = CNST3*A2
311: END IF
312: *
313: IF( A2.LT.CNST1 )
314: $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
315: ELSE
316: *
317: * Case 6, no information to guide us.
318: *
319: IF( TTYPE.EQ.-6 ) THEN
320: G = G + THIRD*( ONE-G )
321: ELSE IF( TTYPE.EQ.-18 ) THEN
322: G = QURTR*THIRD
323: ELSE
324: G = QURTR
325: END IF
326: S = G*DMIN
327: TTYPE = -6
328: END IF
329: *
330: ELSE IF( N0IN.EQ.( N0+1 ) ) THEN
331: *
332: * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
333: *
334: IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN
335: *
336: * Cases 7 and 8.
337: *
338: TTYPE = -7
339: S = THIRD*DMIN1
340: IF( Z( NN-5 ).GT.Z( NN-7 ) )
341: $ RETURN
342: B1 = Z( NN-5 ) / Z( NN-7 )
343: B2 = B1
344: IF( B2.EQ.ZERO )
345: $ GO TO 60
346: DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
347: A2 = B1
348: IF( Z( I4 ).GT.Z( I4-2 ) )
349: $ RETURN
350: B1 = B1*( Z( I4 ) / Z( I4-2 ) )
351: B2 = B2 + B1
352: IF( HUNDRD*MAX( B1, A2 ).LT.B2 )
353: $ GO TO 60
354: 50 CONTINUE
355: 60 CONTINUE
356: B2 = SQRT( CNST3*B2 )
357: A2 = DMIN1 / ( ONE+B2**2 )
358: GAP2 = HALF*DMIN2 - A2
359: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
360: S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
361: ELSE
362: S = MAX( S, A2*( ONE-CNST2*B2 ) )
363: TTYPE = -8
364: END IF
365: ELSE
366: *
367: * Case 9.
368: *
369: S = QURTR*DMIN1
370: IF( DMIN1.EQ.DN1 )
371: $ S = HALF*DMIN1
372: TTYPE = -9
373: END IF
374: *
375: ELSE IF( N0IN.EQ.( N0+2 ) ) THEN
376: *
377: * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
378: *
379: * Cases 10 and 11.
380: *
381: IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN
382: TTYPE = -10
383: S = THIRD*DMIN2
384: IF( Z( NN-5 ).GT.Z( NN-7 ) )
385: $ RETURN
386: B1 = Z( NN-5 ) / Z( NN-7 )
387: B2 = B1
388: IF( B2.EQ.ZERO )
389: $ GO TO 80
390: DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
391: IF( Z( I4 ).GT.Z( I4-2 ) )
392: $ RETURN
393: B1 = B1*( Z( I4 ) / Z( I4-2 ) )
394: B2 = B2 + B1
395: IF( HUNDRD*B1.LT.B2 )
396: $ GO TO 80
397: 70 CONTINUE
398: 80 CONTINUE
399: B2 = SQRT( CNST3*B2 )
400: A2 = DMIN2 / ( ONE+B2**2 )
401: GAP2 = Z( NN-7 ) + Z( NN-9 ) -
402: $ SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2
403: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
404: S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
405: ELSE
406: S = MAX( S, A2*( ONE-CNST2*B2 ) )
407: END IF
408: ELSE
409: S = QURTR*DMIN2
410: TTYPE = -11
411: END IF
412: ELSE IF( N0IN.GT.( N0+2 ) ) THEN
413: *
414: * Case 12, more than two eigenvalues deflated. No information.
415: *
416: S = ZERO
417: TTYPE = -12
418: END IF
419: *
420: TAU = S
421: RETURN
422: *
423: * End of DLASQ4
424: *
425: END