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 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
422