LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*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 lasq4
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
dlasq4
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