![[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|
Return to dlasq3.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 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: *
1.5 bertrand 5: * -- LAPACK routine (version 3.2.2) --
1.1 bertrand 6: *
7: * -- Contributed by Osni Marques of the Lawrence Berkeley National --
8: * -- Laboratory and Beresford Parlett of the Univ. of California at --
9: * -- Berkeley --
1.5 bertrand 10: * -- June 2010 --
1.1 bertrand 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: *
1.5 bertrand 38: * N0 (input/output) INTEGER
1.1 bertrand 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: *
1.5 bertrand 77: * DMIN1 (input/output) DOUBLE PRECISION
78: *
79: * DMIN2 (input/output) DOUBLE PRECISION
80: *
81: * DN (input/output) DOUBLE PRECISION
82: *
83: * DN1 (input/output) DOUBLE PRECISION
84: *
85: * DN2 (input/output) DOUBLE PRECISION
86: *
87: * G (input/output) DOUBLE PRECISION
88: *
89: * TAU (input/output) DOUBLE PRECISION
90: *
1.1 bertrand 91: * These are passed as arguments in order to save their values
92: * between calls to DLASQ3.
93: *
94: * =====================================================================
95: *
96: * .. Parameters ..
97: DOUBLE PRECISION CBIAS
98: PARAMETER ( CBIAS = 1.50D0 )
99: DOUBLE PRECISION ZERO, QURTR, HALF, ONE, TWO, HUNDRD
100: PARAMETER ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0,
101: $ ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 )
102: * ..
103: * .. Local Scalars ..
104: INTEGER IPN4, J4, N0IN, NN, TTYPE
105: DOUBLE PRECISION EPS, S, T, TEMP, TOL, TOL2
106: * ..
107: * .. External Subroutines ..
108: EXTERNAL DLASQ4, DLASQ5, DLASQ6
109: * ..
110: * .. External Function ..
111: DOUBLE PRECISION DLAMCH
112: LOGICAL DISNAN
113: EXTERNAL DISNAN, DLAMCH
114: * ..
115: * .. Intrinsic Functions ..
116: INTRINSIC ABS, MAX, MIN, SQRT
117: * ..
118: * .. Executable Statements ..
119: *
120: N0IN = N0
121: EPS = DLAMCH( 'Precision' )
122: TOL = EPS*HUNDRD
123: TOL2 = TOL**2
124: *
125: * Check for deflation.
126: *
127: 10 CONTINUE
128: *
129: IF( N0.LT.I0 )
130: $ RETURN
131: IF( N0.EQ.I0 )
132: $ GO TO 20
133: NN = 4*N0 + PP
134: IF( N0.EQ.( I0+1 ) )
135: $ GO TO 40
136: *
137: * Check whether E(N0-1) is negligible, 1 eigenvalue.
138: *
139: IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
140: $ Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
141: $ GO TO 30
142: *
143: 20 CONTINUE
144: *
145: Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
146: N0 = N0 - 1
147: GO TO 10
148: *
149: * Check whether E(N0-2) is negligible, 2 eigenvalues.
150: *
151: 30 CONTINUE
152: *
153: IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
154: $ Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
155: $ GO TO 50
156: *
157: 40 CONTINUE
158: *
159: IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
160: S = Z( NN-3 )
161: Z( NN-3 ) = Z( NN-7 )
162: Z( NN-7 ) = S
163: END IF
164: IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN
165: T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
166: S = Z( NN-3 )*( Z( NN-5 ) / T )
167: IF( S.LE.T ) THEN
168: S = Z( NN-3 )*( Z( NN-5 ) /
169: $ ( T*( ONE+SQRT( ONE+S / T ) ) ) )
170: ELSE
171: S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
172: END IF
173: T = Z( NN-7 ) + ( S+Z( NN-5 ) )
174: Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
175: Z( NN-7 ) = T
176: END IF
177: Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
178: Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
179: N0 = N0 - 2
180: GO TO 10
181: *
182: 50 CONTINUE
183: IF( PP.EQ.2 )
184: $ PP = 0
185: *
186: * Reverse the qd-array, if warranted.
187: *
188: IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
189: IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
190: IPN4 = 4*( I0+N0 )
191: DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
192: TEMP = Z( J4-3 )
193: Z( J4-3 ) = Z( IPN4-J4-3 )
194: Z( IPN4-J4-3 ) = TEMP
195: TEMP = Z( J4-2 )
196: Z( J4-2 ) = Z( IPN4-J4-2 )
197: Z( IPN4-J4-2 ) = TEMP
198: TEMP = Z( J4-1 )
199: Z( J4-1 ) = Z( IPN4-J4-5 )
200: Z( IPN4-J4-5 ) = TEMP
201: TEMP = Z( J4 )
202: Z( J4 ) = Z( IPN4-J4-4 )
203: Z( IPN4-J4-4 ) = TEMP
204: 60 CONTINUE
205: IF( N0-I0.LE.4 ) THEN
206: Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
207: Z( 4*N0-PP ) = Z( 4*I0-PP )
208: END IF
209: DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
210: Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
211: $ Z( 4*I0+PP+3 ) )
212: Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
213: $ Z( 4*I0-PP+4 ) )
214: QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
215: DMIN = -ZERO
216: END IF
217: END IF
218: *
219: * Choose a shift.
220: *
221: CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
222: $ DN2, TAU, TTYPE, G )
223: *
224: * Call dqds until DMIN > 0.
225: *
226: 70 CONTINUE
227: *
228: CALL DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
229: $ DN1, DN2, IEEE )
230: *
231: NDIV = NDIV + ( N0-I0+2 )
232: ITER = ITER + 1
233: *
234: * Check status.
235: *
236: IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN
237: *
238: * Success.
239: *
240: GO TO 90
241: *
242: ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND.
243: $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
244: $ ABS( DN ).LT.TOL*SIGMA ) THEN
245: *
246: * Convergence hidden by negative DN.
247: *
248: Z( 4*( N0-1 )-PP+2 ) = ZERO
249: DMIN = ZERO
250: GO TO 90
251: ELSE IF( DMIN.LT.ZERO ) THEN
252: *
253: * TAU too big. Select new TAU and try again.
254: *
255: NFAIL = NFAIL + 1
256: IF( TTYPE.LT.-22 ) THEN
257: *
258: * Failed twice. Play it safe.
259: *
260: TAU = ZERO
261: ELSE IF( DMIN1.GT.ZERO ) THEN
262: *
263: * Late failure. Gives excellent shift.
264: *
265: TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
266: TTYPE = TTYPE - 11
267: ELSE
268: *
269: * Early failure. Divide by 4.
270: *
271: TAU = QURTR*TAU
272: TTYPE = TTYPE - 12
273: END IF
274: GO TO 70
275: ELSE IF( DISNAN( DMIN ) ) THEN
276: *
277: * NaN.
278: *
279: IF( TAU.EQ.ZERO ) THEN
280: GO TO 80
281: ELSE
282: TAU = ZERO
283: GO TO 70
284: END IF
285: ELSE
286: *
287: * Possible underflow. Play it safe.
288: *
289: GO TO 80
290: END IF
291: *
292: * Risk of underflow.
293: *
294: 80 CONTINUE
295: CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
296: NDIV = NDIV + ( N0-I0+2 )
297: ITER = ITER + 1
298: TAU = ZERO
299: *
300: 90 CONTINUE
301: IF( TAU.LT.SIGMA ) THEN
302: DESIG = DESIG + TAU
303: T = SIGMA + DESIG
304: DESIG = DESIG - ( T-SIGMA )
305: ELSE
306: T = SIGMA + TAU
307: DESIG = SIGMA - ( T-TAU ) + DESIG
308: END IF
309: SIGMA = T
310: *
311: RETURN
312: *
313: * End of DLASQ3
314: *
315: END