KTH framework for Nek5000 toolboxes; testing version  0.0.1
dlasq4.f
Go to the documentation of this file.
1  SUBROUTINE dlasq4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
2  $ DN1, DN2, TAU, TTYPE )
3 *
4 * -- LAPACK auxiliary routine (version 3.0) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6 * Courant Institute, Argonne National Lab, and Rice University
7 * October 31, 1999
8 *
9 * .. Scalar Arguments ..
10  INTEGER I0, N0, N0IN, PP, TTYPE
11  DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DN1, DN2, TAU
12 * ..
13 * .. Array Arguments ..
14  DOUBLE PRECISION Z( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DLASQ4 computes an approximation TAU to the smallest eigenvalue
21 * using values of d from the previous transform.
22 *
23 * I0 (input) INTEGER
24 * First index.
25 *
26 * N0 (input) INTEGER
27 * Last index.
28 *
29 * Z (input) DOUBLE PRECISION array, dimension ( 4*N )
30 * Z holds the qd array.
31 *
32 * PP (input) INTEGER
33 * PP=0 for ping, PP=1 for pong.
34 *
35 * NOIN (input) INTEGER
36 * The value of N0 at start of EIGTEST.
37 *
38 * DMIN (input) DOUBLE PRECISION
39 * Minimum value of d.
40 *
41 * DMIN1 (input) DOUBLE PRECISION
42 * Minimum value of d, excluding D( N0 ).
43 *
44 * DMIN2 (input) DOUBLE PRECISION
45 * Minimum value of d, excluding D( N0 ) and D( N0-1 ).
46 *
47 * DN (input) DOUBLE PRECISION
48 * d(N)
49 *
50 * DN1 (input) DOUBLE PRECISION
51 * d(N-1)
52 *
53 * DN2 (input) DOUBLE PRECISION
54 * d(N-2)
55 *
56 * TAU (output) DOUBLE PRECISION
57 * This is the shift.
58 *
59 * TTYPE (output) INTEGER
60 * Shift type.
61 *
62 * Further Details
63 * ===============
64 * CNST1 = 9/16
65 *
66 * =====================================================================
67 *
68 * .. Parameters ..
69  DOUBLE PRECISION CNST1, CNST2, CNST3
70  parameter( cnst1 = 0.5630d0, cnst2 = 1.010d0,
71  $ cnst3 = 1.050d0 )
72  DOUBLE PRECISION QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
73  parameter( qurtr = 0.250d0, third = 0.3330d0,
74  $ half = 0.50d0, zero = 0.0d0, one = 1.0d0,
75  $ two = 2.0d0, hundrd = 100.0d0 )
76 * ..
77 * .. Local Scalars ..
78  INTEGER I4, NN, NP
79  DOUBLE PRECISION A2, B1, B2, G, GAM, GAP1, GAP2, S
80 * ..
81 * .. Intrinsic Functions ..
82  INTRINSIC max, min, sqrt
83 * ..
84 * .. Save statement ..
85  SAVE g
86 * ..
87 * .. Data statement ..
88  DATA g / zero /
89 * ..
90 * .. Executable Statements ..
91 *
92 * A negative DMIN forces the shift to take that absolute value
93 * TTYPE records the type of shift.
94 *
95  IF( dmin.LE.zero ) THEN
96  tau = -dmin
97  ttype = -1
98  RETURN
99  END IF
100 *
101  nn = 4*n0 + pp
102  IF( n0in.EQ.n0 ) THEN
103 *
104 * No eigenvalues deflated.
105 *
106  IF( dmin.EQ.dn .OR. dmin.EQ.dn1 ) THEN
107 *
108  b1 = sqrt( z( nn-3 ) )*sqrt( z( nn-5 ) )
109  b2 = sqrt( z( nn-7 ) )*sqrt( z( nn-9 ) )
110  a2 = z( nn-7 ) + z( nn-5 )
111 *
112 * Cases 2 and 3.
113 *
114  IF( dmin.EQ.dn .AND. dmin1.EQ.dn1 ) THEN
115  gap2 = dmin2 - a2 - dmin2*qurtr
116  IF( gap2.GT.zero .AND. gap2.GT.b2 ) THEN
117  gap1 = a2 - dn - ( b2 / gap2 )*b2
118  ELSE
119  gap1 = a2 - dn - ( b1+b2 )
120  END IF
121  IF( gap1.GT.zero .AND. gap1.GT.b1 ) THEN
122  s = max( dn-( b1 / gap1 )*b1, half*dmin )
123  ttype = -2
124  ELSE
125  s = zero
126  IF( dn.GT.b1 )
127  $ s = dn - b1
128  IF( a2.GT.( b1+b2 ) )
129  $ s = min( s, a2-( b1+b2 ) )
130  s = max( s, third*dmin )
131  ttype = -3
132  END IF
133  ELSE
134 *
135 * Case 4.
136 *
137  ttype = -4
138  s = qurtr*dmin
139  IF( dmin.EQ.dn ) THEN
140  gam = dn
141  a2 = zero
142  IF( z( nn-5 ) .GT. z( nn-7 ) )
143  $ RETURN
144  b2 = z( nn-5 ) / z( nn-7 )
145  np = nn - 9
146  ELSE
147  np = nn - 2*pp
148  b2 = z( np-2 )
149  gam = dn1
150  IF( z( np-4 ) .GT. z( np-2 ) )
151  $ RETURN
152  a2 = z( np-4 ) / z( np-2 )
153  IF( z( nn-9 ) .GT. z( nn-11 ) )
154  $ RETURN
155  b2 = z( nn-9 ) / z( nn-11 )
156  np = nn - 13
157  END IF
158 *
159 * Approximate contribution to norm squared from I < NN-1.
160 *
161  a2 = a2 + b2
162  DO 10 i4 = np, 4*i0 - 1 + pp, -4
163  IF( b2.EQ.zero )
164  $ GO TO 20
165  b1 = b2
166  IF( z( i4 ) .GT. z( i4-2 ) )
167  $ RETURN
168  b2 = b2*( z( i4 ) / z( i4-2 ) )
169  a2 = a2 + b2
170  IF( hundrd*max( b2, b1 ).LT.a2 .OR. cnst1.LT.a2 )
171  $ GO TO 20
172  10 CONTINUE
173  20 CONTINUE
174  a2 = cnst3*a2
175 *
176 * Rayleigh quotient residual bound.
177 *
178  IF( a2.LT.cnst1 )
179  $ s = gam*( one-sqrt( a2 ) ) / ( one+a2 )
180  END IF
181  ELSE IF( dmin.EQ.dn2 ) THEN
182 *
183 * Case 5.
184 *
185  ttype = -5
186  s = qurtr*dmin
187 *
188 * Compute contribution to norm squared from I > NN-2.
189 *
190  np = nn - 2*pp
191  b1 = z( np-2 )
192  b2 = z( np-6 )
193  gam = dn2
194  IF( z( np-8 ).GT.b2 .OR. z( np-4 ).GT.b1 )
195  $ RETURN
196  a2 = ( z( np-8 ) / b2 )*( one+z( np-4 ) / b1 )
197 *
198 * Approximate contribution to norm squared from I < NN-2.
199 *
200  IF( n0-i0.GT.2 ) THEN
201  b2 = z( nn-13 ) / z( nn-15 )
202  a2 = a2 + b2
203  DO 30 i4 = nn - 17, 4*i0 - 1 + pp, -4
204  IF( b2.EQ.zero )
205  $ GO TO 40
206  b1 = b2
207  IF( z( i4 ) .GT. z( i4-2 ) )
208  $ RETURN
209  b2 = b2*( z( i4 ) / z( i4-2 ) )
210  a2 = a2 + b2
211  IF( hundrd*max( b2, b1 ).LT.a2 .OR. cnst1.LT.a2 )
212  $ GO TO 40
213  30 CONTINUE
214  40 CONTINUE
215  a2 = cnst3*a2
216  END IF
217 *
218  IF( a2.LT.cnst1 )
219  $ s = gam*( one-sqrt( a2 ) ) / ( one+a2 )
220  ELSE
221 *
222 * Case 6, no information to guide us.
223 *
224  IF( ttype.EQ.-6 ) THEN
225  g = g + third*( one-g )
226  ELSE IF( ttype.EQ.-18 ) THEN
227  g = qurtr*third
228  ELSE
229  g = qurtr
230  END IF
231  s = g*dmin
232  ttype = -6
233  END IF
234 *
235  ELSE IF( n0in.EQ.( n0+1 ) ) THEN
236 *
237 * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
238 *
239  IF( dmin1.EQ.dn1 .AND. dmin2.EQ.dn2 ) THEN
240 *
241 * Cases 7 and 8.
242 *
243  ttype = -7
244  s = third*dmin1
245  IF( z( nn-5 ).GT.z( nn-7 ) )
246  $ RETURN
247  b1 = z( nn-5 ) / z( nn-7 )
248  b2 = b1
249  IF( b2.EQ.zero )
250  $ GO TO 60
251  DO 50 i4 = 4*n0 - 9 + pp, 4*i0 - 1 + pp, -4
252  a2 = b1
253  IF( z( i4 ).GT.z( i4-2 ) )
254  $ RETURN
255  b1 = b1*( z( i4 ) / z( i4-2 ) )
256  b2 = b2 + b1
257  IF( hundrd*max( b1, a2 ).LT.b2 )
258  $ GO TO 60
259  50 CONTINUE
260  60 CONTINUE
261  b2 = sqrt( cnst3*b2 )
262  a2 = dmin1 / ( one+b2**2 )
263  gap2 = half*dmin2 - a2
264  IF( gap2.GT.zero .AND. gap2.GT.b2*a2 ) THEN
265  s = max( s, a2*( one-cnst2*a2*( b2 / gap2 )*b2 ) )
266  ELSE
267  s = max( s, a2*( one-cnst2*b2 ) )
268  ttype = -8
269  END IF
270  ELSE
271 *
272 * Case 9.
273 *
274  s = qurtr*dmin1
275  IF( dmin1.EQ.dn1 )
276  $ s = half*dmin1
277  ttype = -9
278  END IF
279 *
280  ELSE IF( n0in.EQ.( n0+2 ) ) THEN
281 *
282 * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
283 *
284 * Cases 10 and 11.
285 *
286  IF( dmin2.EQ.dn2 .AND. two*z( nn-5 ).LT.z( nn-7 ) ) THEN
287  ttype = -10
288  s = third*dmin2
289  IF( z( nn-5 ).GT.z( nn-7 ) )
290  $ RETURN
291  b1 = z( nn-5 ) / z( nn-7 )
292  b2 = b1
293  IF( b2.EQ.zero )
294  $ GO TO 80
295  DO 70 i4 = 4*n0 - 9 + pp, 4*i0 - 1 + pp, -4
296  IF( z( i4 ).GT.z( i4-2 ) )
297  $ RETURN
298  b1 = b1*( z( i4 ) / z( i4-2 ) )
299  b2 = b2 + b1
300  IF( hundrd*b1.LT.b2 )
301  $ GO TO 80
302  70 CONTINUE
303  80 CONTINUE
304  b2 = sqrt( cnst3*b2 )
305  a2 = dmin2 / ( one+b2**2 )
306  gap2 = z( nn-7 ) + z( nn-9 ) -
307  $ sqrt( z( nn-11 ) )*sqrt( z( nn-9 ) ) - a2
308  IF( gap2.GT.zero .AND. gap2.GT.b2*a2 ) THEN
309  s = max( s, a2*( one-cnst2*a2*( b2 / gap2 )*b2 ) )
310  ELSE
311  s = max( s, a2*( one-cnst2*b2 ) )
312  END IF
313  ELSE
314  s = qurtr*dmin2
315  ttype = -11
316  END IF
317  ELSE IF( n0in.GT.( n0+2 ) ) THEN
318 *
319 * Case 12, more than two eigenvalues deflated. No information.
320 *
321  s = zero
322  ttype = -12
323  END IF
324 *
325  tau = s
326  RETURN
327 *
328 * End of DLASQ4
329 *
330  END
dlasq4
subroutine dlasq4(I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, DN2, TAU, TTYPE)
Definition: dlasq4.f:3