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: *> \date June 2016
139: *
140: *> \ingroup auxOTHERcomputational
141: *
142: *> \par Further Details:
143: * =====================
144: *>
145: *> \verbatim
146: *>
147: *> CNST1 = 9/16
148: *> \endverbatim
149: *>
150: * =====================================================================
151: SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
152: $ DN1, DN2, TAU, TTYPE, G )
153: *
154: * -- LAPACK computational routine (version 3.7.1) --
155: * -- LAPACK is a software package provided by Univ. of Tennessee, --
156: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
157: * June 2016
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: GAM = DN1
244: IF( Z( NP-4 ) .GT. Z( NP-2 ) )
245: $ RETURN
246: A2 = Z( NP-4 ) / Z( NP-2 )
247: IF( Z( NN-9 ) .GT. Z( NN-11 ) )
248: $ RETURN
249: B2 = Z( NN-9 ) / Z( NN-11 )
250: NP = NN - 13
251: END IF
252: *
253: * Approximate contribution to norm squared from I < NN-1.
254: *
255: A2 = A2 + B2
256: DO 10 I4 = NP, 4*I0 - 1 + PP, -4
257: IF( B2.EQ.ZERO )
258: $ GO TO 20
259: B1 = B2
260: IF( Z( I4 ) .GT. Z( I4-2 ) )
261: $ RETURN
262: B2 = B2*( Z( I4 ) / Z( I4-2 ) )
263: A2 = A2 + B2
264: IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
265: $ GO TO 20
266: 10 CONTINUE
267: 20 CONTINUE
268: A2 = CNST3*A2
269: *
270: * Rayleigh quotient residual bound.
271: *
272: IF( A2.LT.CNST1 )
273: $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
274: END IF
275: ELSE IF( DMIN.EQ.DN2 ) THEN
276: *
277: * Case 5.
278: *
279: TTYPE = -5
280: S = QURTR*DMIN
281: *
282: * Compute contribution to norm squared from I > NN-2.
283: *
284: NP = NN - 2*PP
285: B1 = Z( NP-2 )
286: B2 = Z( NP-6 )
287: GAM = DN2
288: IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 )
289: $ RETURN
290: A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 )
291: *
292: * Approximate contribution to norm squared from I < NN-2.
293: *
294: IF( N0-I0.GT.2 ) THEN
295: B2 = Z( NN-13 ) / Z( NN-15 )
296: A2 = A2 + B2
297: DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4
298: IF( B2.EQ.ZERO )
299: $ GO TO 40
300: B1 = B2
301: IF( Z( I4 ) .GT. Z( I4-2 ) )
302: $ RETURN
303: B2 = B2*( Z( I4 ) / Z( I4-2 ) )
304: A2 = A2 + B2
305: IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
306: $ GO TO 40
307: 30 CONTINUE
308: 40 CONTINUE
309: A2 = CNST3*A2
310: END IF
311: *
312: IF( A2.LT.CNST1 )
313: $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
314: ELSE
315: *
316: * Case 6, no information to guide us.
317: *
318: IF( TTYPE.EQ.-6 ) THEN
319: G = G + THIRD*( ONE-G )
320: ELSE IF( TTYPE.EQ.-18 ) THEN
321: G = QURTR*THIRD
322: ELSE
323: G = QURTR
324: END IF
325: S = G*DMIN
326: TTYPE = -6
327: END IF
328: *
329: ELSE IF( N0IN.EQ.( N0+1 ) ) THEN
330: *
331: * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
332: *
333: IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN
334: *
335: * Cases 7 and 8.
336: *
337: TTYPE = -7
338: S = THIRD*DMIN1
339: IF( Z( NN-5 ).GT.Z( NN-7 ) )
340: $ RETURN
341: B1 = Z( NN-5 ) / Z( NN-7 )
342: B2 = B1
343: IF( B2.EQ.ZERO )
344: $ GO TO 60
345: DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
346: A2 = B1
347: IF( Z( I4 ).GT.Z( I4-2 ) )
348: $ RETURN
349: B1 = B1*( Z( I4 ) / Z( I4-2 ) )
350: B2 = B2 + B1
351: IF( HUNDRD*MAX( B1, A2 ).LT.B2 )
352: $ GO TO 60
353: 50 CONTINUE
354: 60 CONTINUE
355: B2 = SQRT( CNST3*B2 )
356: A2 = DMIN1 / ( ONE+B2**2 )
357: GAP2 = HALF*DMIN2 - A2
358: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
359: S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
360: ELSE
361: S = MAX( S, A2*( ONE-CNST2*B2 ) )
362: TTYPE = -8
363: END IF
364: ELSE
365: *
366: * Case 9.
367: *
368: S = QURTR*DMIN1
369: IF( DMIN1.EQ.DN1 )
370: $ S = HALF*DMIN1
371: TTYPE = -9
372: END IF
373: *
374: ELSE IF( N0IN.EQ.( N0+2 ) ) THEN
375: *
376: * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
377: *
378: * Cases 10 and 11.
379: *
380: IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN
381: TTYPE = -10
382: S = THIRD*DMIN2
383: IF( Z( NN-5 ).GT.Z( NN-7 ) )
384: $ RETURN
385: B1 = Z( NN-5 ) / Z( NN-7 )
386: B2 = B1
387: IF( B2.EQ.ZERO )
388: $ GO TO 80
389: DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
390: IF( Z( I4 ).GT.Z( I4-2 ) )
391: $ RETURN
392: B1 = B1*( Z( I4 ) / Z( I4-2 ) )
393: B2 = B2 + B1
394: IF( HUNDRD*B1.LT.B2 )
395: $ GO TO 80
396: 70 CONTINUE
397: 80 CONTINUE
398: B2 = SQRT( CNST3*B2 )
399: A2 = DMIN2 / ( ONE+B2**2 )
400: GAP2 = Z( NN-7 ) + Z( NN-9 ) -
401: $ SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2
402: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
403: S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
404: ELSE
405: S = MAX( S, A2*( ONE-CNST2*B2 ) )
406: END IF
407: ELSE
408: S = QURTR*DMIN2
409: TTYPE = -11
410: END IF
411: ELSE IF( N0IN.GT.( N0+2 ) ) THEN
412: *
413: * Case 12, more than two eigenvalues deflated. No information.
414: *
415: S = ZERO
416: TTYPE = -12
417: END IF
418: *
419: TAU = S
420: RETURN
421: *
422: * End of DLASQ4
423: *
424: END
CVSweb interface <joel.bertrand@systella.fr>