Annotation of rpl/lapack/lapack/dlasq4.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLASQ4
! 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: *
! 138: *> \date November 2011
! 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.9 ! bertrand 154: * -- LAPACK computational routine (version 3.4.0) --
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.9 ! bertrand 157: * November 2011
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
CVSweb interface <joel.bertrand@systella.fr>