LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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
9 *> Download DLASQ4 + dependencies
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*N )
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 REAL
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 *> \date September 2012
139 *
140 *> \ingroup auxOTHERcomputational
141 *
142 *> \par Further Details:
143 * =====================
144 *>
145 *> \verbatim
146 *>
147 *> CNST1 = 9/16
148 *> \endverbatim
149 *>
150 * =====================================================================
151  SUBROUTINE dlasq4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
152  $ dn1, dn2, tau, ttype, g )
153 *
154 * -- LAPACK computational routine (version 3.4.2) --
155 * -- LAPACK is a software package provided by Univ. of Tennessee, --
156 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
157 * September 2012
158 *
159 * .. Scalar Arguments ..
160  INTEGER i0, n0, n0in, pp, ttype
161  DOUBLE PRECISION dmin, dmin1, dmin2, dn, dn1, dn2, g, tau
162 * ..
163 * .. Array Arguments ..
164  DOUBLE PRECISION z( * )
165 * ..
166 *
167 * =====================================================================
168 *
169 * .. Parameters ..
170  DOUBLE PRECISION cnst1, cnst2, cnst3
171  parameter( cnst1 = 0.5630d0, cnst2 = 1.010d0,
172  $ cnst3 = 1.050d0 )
173  DOUBLE PRECISION qurtr, third, half, zero, one, two, hundrd
174  parameter( qurtr = 0.250d0, third = 0.3330d0,
175  $ half = 0.50d0, zero = 0.0d0, one = 1.0d0,
176  $ two = 2.0d0, hundrd = 100.0d0 )
177 * ..
178 * .. Local Scalars ..
179  INTEGER i4, nn, np
180  DOUBLE PRECISION a2, b1, b2, gam, gap1, gap2, s
181 * ..
182 * .. Intrinsic Functions ..
183  INTRINSIC max, min, sqrt
184 * ..
185 * .. Executable Statements ..
186 *
187 * A negative DMIN forces the shift to take that absolute value
188 * TTYPE records the type of shift.
189 *
190  IF( dmin.LE.zero ) THEN
191  tau = -dmin
192  ttype = -1
193  return
194  END IF
195 *
196  nn = 4*n0 + pp
197  IF( n0in.EQ.n0 ) THEN
198 *
199 * No eigenvalues deflated.
200 *
201  IF( dmin.EQ.dn .OR. dmin.EQ.dn1 ) THEN
202 *
203  b1 = sqrt( z( nn-3 ) )*sqrt( z( nn-5 ) )
204  b2 = sqrt( z( nn-7 ) )*sqrt( z( nn-9 ) )
205  a2 = z( nn-7 ) + z( nn-5 )
206 *
207 * Cases 2 and 3.
208 *
209  IF( dmin.EQ.dn .AND. dmin1.EQ.dn1 ) THEN
210  gap2 = dmin2 - a2 - dmin2*qurtr
211  IF( gap2.GT.zero .AND. gap2.GT.b2 ) THEN
212  gap1 = a2 - dn - ( b2 / gap2 )*b2
213  ELSE
214  gap1 = a2 - dn - ( b1+b2 )
215  END IF
216  IF( gap1.GT.zero .AND. gap1.GT.b1 ) THEN
217  s = max( dn-( b1 / gap1 )*b1, half*dmin )
218  ttype = -2
219  ELSE
220  s = zero
221  IF( dn.GT.b1 )
222  $ s = dn - b1
223  IF( a2.GT.( b1+b2 ) )
224  $ s = min( s, a2-( b1+b2 ) )
225  s = max( s, third*dmin )
226  ttype = -3
227  END IF
228  ELSE
229 *
230 * Case 4.
231 *
232  ttype = -4
233  s = qurtr*dmin
234  IF( dmin.EQ.dn ) THEN
235  gam = dn
236  a2 = zero
237  IF( z( nn-5 ) .GT. z( nn-7 ) )
238  $ return
239  b2 = z( nn-5 ) / z( nn-7 )
240  np = nn - 9
241  ELSE
242  np = nn - 2*pp
243  b2 = z( np-2 )
244  gam = dn1
245  IF( z( np-4 ) .GT. z( np-2 ) )
246  $ return
247  a2 = z( np-4 ) / z( np-2 )
248  IF( z( nn-9 ) .GT. z( nn-11 ) )
249  $ return
250  b2 = z( nn-9 ) / z( nn-11 )
251  np = nn - 13
252  END IF
253 *
254 * Approximate contribution to norm squared from I < NN-1.
255 *
256  a2 = a2 + b2
257  DO 10 i4 = np, 4*i0 - 1 + pp, -4
258  IF( b2.EQ.zero )
259  $ go to 20
260  b1 = b2
261  IF( z( i4 ) .GT. z( i4-2 ) )
262  $ return
263  b2 = b2*( z( i4 ) / z( i4-2 ) )
264  a2 = a2 + b2
265  IF( hundrd*max( b2, b1 ).LT.a2 .OR. cnst1.LT.a2 )
266  $ go to 20
267  10 continue
268  20 continue
269  a2 = cnst3*a2
270 *
271 * Rayleigh quotient residual bound.
272 *
273  IF( a2.LT.cnst1 )
274  $ s = gam*( one-sqrt( a2 ) ) / ( one+a2 )
275  END IF
276  ELSE IF( dmin.EQ.dn2 ) THEN
277 *
278 * Case 5.
279 *
280  ttype = -5
281  s = qurtr*dmin
282 *
283 * Compute contribution to norm squared from I > NN-2.
284 *
285  np = nn - 2*pp
286  b1 = z( np-2 )
287  b2 = z( np-6 )
288  gam = dn2
289  IF( z( np-8 ).GT.b2 .OR. z( np-4 ).GT.b1 )
290  $ return
291  a2 = ( z( np-8 ) / b2 )*( one+z( np-4 ) / b1 )
292 *
293 * Approximate contribution to norm squared from I < NN-2.
294 *
295  IF( n0-i0.GT.2 ) THEN
296  b2 = z( nn-13 ) / z( nn-15 )
297  a2 = a2 + b2
298  DO 30 i4 = nn - 17, 4*i0 - 1 + pp, -4
299  IF( b2.EQ.zero )
300  $ go to 40
301  b1 = b2
302  IF( z( i4 ) .GT. z( i4-2 ) )
303  $ return
304  b2 = b2*( z( i4 ) / z( i4-2 ) )
305  a2 = a2 + b2
306  IF( hundrd*max( b2, b1 ).LT.a2 .OR. cnst1.LT.a2 )
307  $ go to 40
308  30 continue
309  40 continue
310  a2 = cnst3*a2
311  END IF
312 *
313  IF( a2.LT.cnst1 )
314  $ s = gam*( one-sqrt( a2 ) ) / ( one+a2 )
315  ELSE
316 *
317 * Case 6, no information to guide us.
318 *
319  IF( ttype.EQ.-6 ) THEN
320  g = g + third*( one-g )
321  ELSE IF( ttype.EQ.-18 ) THEN
322  g = qurtr*third
323  ELSE
324  g = qurtr
325  END IF
326  s = g*dmin
327  ttype = -6
328  END IF
329 *
330  ELSE IF( n0in.EQ.( n0+1 ) ) THEN
331 *
332 * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
333 *
334  IF( dmin1.EQ.dn1 .AND. dmin2.EQ.dn2 ) THEN
335 *
336 * Cases 7 and 8.
337 *
338  ttype = -7
339  s = third*dmin1
340  IF( z( nn-5 ).GT.z( nn-7 ) )
341  $ return
342  b1 = z( nn-5 ) / z( nn-7 )
343  b2 = b1
344  IF( b2.EQ.zero )
345  $ go to 60
346  DO 50 i4 = 4*n0 - 9 + pp, 4*i0 - 1 + pp, -4
347  a2 = b1
348  IF( z( i4 ).GT.z( i4-2 ) )
349  $ return
350  b1 = b1*( z( i4 ) / z( i4-2 ) )
351  b2 = b2 + b1
352  IF( hundrd*max( b1, a2 ).LT.b2 )
353  $ go to 60
354  50 continue
355  60 continue
356  b2 = sqrt( cnst3*b2 )
357  a2 = dmin1 / ( one+b2**2 )
358  gap2 = half*dmin2 - a2
359  IF( gap2.GT.zero .AND. gap2.GT.b2*a2 ) THEN
360  s = max( s, a2*( one-cnst2*a2*( b2 / gap2 )*b2 ) )
361  ELSE
362  s = max( s, a2*( one-cnst2*b2 ) )
363  ttype = -8
364  END IF
365  ELSE
366 *
367 * Case 9.
368 *
369  s = qurtr*dmin1
370  IF( dmin1.EQ.dn1 )
371  $ s = half*dmin1
372  ttype = -9
373  END IF
374 *
375  ELSE IF( n0in.EQ.( n0+2 ) ) THEN
376 *
377 * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
378 *
379 * Cases 10 and 11.
380 *
381  IF( dmin2.EQ.dn2 .AND. two*z( nn-5 ).LT.z( nn-7 ) ) THEN
382  ttype = -10
383  s = third*dmin2
384  IF( z( nn-5 ).GT.z( nn-7 ) )
385  $ return
386  b1 = z( nn-5 ) / z( nn-7 )
387  b2 = b1
388  IF( b2.EQ.zero )
389  $ go to 80
390  DO 70 i4 = 4*n0 - 9 + pp, 4*i0 - 1 + pp, -4
391  IF( z( i4 ).GT.z( i4-2 ) )
392  $ return
393  b1 = b1*( z( i4 ) / z( i4-2 ) )
394  b2 = b2 + b1
395  IF( hundrd*b1.LT.b2 )
396  $ go to 80
397  70 continue
398  80 continue
399  b2 = sqrt( cnst3*b2 )
400  a2 = dmin2 / ( one+b2**2 )
401  gap2 = z( nn-7 ) + z( nn-9 ) -
402  $ sqrt( z( nn-11 ) )*sqrt( z( nn-9 ) ) - a2
403  IF( gap2.GT.zero .AND. gap2.GT.b2*a2 ) THEN
404  s = max( s, a2*( one-cnst2*a2*( b2 / gap2 )*b2 ) )
405  ELSE
406  s = max( s, a2*( one-cnst2*b2 ) )
407  END IF
408  ELSE
409  s = qurtr*dmin2
410  ttype = -11
411  END IF
412  ELSE IF( n0in.GT.( n0+2 ) ) THEN
413 *
414 * Case 12, more than two eigenvalues deflated. No information.
415 *
416  s = zero
417  ttype = -12
418  END IF
419 *
420  tau = s
421  return
422 *
423 * End of DLASQ4
424 *
425  END