LAPACK  3.10.1 LAPACK: Linear Algebra PACKage
dlasq4.f
Go to the documentation of this file.
1 *> \brief \b DLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq4.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq4.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq4.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
22 * DN1, DN2, TAU, TTYPE, G )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER I0, N0, N0IN, PP, TTYPE
26 * DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
27 * ..
28 * .. Array Arguments ..
29 * DOUBLE PRECISION Z( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> DLASQ4 computes an approximation TAU to the smallest eigenvalue
39 *> using values of d from the previous transform.
40 *> \endverbatim
41 *
42 * Arguments:
43 * ==========
44 *
45 *> \param[in] I0
46 *> \verbatim
47 *> I0 is INTEGER
48 *> First index.
49 *> \endverbatim
50 *>
51 *> \param[in] N0
52 *> \verbatim
53 *> N0 is INTEGER
54 *> Last index.
55 *> \endverbatim
56 *>
57 *> \param[in] Z
58 *> \verbatim
59 *> Z is DOUBLE PRECISION array, dimension ( 4*N0 )
60 *> Z holds the qd array.
61 *> \endverbatim
62 *>
63 *> \param[in] PP
64 *> \verbatim
65 *> PP is INTEGER
66 *> PP=0 for ping, PP=1 for pong.
67 *> \endverbatim
68 *>
69 *> \param[in] N0IN
70 *> \verbatim
71 *> N0IN is INTEGER
72 *> The value of N0 at start of EIGTEST.
73 *> \endverbatim
74 *>
75 *> \param[in] DMIN
76 *> \verbatim
77 *> DMIN is DOUBLE PRECISION
78 *> Minimum value of d.
79 *> \endverbatim
80 *>
81 *> \param[in] DMIN1
82 *> \verbatim
83 *> DMIN1 is DOUBLE PRECISION
84 *> Minimum value of d, excluding D( N0 ).
85 *> \endverbatim
86 *>
87 *> \param[in] DMIN2
88 *> \verbatim
89 *> DMIN2 is DOUBLE PRECISION
90 *> Minimum value of d, excluding D( N0 ) and D( N0-1 ).
91 *> \endverbatim
92 *>
93 *> \param[in] DN
94 *> \verbatim
95 *> DN is DOUBLE PRECISION
96 *> d(N)
97 *> \endverbatim
98 *>
99 *> \param[in] DN1
100 *> \verbatim
101 *> DN1 is DOUBLE PRECISION
102 *> d(N-1)
103 *> \endverbatim
104 *>
105 *> \param[in] DN2
106 *> \verbatim
107 *> DN2 is DOUBLE PRECISION
108 *> d(N-2)
109 *> \endverbatim
110 *>
111 *> \param[out] TAU
112 *> \verbatim
113 *> TAU is DOUBLE PRECISION
114 *> This is the shift.
115 *> \endverbatim
116 *>
117 *> \param[out] TTYPE
118 *> \verbatim
119 *> TTYPE is INTEGER
120 *> Shift type.
121 *> \endverbatim
122 *>
123 *> \param[in,out] G
124 *> \verbatim
125 *> G is DOUBLE PRECISION
126 *> G is passed as an argument in order to save its value between
127 *> calls to DLASQ4.
128 *> \endverbatim
129 *
130 * Authors:
131 * ========
132 *
133 *> \author Univ. of Tennessee
134 *> \author Univ. of California Berkeley
135 *> \author Univ. of Colorado Denver
136 *> \author NAG Ltd.
137 *
138 *> \ingroup auxOTHERcomputational
139 *
140 *> \par Further Details:
141 * =====================
142 *>
143 *> \verbatim
144 *>
145 *> CNST1 = 9/16
146 *> \endverbatim
147 *>
148 * =====================================================================
149  SUBROUTINE dlasq4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
150  \$ DN1, DN2, TAU, TTYPE, G )
151 *
152 * -- LAPACK computational routine --
153 * -- LAPACK is a software package provided by Univ. of Tennessee, --
154 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155 *
156 * .. Scalar Arguments ..
157  INTEGER I0, N0, N0IN, PP, TTYPE
158  DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
159 * ..
160 * .. Array Arguments ..
161  DOUBLE PRECISION Z( * )
162 * ..
163 *
164 * =====================================================================
165 *
166 * .. Parameters ..
167  DOUBLE PRECISION CNST1, CNST2, CNST3
168  parameter( cnst1 = 0.5630d0, cnst2 = 1.010d0,
169  \$ cnst3 = 1.050d0 )
170  DOUBLE PRECISION QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
171  parameter( qurtr = 0.250d0, third = 0.3330d0,
172  \$ half = 0.50d0, zero = 0.0d0, one = 1.0d0,
173  \$ two = 2.0d0, hundrd = 100.0d0 )
174 * ..
175 * .. Local Scalars ..
176  INTEGER I4, NN, NP
177  DOUBLE PRECISION A2, B1, B2, GAM, GAP1, GAP2, S
178 * ..
179 * .. Intrinsic Functions ..
180  INTRINSIC max, min, sqrt
181 * ..
182 * .. Executable Statements ..
183 *
184 * A negative DMIN forces the shift to take that absolute value
185 * TTYPE records the type of shift.
186 *
187  IF( dmin.LE.zero ) THEN
188  tau = -dmin
189  ttype = -1
190  RETURN
191  END IF
192 *
193  nn = 4*n0 + pp
194  IF( n0in.EQ.n0 ) THEN
195 *
196 * No eigenvalues deflated.
197 *
198  IF( dmin.EQ.dn .OR. dmin.EQ.dn1 ) THEN
199 *
200  b1 = sqrt( z( nn-3 ) )*sqrt( z( nn-5 ) )
201  b2 = sqrt( z( nn-7 ) )*sqrt( z( nn-9 ) )
202  a2 = z( nn-7 ) + z( nn-5 )
203 *
204 * Cases 2 and 3.
205 *
206  IF( dmin.EQ.dn .AND. dmin1.EQ.dn1 ) THEN
207  gap2 = dmin2 - a2 - dmin2*qurtr
208  IF( gap2.GT.zero .AND. gap2.GT.b2 ) THEN
209  gap1 = a2 - dn - ( b2 / gap2 )*b2
210  ELSE
211  gap1 = a2 - dn - ( b1+b2 )
212  END IF
213  IF( gap1.GT.zero .AND. gap1.GT.b1 ) THEN
214  s = max( dn-( b1 / gap1 )*b1, half*dmin )
215  ttype = -2
216  ELSE
217  s = zero
218  IF( dn.GT.b1 )
219  \$ s = dn - b1
220  IF( a2.GT.( b1+b2 ) )
221  \$ s = min( s, a2-( b1+b2 ) )
222  s = max( s, third*dmin )
223  ttype = -3
224  END IF
225  ELSE
226 *
227 * Case 4.
228 *
229  ttype = -4
230  s = qurtr*dmin
231  IF( dmin.EQ.dn ) THEN
232  gam = dn
233  a2 = zero
234  IF( z( nn-5 ) .GT. z( nn-7 ) )
235  \$ RETURN
236  b2 = z( nn-5 ) / z( nn-7 )
237  np = nn - 9
238  ELSE
239  np = nn - 2*pp
240  gam = dn1
241  IF( z( np-4 ) .GT. z( np-2 ) )
242  \$ RETURN
243  a2 = z( np-4 ) / z( np-2 )
244  IF( z( nn-9 ) .GT. z( nn-11 ) )
245  \$ RETURN
246  b2 = z( nn-9 ) / z( nn-11 )
247  np = nn - 13
248  END IF
249 *
250 * Approximate contribution to norm squared from I < NN-1.
251 *
252  a2 = a2 + b2
253  DO 10 i4 = np, 4*i0 - 1 + pp, -4
254  IF( b2.EQ.zero )
255  \$ GO TO 20
256  b1 = b2
257  IF( z( i4 ) .GT. z( i4-2 ) )
258  \$ RETURN
259  b2 = b2*( z( i4 ) / z( i4-2 ) )
260  a2 = a2 + b2
261  IF( hundrd*max( b2, b1 ).LT.a2 .OR. cnst1.LT.a2 )
262  \$ GO TO 20
263  10 CONTINUE
264  20 CONTINUE
265  a2 = cnst3*a2
266 *
267 * Rayleigh quotient residual bound.
268 *
269  IF( a2.LT.cnst1 )
270  \$ s = gam*( one-sqrt( a2 ) ) / ( one+a2 )
271  END IF
272  ELSE IF( dmin.EQ.dn2 ) THEN
273 *
274 * Case 5.
275 *
276  ttype = -5
277  s = qurtr*dmin
278 *
279 * Compute contribution to norm squared from I > NN-2.
280 *
281  np = nn - 2*pp
282  b1 = z( np-2 )
283  b2 = z( np-6 )
284  gam = dn2
285  IF( z( np-8 ).GT.b2 .OR. z( np-4 ).GT.b1 )
286  \$ RETURN
287  a2 = ( z( np-8 ) / b2 )*( one+z( np-4 ) / b1 )
288 *
289 * Approximate contribution to norm squared from I < NN-2.
290 *
291  IF( n0-i0.GT.2 ) THEN
292  b2 = z( nn-13 ) / z( nn-15 )
293  a2 = a2 + b2
294  DO 30 i4 = nn - 17, 4*i0 - 1 + pp, -4
295  IF( b2.EQ.zero )
296  \$ GO TO 40
297  b1 = b2
298  IF( z( i4 ) .GT. z( i4-2 ) )
299  \$ RETURN
300  b2 = b2*( z( i4 ) / z( i4-2 ) )
301  a2 = a2 + b2
302  IF( hundrd*max( b2, b1 ).LT.a2 .OR. cnst1.LT.a2 )
303  \$ GO TO 40
304  30 CONTINUE
305  40 CONTINUE
306  a2 = cnst3*a2
307  END IF
308 *
309  IF( a2.LT.cnst1 )
310  \$ s = gam*( one-sqrt( a2 ) ) / ( one+a2 )
311  ELSE
312 *
313 * Case 6, no information to guide us.
314 *
315  IF( ttype.EQ.-6 ) THEN
316  g = g + third*( one-g )
317  ELSE IF( ttype.EQ.-18 ) THEN
318  g = qurtr*third
319  ELSE
320  g = qurtr
321  END IF
322  s = g*dmin
323  ttype = -6
324  END IF
325 *
326  ELSE IF( n0in.EQ.( n0+1 ) ) THEN
327 *
328 * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
329 *
330  IF( dmin1.EQ.dn1 .AND. dmin2.EQ.dn2 ) THEN
331 *
332 * Cases 7 and 8.
333 *
334  ttype = -7
335  s = third*dmin1
336  IF( z( nn-5 ).GT.z( nn-7 ) )
337  \$ RETURN
338  b1 = z( nn-5 ) / z( nn-7 )
339  b2 = b1
340  IF( b2.EQ.zero )
341  \$ GO TO 60
342  DO 50 i4 = 4*n0 - 9 + pp, 4*i0 - 1 + pp, -4
343  a2 = b1
344  IF( z( i4 ).GT.z( i4-2 ) )
345  \$ RETURN
346  b1 = b1*( z( i4 ) / z( i4-2 ) )
347  b2 = b2 + b1
348  IF( hundrd*max( b1, a2 ).LT.b2 )
349  \$ GO TO 60
350  50 CONTINUE
351  60 CONTINUE
352  b2 = sqrt( cnst3*b2 )
353  a2 = dmin1 / ( one+b2**2 )
354  gap2 = half*dmin2 - a2
355  IF( gap2.GT.zero .AND. gap2.GT.b2*a2 ) THEN
356  s = max( s, a2*( one-cnst2*a2*( b2 / gap2 )*b2 ) )
357  ELSE
358  s = max( s, a2*( one-cnst2*b2 ) )
359  ttype = -8
360  END IF
361  ELSE
362 *
363 * Case 9.
364 *
365  s = qurtr*dmin1
366  IF( dmin1.EQ.dn1 )
367  \$ s = half*dmin1
368  ttype = -9
369  END IF
370 *
371  ELSE IF( n0in.EQ.( n0+2 ) ) THEN
372 *
373 * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
374 *
375 * Cases 10 and 11.
376 *
377  IF( dmin2.EQ.dn2 .AND. two*z( nn-5 ).LT.z( nn-7 ) ) THEN
378  ttype = -10
379  s = third*dmin2
380  IF( z( nn-5 ).GT.z( nn-7 ) )
381  \$ RETURN
382  b1 = z( nn-5 ) / z( nn-7 )
383  b2 = b1
384  IF( b2.EQ.zero )
385  \$ GO TO 80
386  DO 70 i4 = 4*n0 - 9 + pp, 4*i0 - 1 + pp, -4
387  IF( z( i4 ).GT.z( i4-2 ) )
388  \$ RETURN
389  b1 = b1*( z( i4 ) / z( i4-2 ) )
390  b2 = b2 + b1
391  IF( hundrd*b1.LT.b2 )
392  \$ GO TO 80
393  70 CONTINUE
394  80 CONTINUE
395  b2 = sqrt( cnst3*b2 )
396  a2 = dmin2 / ( one+b2**2 )
397  gap2 = z( nn-7 ) + z( nn-9 ) -
398  \$ sqrt( z( nn-11 ) )*sqrt( z( nn-9 ) ) - a2
399  IF( gap2.GT.zero .AND. gap2.GT.b2*a2 ) THEN
400  s = max( s, a2*( one-cnst2*a2*( b2 / gap2 )*b2 ) )
401  ELSE
402  s = max( s, a2*( one-cnst2*b2 ) )
403  END IF
404  ELSE
405  s = qurtr*dmin2
406  ttype = -11
407  END IF
408  ELSE IF( n0in.GT.( n0+2 ) ) THEN
409 *
410 * Case 12, more than two eigenvalues deflated. No information.
411 *
412  s = zero
413  ttype = -12
414  END IF
415 *
416  tau = s
417  RETURN
418 *
419 * End of DLASQ4
420 *
421  END
subroutine dlasq4(I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, DN2, TAU, TTYPE, G)
DLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous trans...
Definition: dlasq4.f:151