1: *> \brief \b DLASQ3
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLASQ3 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq3.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq3.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq3.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
22: * ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
23: * DN2, G, TAU )
24: *
25: * .. Scalar Arguments ..
26: * LOGICAL IEEE
27: * INTEGER I0, ITER, N0, NDIV, NFAIL, PP
28: * DOUBLE PRECISION DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
29: * $ QMAX, SIGMA, TAU
30: * ..
31: * .. Array Arguments ..
32: * DOUBLE PRECISION Z( * )
33: * ..
34: *
35: *
36: *> \par Purpose:
37: * =============
38: *>
39: *> \verbatim
40: *>
41: *> DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
42: *> In case of failure it changes shifts, and tries again until output
43: *> is positive.
44: *> \endverbatim
45: *
46: * Arguments:
47: * ==========
48: *
49: *> \param[in] I0
50: *> \verbatim
51: *> I0 is INTEGER
52: *> First index.
53: *> \endverbatim
54: *>
55: *> \param[in,out] N0
56: *> \verbatim
57: *> N0 is INTEGER
58: *> Last index.
59: *> \endverbatim
60: *>
61: *> \param[in] Z
62: *> \verbatim
63: *> Z is DOUBLE PRECISION array, dimension ( 4*N )
64: *> Z holds the qd array.
65: *> \endverbatim
66: *>
67: *> \param[in,out] PP
68: *> \verbatim
69: *> PP is INTEGER
70: *> PP=0 for ping, PP=1 for pong.
71: *> PP=2 indicates that flipping was applied to the Z array
72: *> and that the initial tests for deflation should not be
73: *> performed.
74: *> \endverbatim
75: *>
76: *> \param[out] DMIN
77: *> \verbatim
78: *> DMIN is DOUBLE PRECISION
79: *> Minimum value of d.
80: *> \endverbatim
81: *>
82: *> \param[out] SIGMA
83: *> \verbatim
84: *> SIGMA is DOUBLE PRECISION
85: *> Sum of shifts used in current segment.
86: *> \endverbatim
87: *>
88: *> \param[in,out] DESIG
89: *> \verbatim
90: *> DESIG is DOUBLE PRECISION
91: *> Lower order part of SIGMA
92: *> \endverbatim
93: *>
94: *> \param[in] QMAX
95: *> \verbatim
96: *> QMAX is DOUBLE PRECISION
97: *> Maximum value of q.
98: *> \endverbatim
99: *>
100: *> \param[out] NFAIL
101: *> \verbatim
102: *> NFAIL is INTEGER
103: *> Number of times shift was too big.
104: *> \endverbatim
105: *>
106: *> \param[out] ITER
107: *> \verbatim
108: *> ITER is INTEGER
109: *> Number of iterations.
110: *> \endverbatim
111: *>
112: *> \param[out] NDIV
113: *> \verbatim
114: *> NDIV is INTEGER
115: *> Number of divisions.
116: *> \endverbatim
117: *>
118: *> \param[in] IEEE
119: *> \verbatim
120: *> IEEE is LOGICAL
121: *> Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
122: *> \endverbatim
123: *>
124: *> \param[in,out] TTYPE
125: *> \verbatim
126: *> TTYPE is INTEGER
127: *> Shift type.
128: *> \endverbatim
129: *>
130: *> \param[in,out] DMIN1
131: *> \verbatim
132: *> DMIN1 is DOUBLE PRECISION
133: *> \endverbatim
134: *>
135: *> \param[in,out] DMIN2
136: *> \verbatim
137: *> DMIN2 is DOUBLE PRECISION
138: *> \endverbatim
139: *>
140: *> \param[in,out] DN
141: *> \verbatim
142: *> DN is DOUBLE PRECISION
143: *> \endverbatim
144: *>
145: *> \param[in,out] DN1
146: *> \verbatim
147: *> DN1 is DOUBLE PRECISION
148: *> \endverbatim
149: *>
150: *> \param[in,out] DN2
151: *> \verbatim
152: *> DN2 is DOUBLE PRECISION
153: *> \endverbatim
154: *>
155: *> \param[in,out] G
156: *> \verbatim
157: *> G is DOUBLE PRECISION
158: *> \endverbatim
159: *>
160: *> \param[in,out] TAU
161: *> \verbatim
162: *> TAU is DOUBLE PRECISION
163: *>
164: *> These are passed as arguments in order to save their values
165: *> between calls to DLASQ3.
166: *> \endverbatim
167: *
168: * Authors:
169: * ========
170: *
171: *> \author Univ. of Tennessee
172: *> \author Univ. of California Berkeley
173: *> \author Univ. of Colorado Denver
174: *> \author NAG Ltd.
175: *
176: *> \date November 2011
177: *
178: *> \ingroup auxOTHERcomputational
179: *
180: * =====================================================================
181: SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
182: $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
183: $ DN2, G, TAU )
184: *
185: * -- LAPACK computational routine (version 3.4.0) --
186: * -- LAPACK is a software package provided by Univ. of Tennessee, --
187: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188: * November 2011
189: *
190: * .. Scalar Arguments ..
191: LOGICAL IEEE
192: INTEGER I0, ITER, N0, NDIV, NFAIL, PP
193: DOUBLE PRECISION DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
194: $ QMAX, SIGMA, TAU
195: * ..
196: * .. Array Arguments ..
197: DOUBLE PRECISION Z( * )
198: * ..
199: *
200: * =====================================================================
201: *
202: * .. Parameters ..
203: DOUBLE PRECISION CBIAS
204: PARAMETER ( CBIAS = 1.50D0 )
205: DOUBLE PRECISION ZERO, QURTR, HALF, ONE, TWO, HUNDRD
206: PARAMETER ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0,
207: $ ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 )
208: * ..
209: * .. Local Scalars ..
210: INTEGER IPN4, J4, N0IN, NN, TTYPE
211: DOUBLE PRECISION EPS, S, T, TEMP, TOL, TOL2
212: * ..
213: * .. External Subroutines ..
214: EXTERNAL DLASQ4, DLASQ5, DLASQ6
215: * ..
216: * .. External Function ..
217: DOUBLE PRECISION DLAMCH
218: LOGICAL DISNAN
219: EXTERNAL DISNAN, DLAMCH
220: * ..
221: * .. Intrinsic Functions ..
222: INTRINSIC ABS, MAX, MIN, SQRT
223: * ..
224: * .. Executable Statements ..
225: *
226: N0IN = N0
227: EPS = DLAMCH( 'Precision' )
228: TOL = EPS*HUNDRD
229: TOL2 = TOL**2
230: *
231: * Check for deflation.
232: *
233: 10 CONTINUE
234: *
235: IF( N0.LT.I0 )
236: $ RETURN
237: IF( N0.EQ.I0 )
238: $ GO TO 20
239: NN = 4*N0 + PP
240: IF( N0.EQ.( I0+1 ) )
241: $ GO TO 40
242: *
243: * Check whether E(N0-1) is negligible, 1 eigenvalue.
244: *
245: IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
246: $ Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
247: $ GO TO 30
248: *
249: 20 CONTINUE
250: *
251: Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
252: N0 = N0 - 1
253: GO TO 10
254: *
255: * Check whether E(N0-2) is negligible, 2 eigenvalues.
256: *
257: 30 CONTINUE
258: *
259: IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
260: $ Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
261: $ GO TO 50
262: *
263: 40 CONTINUE
264: *
265: IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
266: S = Z( NN-3 )
267: Z( NN-3 ) = Z( NN-7 )
268: Z( NN-7 ) = S
269: END IF
270: IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN
271: T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
272: S = Z( NN-3 )*( Z( NN-5 ) / T )
273: IF( S.LE.T ) THEN
274: S = Z( NN-3 )*( Z( NN-5 ) /
275: $ ( T*( ONE+SQRT( ONE+S / T ) ) ) )
276: ELSE
277: S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
278: END IF
279: T = Z( NN-7 ) + ( S+Z( NN-5 ) )
280: Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
281: Z( NN-7 ) = T
282: END IF
283: Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
284: Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
285: N0 = N0 - 2
286: GO TO 10
287: *
288: 50 CONTINUE
289: IF( PP.EQ.2 )
290: $ PP = 0
291: *
292: * Reverse the qd-array, if warranted.
293: *
294: IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
295: IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
296: IPN4 = 4*( I0+N0 )
297: DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
298: TEMP = Z( J4-3 )
299: Z( J4-3 ) = Z( IPN4-J4-3 )
300: Z( IPN4-J4-3 ) = TEMP
301: TEMP = Z( J4-2 )
302: Z( J4-2 ) = Z( IPN4-J4-2 )
303: Z( IPN4-J4-2 ) = TEMP
304: TEMP = Z( J4-1 )
305: Z( J4-1 ) = Z( IPN4-J4-5 )
306: Z( IPN4-J4-5 ) = TEMP
307: TEMP = Z( J4 )
308: Z( J4 ) = Z( IPN4-J4-4 )
309: Z( IPN4-J4-4 ) = TEMP
310: 60 CONTINUE
311: IF( N0-I0.LE.4 ) THEN
312: Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
313: Z( 4*N0-PP ) = Z( 4*I0-PP )
314: END IF
315: DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
316: Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
317: $ Z( 4*I0+PP+3 ) )
318: Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
319: $ Z( 4*I0-PP+4 ) )
320: QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
321: DMIN = -ZERO
322: END IF
323: END IF
324: *
325: * Choose a shift.
326: *
327: CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
328: $ DN2, TAU, TTYPE, G )
329: *
330: * Call dqds until DMIN > 0.
331: *
332: 70 CONTINUE
333: *
334: CALL DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
335: $ DN1, DN2, IEEE )
336: *
337: NDIV = NDIV + ( N0-I0+2 )
338: ITER = ITER + 1
339: *
340: * Check status.
341: *
342: IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN
343: *
344: * Success.
345: *
346: GO TO 90
347: *
348: ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND.
349: $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
350: $ ABS( DN ).LT.TOL*SIGMA ) THEN
351: *
352: * Convergence hidden by negative DN.
353: *
354: Z( 4*( N0-1 )-PP+2 ) = ZERO
355: DMIN = ZERO
356: GO TO 90
357: ELSE IF( DMIN.LT.ZERO ) THEN
358: *
359: * TAU too big. Select new TAU and try again.
360: *
361: NFAIL = NFAIL + 1
362: IF( TTYPE.LT.-22 ) THEN
363: *
364: * Failed twice. Play it safe.
365: *
366: TAU = ZERO
367: ELSE IF( DMIN1.GT.ZERO ) THEN
368: *
369: * Late failure. Gives excellent shift.
370: *
371: TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
372: TTYPE = TTYPE - 11
373: ELSE
374: *
375: * Early failure. Divide by 4.
376: *
377: TAU = QURTR*TAU
378: TTYPE = TTYPE - 12
379: END IF
380: GO TO 70
381: ELSE IF( DISNAN( DMIN ) ) THEN
382: *
383: * NaN.
384: *
385: IF( TAU.EQ.ZERO ) THEN
386: GO TO 80
387: ELSE
388: TAU = ZERO
389: GO TO 70
390: END IF
391: ELSE
392: *
393: * Possible underflow. Play it safe.
394: *
395: GO TO 80
396: END IF
397: *
398: * Risk of underflow.
399: *
400: 80 CONTINUE
401: CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
402: NDIV = NDIV + ( N0-I0+2 )
403: ITER = ITER + 1
404: TAU = ZERO
405: *
406: 90 CONTINUE
407: IF( TAU.LT.SIGMA ) THEN
408: DESIG = DESIG + TAU
409: T = SIGMA + DESIG
410: DESIG = DESIG - ( T-SIGMA )
411: ELSE
412: T = SIGMA + TAU
413: DESIG = SIGMA - ( T-TAU ) + DESIG
414: END IF
415: SIGMA = T
416: *
417: RETURN
418: *
419: * End of DLASQ3
420: *
421: END
CVSweb interface <joel.bertrand@systella.fr>