1: *> \brief \b DLASQ3 checks for deflation, computes a shift and calls dqds. 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 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,out] Z
62: *> \verbatim
63: *> Z is DOUBLE PRECISION array, dimension ( 4*N0 )
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[in,out] NFAIL
101: *> \verbatim
102: *> NFAIL is INTEGER
103: *> Increment NFAIL by 1 each time the shift was too big.
104: *> \endverbatim
105: *>
106: *> \param[in,out] ITER
107: *> \verbatim
108: *> ITER is INTEGER
109: *> Increment ITER by 1 for each iteration.
110: *> \endverbatim
111: *>
112: *> \param[in,out] NDIV
113: *> \verbatim
114: *> NDIV is INTEGER
115: *> Increment NDIV by 1 for each division.
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: *> \ingroup auxOTHERcomputational
177: *
178: * =====================================================================
179: SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
180: $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
181: $ DN2, G, TAU )
182: *
183: * -- LAPACK computational routine --
184: * -- LAPACK is a software package provided by Univ. of Tennessee, --
185: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186: *
187: * .. Scalar Arguments ..
188: LOGICAL IEEE
189: INTEGER I0, ITER, N0, NDIV, NFAIL, PP
190: DOUBLE PRECISION DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
191: $ QMAX, SIGMA, TAU
192: * ..
193: * .. Array Arguments ..
194: DOUBLE PRECISION Z( * )
195: * ..
196: *
197: * =====================================================================
198: *
199: * .. Parameters ..
200: DOUBLE PRECISION CBIAS
201: PARAMETER ( CBIAS = 1.50D0 )
202: DOUBLE PRECISION ZERO, QURTR, HALF, ONE, TWO, HUNDRD
203: PARAMETER ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0,
204: $ ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 )
205: * ..
206: * .. Local Scalars ..
207: INTEGER IPN4, J4, N0IN, NN, TTYPE
208: DOUBLE PRECISION EPS, S, T, TEMP, TOL, TOL2
209: * ..
210: * .. External Subroutines ..
211: EXTERNAL DLASQ4, DLASQ5, DLASQ6
212: * ..
213: * .. External Function ..
214: DOUBLE PRECISION DLAMCH
215: LOGICAL DISNAN
216: EXTERNAL DISNAN, DLAMCH
217: * ..
218: * .. Intrinsic Functions ..
219: INTRINSIC ABS, MAX, MIN, SQRT
220: * ..
221: * .. Executable Statements ..
222: *
223: N0IN = N0
224: EPS = DLAMCH( 'Precision' )
225: TOL = EPS*HUNDRD
226: TOL2 = TOL**2
227: *
228: * Check for deflation.
229: *
230: 10 CONTINUE
231: *
232: IF( N0.LT.I0 )
233: $ RETURN
234: IF( N0.EQ.I0 )
235: $ GO TO 20
236: NN = 4*N0 + PP
237: IF( N0.EQ.( I0+1 ) )
238: $ GO TO 40
239: *
240: * Check whether E(N0-1) is negligible, 1 eigenvalue.
241: *
242: IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
243: $ Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
244: $ GO TO 30
245: *
246: 20 CONTINUE
247: *
248: Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
249: N0 = N0 - 1
250: GO TO 10
251: *
252: * Check whether E(N0-2) is negligible, 2 eigenvalues.
253: *
254: 30 CONTINUE
255: *
256: IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
257: $ Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
258: $ GO TO 50
259: *
260: 40 CONTINUE
261: *
262: IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
263: S = Z( NN-3 )
264: Z( NN-3 ) = Z( NN-7 )
265: Z( NN-7 ) = S
266: END IF
267: T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
268: IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2.AND.T.NE.ZERO ) THEN
269: S = Z( NN-3 )*( Z( NN-5 ) / T )
270: IF( S.LE.T ) THEN
271: S = Z( NN-3 )*( Z( NN-5 ) /
272: $ ( T*( ONE+SQRT( ONE+S / T ) ) ) )
273: ELSE
274: S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
275: END IF
276: T = Z( NN-7 ) + ( S+Z( NN-5 ) )
277: Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
278: Z( NN-7 ) = T
279: END IF
280: Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
281: Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
282: N0 = N0 - 2
283: GO TO 10
284: *
285: 50 CONTINUE
286: IF( PP.EQ.2 )
287: $ PP = 0
288: *
289: * Reverse the qd-array, if warranted.
290: *
291: IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
292: IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
293: IPN4 = 4*( I0+N0 )
294: DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
295: TEMP = Z( J4-3 )
296: Z( J4-3 ) = Z( IPN4-J4-3 )
297: Z( IPN4-J4-3 ) = TEMP
298: TEMP = Z( J4-2 )
299: Z( J4-2 ) = Z( IPN4-J4-2 )
300: Z( IPN4-J4-2 ) = TEMP
301: TEMP = Z( J4-1 )
302: Z( J4-1 ) = Z( IPN4-J4-5 )
303: Z( IPN4-J4-5 ) = TEMP
304: TEMP = Z( J4 )
305: Z( J4 ) = Z( IPN4-J4-4 )
306: Z( IPN4-J4-4 ) = TEMP
307: 60 CONTINUE
308: IF( N0-I0.LE.4 ) THEN
309: Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
310: Z( 4*N0-PP ) = Z( 4*I0-PP )
311: END IF
312: DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
313: Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
314: $ Z( 4*I0+PP+3 ) )
315: Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
316: $ Z( 4*I0-PP+4 ) )
317: QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
318: DMIN = -ZERO
319: END IF
320: END IF
321: *
322: * Choose a shift.
323: *
324: CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
325: $ DN2, TAU, TTYPE, G )
326: *
327: * Call dqds until DMIN > 0.
328: *
329: 70 CONTINUE
330: *
331: CALL DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
332: $ DN1, DN2, IEEE, EPS )
333: *
334: NDIV = NDIV + ( N0-I0+2 )
335: ITER = ITER + 1
336: *
337: * Check status.
338: *
339: IF( DMIN.GE.ZERO .AND. DMIN1.GE.ZERO ) THEN
340: *
341: * Success.
342: *
343: GO TO 90
344: *
345: ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND.
346: $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
347: $ ABS( DN ).LT.TOL*SIGMA ) THEN
348: *
349: * Convergence hidden by negative DN.
350: *
351: Z( 4*( N0-1 )-PP+2 ) = ZERO
352: DMIN = ZERO
353: GO TO 90
354: ELSE IF( DMIN.LT.ZERO ) THEN
355: *
356: * TAU too big. Select new TAU and try again.
357: *
358: NFAIL = NFAIL + 1
359: IF( TTYPE.LT.-22 ) THEN
360: *
361: * Failed twice. Play it safe.
362: *
363: TAU = ZERO
364: ELSE IF( DMIN1.GT.ZERO ) THEN
365: *
366: * Late failure. Gives excellent shift.
367: *
368: TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
369: TTYPE = TTYPE - 11
370: ELSE
371: *
372: * Early failure. Divide by 4.
373: *
374: TAU = QURTR*TAU
375: TTYPE = TTYPE - 12
376: END IF
377: GO TO 70
378: ELSE IF( DISNAN( DMIN ) ) THEN
379: *
380: * NaN.
381: *
382: IF( TAU.EQ.ZERO ) THEN
383: GO TO 80
384: ELSE
385: TAU = ZERO
386: GO TO 70
387: END IF
388: ELSE
389: *
390: * Possible underflow. Play it safe.
391: *
392: GO TO 80
393: END IF
394: *
395: * Risk of underflow.
396: *
397: 80 CONTINUE
398: CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
399: NDIV = NDIV + ( N0-I0+2 )
400: ITER = ITER + 1
401: TAU = ZERO
402: *
403: 90 CONTINUE
404: IF( TAU.LT.SIGMA ) THEN
405: DESIG = DESIG + TAU
406: T = SIGMA + DESIG
407: DESIG = DESIG - ( T-SIGMA )
408: ELSE
409: T = SIGMA + TAU
410: DESIG = SIGMA - ( T-TAU ) + DESIG
411: END IF
412: SIGMA = T
413: *
414: RETURN
415: *
416: * End of DLASQ3
417: *
418: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack