1*> \brief \b DLASQ3 checks for deflation, computes a shift and calls dqds. 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 DLASQ3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
22*                          ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
23*                          DN2, G, TAU )
24*
25*       .. Scalar Arguments ..
26*       LOGICAL            IEEE
27*       INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
28*       DOUBLE PRECISION   DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
29*      $                   QMAX, SIGMA, TAU
30*       ..
31*       .. Array Arguments ..
32*       DOUBLE PRECISION   Z( * )
33*       ..
34*
35*
36*> \par Purpose:
37*  =============
38*>
39*> \verbatim
40*>
41*> DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
42*> In case of failure it changes shifts, and tries again until output
43*> is positive.
44*> \endverbatim
45*
46*  Arguments:
47*  ==========
48*
49*> \param[in] I0
50*> \verbatim
51*>          I0 is INTEGER
52*>         First index.
53*> \endverbatim
54*>
55*> \param[in,out] N0
56*> \verbatim
57*>          N0 is INTEGER
58*>         Last index.
59*> \endverbatim
60*>
61*> \param[in,out] Z
62*> \verbatim
63*>          Z is DOUBLE PRECISION array, dimension ( 4*N0 )
64*>         Z holds the qd array.
65*> \endverbatim
66*>
67*> \param[in,out] PP
68*> \verbatim
69*>          PP is INTEGER
70*>         PP=0 for ping, PP=1 for pong.
71*>         PP=2 indicates that flipping was applied to the Z array
72*>         and that the initial tests for deflation should not be
73*>         performed.
74*> \endverbatim
75*>
76*> \param[out] DMIN
77*> \verbatim
78*>          DMIN is DOUBLE PRECISION
79*>         Minimum value of d.
80*> \endverbatim
81*>
82*> \param[out] SIGMA
83*> \verbatim
84*>          SIGMA is DOUBLE PRECISION
85*>         Sum of shifts used in current segment.
86*> \endverbatim
87*>
88*> \param[in,out] DESIG
89*> \verbatim
90*>          DESIG is DOUBLE PRECISION
91*>         Lower order part of SIGMA
92*> \endverbatim
93*>
94*> \param[in] QMAX
95*> \verbatim
96*>          QMAX is DOUBLE PRECISION
97*>         Maximum value of q.
98*> \endverbatim
99*>
100*> \param[in,out] NFAIL
101*> \verbatim
102*>          NFAIL is INTEGER
103*>         Increment NFAIL by 1 each time the shift was too big.
104*> \endverbatim
105*>
106*> \param[in,out] ITER
107*> \verbatim
108*>          ITER is INTEGER
109*>         Increment ITER by 1 for each iteration.
110*> \endverbatim
111*>
112*> \param[in,out] NDIV
113*> \verbatim
114*>          NDIV is INTEGER
115*>         Increment NDIV by 1 for each division.
116*> \endverbatim
117*>
118*> \param[in] IEEE
119*> \verbatim
120*>          IEEE is LOGICAL
121*>         Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
122*> \endverbatim
123*>
124*> \param[in,out] TTYPE
125*> \verbatim
126*>          TTYPE is INTEGER
127*>         Shift type.
128*> \endverbatim
129*>
130*> \param[in,out] DMIN1
131*> \verbatim
132*>          DMIN1 is DOUBLE PRECISION
133*> \endverbatim
134*>
135*> \param[in,out] DMIN2
136*> \verbatim
137*>          DMIN2 is DOUBLE PRECISION
138*> \endverbatim
139*>
140*> \param[in,out] DN
141*> \verbatim
142*>          DN is DOUBLE PRECISION
143*> \endverbatim
144*>
145*> \param[in,out] DN1
146*> \verbatim
147*>          DN1 is DOUBLE PRECISION
148*> \endverbatim
149*>
150*> \param[in,out] DN2
151*> \verbatim
152*>          DN2 is DOUBLE PRECISION
153*> \endverbatim
154*>
155*> \param[in,out] G
156*> \verbatim
157*>          G is DOUBLE PRECISION
158*> \endverbatim
159*>
160*> \param[in,out] TAU
161*> \verbatim
162*>          TAU is DOUBLE PRECISION
163*>
164*>         These are passed as arguments in order to save their values
165*>         between calls to DLASQ3.
166*> \endverbatim
167*
168*  Authors:
169*  ========
170*
171*> \author Univ. of Tennessee
172*> \author Univ. of California Berkeley
173*> \author Univ. of Colorado Denver
174*> \author NAG Ltd.
175*
176*> \date June 2016
177*
178*> \ingroup auxOTHERcomputational
179*
180*  =====================================================================
181      SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
182     $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
183     $                   DN2, G, TAU )
184*
185*  -- LAPACK computational routine (version 3.7.0) --
186*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
187*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188*     June 2016
189*
190*     .. Scalar Arguments ..
191      LOGICAL            IEEE
192      INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
193      DOUBLE PRECISION   DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
194     $                   QMAX, SIGMA, TAU
195*     ..
196*     .. Array Arguments ..
197      DOUBLE PRECISION   Z( * )
198*     ..
199*
200*  =====================================================================
201*
202*     .. Parameters ..
203      DOUBLE PRECISION   CBIAS
204      PARAMETER          ( CBIAS = 1.50D0 )
205      DOUBLE PRECISION   ZERO, QURTR, HALF, ONE, TWO, HUNDRD
206      PARAMETER          ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0,
207     $                     ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 )
208*     ..
209*     .. Local Scalars ..
210      INTEGER            IPN4, J4, N0IN, NN, TTYPE
211      DOUBLE PRECISION   EPS, S, T, TEMP, TOL, TOL2
212*     ..
213*     .. External Subroutines ..
214      EXTERNAL           DLASQ4, DLASQ5, DLASQ6
215*     ..
216*     .. External Function ..
217      DOUBLE PRECISION   DLAMCH
218      LOGICAL            DISNAN
219      EXTERNAL           DISNAN, DLAMCH
220*     ..
221*     .. Intrinsic Functions ..
222      INTRINSIC          ABS, MAX, MIN, SQRT
223*     ..
224*     .. Executable Statements ..
225*
226      N0IN = N0
227      EPS = DLAMCH( 'Precision' )
228      TOL = EPS*HUNDRD
229      TOL2 = TOL**2
230*
231*     Check for deflation.
232*
233   10 CONTINUE
234*
235      IF( N0.LT.I0 )
236     $   RETURN
237      IF( N0.EQ.I0 )
238     $   GO TO 20
239      NN = 4*N0 + PP
240      IF( N0.EQ.( I0+1 ) )
241     $   GO TO 40
242*
243*     Check whether E(N0-1) is negligible, 1 eigenvalue.
244*
245      IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
246     $    Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
247     $   GO TO 30
248*
249   20 CONTINUE
250*
251      Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
252      N0 = N0 - 1
253      GO TO 10
254*
255*     Check  whether E(N0-2) is negligible, 2 eigenvalues.
256*
257   30 CONTINUE
258*
259      IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
260     $    Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
261     $   GO TO 50
262*
263   40 CONTINUE
264*
265      IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
266         S = Z( NN-3 )
267         Z( NN-3 ) = Z( NN-7 )
268         Z( NN-7 ) = S
269      END IF
270      T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
271      IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2.AND.T.NE.ZERO ) THEN
272         S = Z( NN-3 )*( Z( NN-5 ) / T )
273         IF( S.LE.T ) THEN
274            S = Z( NN-3 )*( Z( NN-5 ) /
275     $          ( T*( ONE+SQRT( ONE+S / T ) ) ) )
276         ELSE
277            S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
278         END IF
279         T = Z( NN-7 ) + ( S+Z( NN-5 ) )
280         Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
281         Z( NN-7 ) = T
282      END IF
283      Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
284      Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
285      N0 = N0 - 2
286      GO TO 10
287*
288   50 CONTINUE
289      IF( PP.EQ.2 )
290     $   PP = 0
291*
292*     Reverse the qd-array, if warranted.
293*
294      IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
295         IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
296            IPN4 = 4*( I0+N0 )
297            DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
298               TEMP = Z( J4-3 )
299               Z( J4-3 ) = Z( IPN4-J4-3 )
300               Z( IPN4-J4-3 ) = TEMP
301               TEMP = Z( J4-2 )
302               Z( J4-2 ) = Z( IPN4-J4-2 )
303               Z( IPN4-J4-2 ) = TEMP
304               TEMP = Z( J4-1 )
305               Z( J4-1 ) = Z( IPN4-J4-5 )
306               Z( IPN4-J4-5 ) = TEMP
307               TEMP = Z( J4 )
308               Z( J4 ) = Z( IPN4-J4-4 )
309               Z( IPN4-J4-4 ) = TEMP
310   60       CONTINUE
311            IF( N0-I0.LE.4 ) THEN
312               Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
313               Z( 4*N0-PP ) = Z( 4*I0-PP )
314            END IF
315            DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
316            Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
317     $                            Z( 4*I0+PP+3 ) )
318            Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
319     $                          Z( 4*I0-PP+4 ) )
320            QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
321            DMIN = -ZERO
322         END IF
323      END IF
324*
325*     Choose a shift.
326*
327      CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
328     $             DN2, TAU, TTYPE, G )
329*
330*     Call dqds until DMIN > 0.
331*
332   70 CONTINUE
333*
334      CALL DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
335     $             DN1, DN2, IEEE, EPS )
336*
337      NDIV = NDIV + ( N0-I0+2 )
338      ITER = ITER + 1
339*
340*     Check status.
341*
342      IF( DMIN.GE.ZERO .AND. DMIN1.GE.ZERO ) THEN
343*
344*        Success.
345*
346         GO TO 90
347*
348      ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND.
349     $         Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
350     $         ABS( DN ).LT.TOL*SIGMA ) THEN
351*
352*        Convergence hidden by negative DN.
353*
354         Z( 4*( N0-1 )-PP+2 ) = ZERO
355         DMIN = ZERO
356         GO TO 90
357      ELSE IF( DMIN.LT.ZERO ) THEN
358*
359*        TAU too big. Select new TAU and try again.
360*
361         NFAIL = NFAIL + 1
362         IF( TTYPE.LT.-22 ) THEN
363*
364*           Failed twice. Play it safe.
365*
366            TAU = ZERO
367         ELSE IF( DMIN1.GT.ZERO ) THEN
368*
369*           Late failure. Gives excellent shift.
370*
371            TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
372            TTYPE = TTYPE - 11
373         ELSE
374*
375*           Early failure. Divide by 4.
376*
377            TAU = QURTR*TAU
378            TTYPE = TTYPE - 12
379         END IF
380         GO TO 70
381      ELSE IF( DISNAN( DMIN ) ) THEN
382*
383*        NaN.
384*
385         IF( TAU.EQ.ZERO ) THEN
386            GO TO 80
387         ELSE
388            TAU = ZERO
389            GO TO 70
390         END IF
391      ELSE
392*
393*        Possible underflow. Play it safe.
394*
395         GO TO 80
396      END IF
397*
398*     Risk of underflow.
399*
400   80 CONTINUE
401      CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
402      NDIV = NDIV + ( N0-I0+2 )
403      ITER = ITER + 1
404      TAU = ZERO
405*
406   90 CONTINUE
407      IF( TAU.LT.SIGMA ) THEN
408         DESIG = DESIG + TAU
409         T = SIGMA + DESIG
410         DESIG = DESIG - ( T-SIGMA )
411      ELSE
412         T = SIGMA + TAU
413         DESIG = SIGMA - ( T-TAU ) + DESIG
414      END IF
415      SIGMA = T
416*
417      RETURN
418*
419*     End of DLASQ3
420*
421      END
422