1: SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
2: $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
3: $ DN2, G, TAU )
4: *
5: * -- LAPACK routine (version 3.2) --
6: *
7: * -- Contributed by Osni Marques of the Lawrence Berkeley National --
8: * -- Laboratory and Beresford Parlett of the Univ. of California at --
9: * -- Berkeley --
10: * -- November 2008 --
11: *
12: * -- LAPACK is a software package provided by Univ. of Tennessee, --
13: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
14: *
15: * .. Scalar Arguments ..
16: LOGICAL IEEE
17: INTEGER I0, ITER, N0, NDIV, NFAIL, PP
18: DOUBLE PRECISION DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
19: $ QMAX, SIGMA, TAU
20: * ..
21: * .. Array Arguments ..
22: DOUBLE PRECISION Z( * )
23: * ..
24: *
25: * Purpose
26: * =======
27: *
28: * DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
29: * In case of failure it changes shifts, and tries again until output
30: * is positive.
31: *
32: * Arguments
33: * =========
34: *
35: * I0 (input) INTEGER
36: * First index.
37: *
38: * N0 (input) INTEGER
39: * Last index.
40: *
41: * Z (input) DOUBLE PRECISION array, dimension ( 4*N )
42: * Z holds the qd array.
43: *
44: * PP (input/output) INTEGER
45: * PP=0 for ping, PP=1 for pong.
46: * PP=2 indicates that flipping was applied to the Z array
47: * and that the initial tests for deflation should not be
48: * performed.
49: *
50: * DMIN (output) DOUBLE PRECISION
51: * Minimum value of d.
52: *
53: * SIGMA (output) DOUBLE PRECISION
54: * Sum of shifts used in current segment.
55: *
56: * DESIG (input/output) DOUBLE PRECISION
57: * Lower order part of SIGMA
58: *
59: * QMAX (input) DOUBLE PRECISION
60: * Maximum value of q.
61: *
62: * NFAIL (output) INTEGER
63: * Number of times shift was too big.
64: *
65: * ITER (output) INTEGER
66: * Number of iterations.
67: *
68: * NDIV (output) INTEGER
69: * Number of divisions.
70: *
71: * IEEE (input) LOGICAL
72: * Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
73: *
74: * TTYPE (input/output) INTEGER
75: * Shift type.
76: *
77: * DMIN1, DMIN2, DN, DN1, DN2, G, TAU (input/output) DOUBLE PRECISION
78: * These are passed as arguments in order to save their values
79: * between calls to DLASQ3.
80: *
81: * =====================================================================
82: *
83: * .. Parameters ..
84: DOUBLE PRECISION CBIAS
85: PARAMETER ( CBIAS = 1.50D0 )
86: DOUBLE PRECISION ZERO, QURTR, HALF, ONE, TWO, HUNDRD
87: PARAMETER ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0,
88: $ ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 )
89: * ..
90: * .. Local Scalars ..
91: INTEGER IPN4, J4, N0IN, NN, TTYPE
92: DOUBLE PRECISION EPS, S, T, TEMP, TOL, TOL2
93: * ..
94: * .. External Subroutines ..
95: EXTERNAL DLASQ4, DLASQ5, DLASQ6
96: * ..
97: * .. External Function ..
98: DOUBLE PRECISION DLAMCH
99: LOGICAL DISNAN
100: EXTERNAL DISNAN, DLAMCH
101: * ..
102: * .. Intrinsic Functions ..
103: INTRINSIC ABS, MAX, MIN, SQRT
104: * ..
105: * .. Executable Statements ..
106: *
107: N0IN = N0
108: EPS = DLAMCH( 'Precision' )
109: TOL = EPS*HUNDRD
110: TOL2 = TOL**2
111: *
112: * Check for deflation.
113: *
114: 10 CONTINUE
115: *
116: IF( N0.LT.I0 )
117: $ RETURN
118: IF( N0.EQ.I0 )
119: $ GO TO 20
120: NN = 4*N0 + PP
121: IF( N0.EQ.( I0+1 ) )
122: $ GO TO 40
123: *
124: * Check whether E(N0-1) is negligible, 1 eigenvalue.
125: *
126: IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
127: $ Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
128: $ GO TO 30
129: *
130: 20 CONTINUE
131: *
132: Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
133: N0 = N0 - 1
134: GO TO 10
135: *
136: * Check whether E(N0-2) is negligible, 2 eigenvalues.
137: *
138: 30 CONTINUE
139: *
140: IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
141: $ Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
142: $ GO TO 50
143: *
144: 40 CONTINUE
145: *
146: IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
147: S = Z( NN-3 )
148: Z( NN-3 ) = Z( NN-7 )
149: Z( NN-7 ) = S
150: END IF
151: IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN
152: T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
153: S = Z( NN-3 )*( Z( NN-5 ) / T )
154: IF( S.LE.T ) THEN
155: S = Z( NN-3 )*( Z( NN-5 ) /
156: $ ( T*( ONE+SQRT( ONE+S / T ) ) ) )
157: ELSE
158: S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
159: END IF
160: T = Z( NN-7 ) + ( S+Z( NN-5 ) )
161: Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
162: Z( NN-7 ) = T
163: END IF
164: Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
165: Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
166: N0 = N0 - 2
167: GO TO 10
168: *
169: 50 CONTINUE
170: IF( PP.EQ.2 )
171: $ PP = 0
172: *
173: * Reverse the qd-array, if warranted.
174: *
175: IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
176: IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
177: IPN4 = 4*( I0+N0 )
178: DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
179: TEMP = Z( J4-3 )
180: Z( J4-3 ) = Z( IPN4-J4-3 )
181: Z( IPN4-J4-3 ) = TEMP
182: TEMP = Z( J4-2 )
183: Z( J4-2 ) = Z( IPN4-J4-2 )
184: Z( IPN4-J4-2 ) = TEMP
185: TEMP = Z( J4-1 )
186: Z( J4-1 ) = Z( IPN4-J4-5 )
187: Z( IPN4-J4-5 ) = TEMP
188: TEMP = Z( J4 )
189: Z( J4 ) = Z( IPN4-J4-4 )
190: Z( IPN4-J4-4 ) = TEMP
191: 60 CONTINUE
192: IF( N0-I0.LE.4 ) THEN
193: Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
194: Z( 4*N0-PP ) = Z( 4*I0-PP )
195: END IF
196: DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
197: Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
198: $ Z( 4*I0+PP+3 ) )
199: Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
200: $ Z( 4*I0-PP+4 ) )
201: QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
202: DMIN = -ZERO
203: END IF
204: END IF
205: *
206: * Choose a shift.
207: *
208: CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
209: $ DN2, TAU, TTYPE, G )
210: *
211: * Call dqds until DMIN > 0.
212: *
213: 70 CONTINUE
214: *
215: CALL DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
216: $ DN1, DN2, IEEE )
217: *
218: NDIV = NDIV + ( N0-I0+2 )
219: ITER = ITER + 1
220: *
221: * Check status.
222: *
223: IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN
224: *
225: * Success.
226: *
227: GO TO 90
228: *
229: ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND.
230: $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
231: $ ABS( DN ).LT.TOL*SIGMA ) THEN
232: *
233: * Convergence hidden by negative DN.
234: *
235: Z( 4*( N0-1 )-PP+2 ) = ZERO
236: DMIN = ZERO
237: GO TO 90
238: ELSE IF( DMIN.LT.ZERO ) THEN
239: *
240: * TAU too big. Select new TAU and try again.
241: *
242: NFAIL = NFAIL + 1
243: IF( TTYPE.LT.-22 ) THEN
244: *
245: * Failed twice. Play it safe.
246: *
247: TAU = ZERO
248: ELSE IF( DMIN1.GT.ZERO ) THEN
249: *
250: * Late failure. Gives excellent shift.
251: *
252: TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
253: TTYPE = TTYPE - 11
254: ELSE
255: *
256: * Early failure. Divide by 4.
257: *
258: TAU = QURTR*TAU
259: TTYPE = TTYPE - 12
260: END IF
261: GO TO 70
262: ELSE IF( DISNAN( DMIN ) ) THEN
263: *
264: * NaN.
265: *
266: IF( TAU.EQ.ZERO ) THEN
267: GO TO 80
268: ELSE
269: TAU = ZERO
270: GO TO 70
271: END IF
272: ELSE
273: *
274: * Possible underflow. Play it safe.
275: *
276: GO TO 80
277: END IF
278: *
279: * Risk of underflow.
280: *
281: 80 CONTINUE
282: CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
283: NDIV = NDIV + ( N0-I0+2 )
284: ITER = ITER + 1
285: TAU = ZERO
286: *
287: 90 CONTINUE
288: IF( TAU.LT.SIGMA ) THEN
289: DESIG = DESIG + TAU
290: T = SIGMA + DESIG
291: DESIG = DESIG - ( T-SIGMA )
292: ELSE
293: T = SIGMA + TAU
294: DESIG = SIGMA - ( T-TAU ) + DESIG
295: END IF
296: SIGMA = T
297: *
298: RETURN
299: *
300: * End of DLASQ3
301: *
302: 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