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*> \date June 2016
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.7.1) --
155*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
156*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
157*     June 2016
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                  GAM = DN1
244                  IF( Z( NP-4 ) .GT. Z( NP-2 ) )
245     $               RETURN
246                  A2 = Z( NP-4 ) / Z( NP-2 )
247                  IF( Z( NN-9 ) .GT. Z( NN-11 ) )
248     $               RETURN
249                  B2 = Z( NN-9 ) / Z( NN-11 )
250                  NP = NN - 13
251               END IF
252*
253*              Approximate contribution to norm squared from I < NN-1.
254*
255               A2 = A2 + B2
256               DO 10 I4 = NP, 4*I0 - 1 + PP, -4
257                  IF( B2.EQ.ZERO )
258     $               GO TO 20
259                  B1 = B2
260                  IF( Z( I4 ) .GT. Z( I4-2 ) )
261     $               RETURN
262                  B2 = B2*( Z( I4 ) / Z( I4-2 ) )
263                  A2 = A2 + B2
264                  IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
265     $               GO TO 20
266   10          CONTINUE
267   20          CONTINUE
268               A2 = CNST3*A2
269*
270*              Rayleigh quotient residual bound.
271*
272               IF( A2.LT.CNST1 )
273     $            S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
274            END IF
275         ELSE IF( DMIN.EQ.DN2 ) THEN
276*
277*           Case 5.
278*
279            TTYPE = -5
280            S = QURTR*DMIN
281*
282*           Compute contribution to norm squared from I > NN-2.
283*
284            NP = NN - 2*PP
285            B1 = Z( NP-2 )
286            B2 = Z( NP-6 )
287            GAM = DN2
288            IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 )
289     $         RETURN
290            A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 )
291*
292*           Approximate contribution to norm squared from I < NN-2.
293*
294            IF( N0-I0.GT.2 ) THEN
295               B2 = Z( NN-13 ) / Z( NN-15 )
296               A2 = A2 + B2
297               DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4
298                  IF( B2.EQ.ZERO )
299     $               GO TO 40
300                  B1 = B2
301                  IF( Z( I4 ) .GT. Z( I4-2 ) )
302     $               RETURN
303                  B2 = B2*( Z( I4 ) / Z( I4-2 ) )
304                  A2 = A2 + B2
305                  IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
306     $               GO TO 40
307   30          CONTINUE
308   40          CONTINUE
309               A2 = CNST3*A2
310            END IF
311*
312            IF( A2.LT.CNST1 )
313     $         S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
314         ELSE
315*
316*           Case 6, no information to guide us.
317*
318            IF( TTYPE.EQ.-6 ) THEN
319               G = G + THIRD*( ONE-G )
320            ELSE IF( TTYPE.EQ.-18 ) THEN
321               G = QURTR*THIRD
322            ELSE
323               G = QURTR
324            END IF
325            S = G*DMIN
326            TTYPE = -6
327         END IF
328*
329      ELSE IF( N0IN.EQ.( N0+1 ) ) THEN
330*
331*        One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
332*
333         IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN
334*
335*           Cases 7 and 8.
336*
337            TTYPE = -7
338            S = THIRD*DMIN1
339            IF( Z( NN-5 ).GT.Z( NN-7 ) )
340     $         RETURN
341            B1 = Z( NN-5 ) / Z( NN-7 )
342            B2 = B1
343            IF( B2.EQ.ZERO )
344     $         GO TO 60
345            DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
346               A2 = B1
347               IF( Z( I4 ).GT.Z( I4-2 ) )
348     $            RETURN
349               B1 = B1*( Z( I4 ) / Z( I4-2 ) )
350               B2 = B2 + B1
351               IF( HUNDRD*MAX( B1, A2 ).LT.B2 )
352     $            GO TO 60
353   50       CONTINUE
354   60       CONTINUE
355            B2 = SQRT( CNST3*B2 )
356            A2 = DMIN1 / ( ONE+B2**2 )
357            GAP2 = HALF*DMIN2 - A2
358            IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
359               S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
360            ELSE
361               S = MAX( S, A2*( ONE-CNST2*B2 ) )
362               TTYPE = -8
363            END IF
364         ELSE
365*
366*           Case 9.
367*
368            S = QURTR*DMIN1
369            IF( DMIN1.EQ.DN1 )
370     $         S = HALF*DMIN1
371            TTYPE = -9
372         END IF
373*
374      ELSE IF( N0IN.EQ.( N0+2 ) ) THEN
375*
376*        Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
377*
378*        Cases 10 and 11.
379*
380         IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN
381            TTYPE = -10
382            S = THIRD*DMIN2
383            IF( Z( NN-5 ).GT.Z( NN-7 ) )
384     $         RETURN
385            B1 = Z( NN-5 ) / Z( NN-7 )
386            B2 = B1
387            IF( B2.EQ.ZERO )
388     $         GO TO 80
389            DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
390               IF( Z( I4 ).GT.Z( I4-2 ) )
391     $            RETURN
392               B1 = B1*( Z( I4 ) / Z( I4-2 ) )
393               B2 = B2 + B1
394               IF( HUNDRD*B1.LT.B2 )
395     $            GO TO 80
396   70       CONTINUE
397   80       CONTINUE
398            B2 = SQRT( CNST3*B2 )
399            A2 = DMIN2 / ( ONE+B2**2 )
400            GAP2 = Z( NN-7 ) + Z( NN-9 ) -
401     $             SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2
402            IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
403               S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
404            ELSE
405               S = MAX( S, A2*( ONE-CNST2*B2 ) )
406            END IF
407         ELSE
408            S = QURTR*DMIN2
409            TTYPE = -11
410         END IF
411      ELSE IF( N0IN.GT.( N0+2 ) ) THEN
412*
413*        Case 12, more than two eigenvalues deflated. No information.
414*
415         S = ZERO
416         TTYPE = -12
417      END IF
418*
419      TAU = S
420      RETURN
421*
422*     End of DLASQ4
423*
424      END
425