1*> \brief \b SLASQ3 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 SLASQ3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SLASQ3( 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*       REAL               DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
29*      $                   QMAX, SIGMA, TAU
30*       ..
31*       .. Array Arguments ..
32*       REAL               Z( * )
33*       ..
34*
35*
36*> \par Purpose:
37*  =============
38*>
39*> \verbatim
40*>
41*> SLASQ3 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 REAL 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 REAL
79*>         Minimum value of d.
80*> \endverbatim
81*>
82*> \param[out] SIGMA
83*> \verbatim
84*>          SIGMA is REAL
85*>         Sum of shifts used in current segment.
86*> \endverbatim
87*>
88*> \param[in,out] DESIG
89*> \verbatim
90*>          DESIG is REAL
91*>         Lower order part of SIGMA
92*> \endverbatim
93*>
94*> \param[in] QMAX
95*> \verbatim
96*>          QMAX is REAL
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 SLASQ5).
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 REAL
133*> \endverbatim
134*>
135*> \param[in,out] DMIN2
136*> \verbatim
137*>          DMIN2 is REAL
138*> \endverbatim
139*>
140*> \param[in,out] DN
141*> \verbatim
142*>          DN is REAL
143*> \endverbatim
144*>
145*> \param[in,out] DN1
146*> \verbatim
147*>          DN1 is REAL
148*> \endverbatim
149*>
150*> \param[in,out] DN2
151*> \verbatim
152*>          DN2 is REAL
153*> \endverbatim
154*>
155*> \param[in,out] G
156*> \verbatim
157*>          G is REAL
158*> \endverbatim
159*>
160*> \param[in,out] TAU
161*> \verbatim
162*>          TAU is REAL
163*>
164*>         These are passed as arguments in order to save their values
165*>         between calls to SLASQ3.
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*> \ingroup auxOTHERcomputational
177*
178*  =====================================================================
179      SUBROUTINE SLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
180     $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
181     $                   DN2, G, TAU )
182*
183*  -- LAPACK computational routine --
184*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
185*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186*
187*     .. Scalar Arguments ..
188      LOGICAL            IEEE
189      INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
190      REAL               DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
191     $                   QMAX, SIGMA, TAU
192*     ..
193*     .. Array Arguments ..
194      REAL               Z( * )
195*     ..
196*
197*  =====================================================================
198*
199*     .. Parameters ..
200      REAL               CBIAS
201      PARAMETER          ( CBIAS = 1.50E0 )
202      REAL               ZERO, QURTR, HALF, ONE, TWO, HUNDRD
203      PARAMETER          ( ZERO = 0.0E0, QURTR = 0.250E0, HALF = 0.5E0,
204     $                     ONE = 1.0E0, TWO = 2.0E0, HUNDRD = 100.0E0 )
205*     ..
206*     .. Local Scalars ..
207      INTEGER            IPN4, J4, N0IN, NN, TTYPE
208      REAL               EPS, S, T, TEMP, TOL, TOL2
209*     ..
210*     .. External Subroutines ..
211      EXTERNAL           SLASQ4, SLASQ5, SLASQ6
212*     ..
213*     .. External Function ..
214      REAL               SLAMCH
215      LOGICAL            SISNAN
216      EXTERNAL           SISNAN, SLAMCH
217*     ..
218*     .. Intrinsic Functions ..
219      INTRINSIC          ABS, MAX, MIN, SQRT
220*     ..
221*     .. Executable Statements ..
222*
223      N0IN = N0
224      EPS = SLAMCH( 'Precision' )
225      TOL = EPS*HUNDRD
226      TOL2 = TOL**2
227*
228*     Check for deflation.
229*
230   10 CONTINUE
231*
232      IF( N0.LT.I0 )
233     $   RETURN
234      IF( N0.EQ.I0 )
235     $   GO TO 20
236      NN = 4*N0 + PP
237      IF( N0.EQ.( I0+1 ) )
238     $   GO TO 40
239*
240*     Check whether E(N0-1) is negligible, 1 eigenvalue.
241*
242      IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
243     $    Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
244     $   GO TO 30
245*
246   20 CONTINUE
247*
248      Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
249      N0 = N0 - 1
250      GO TO 10
251*
252*     Check  whether E(N0-2) is negligible, 2 eigenvalues.
253*
254   30 CONTINUE
255*
256      IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
257     $    Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
258     $   GO TO 50
259*
260   40 CONTINUE
261*
262      IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
263         S = Z( NN-3 )
264         Z( NN-3 ) = Z( NN-7 )
265         Z( NN-7 ) = S
266      END IF
267      T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
268      IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2.AND.T.NE.ZERO ) THEN
269         S = Z( NN-3 )*( Z( NN-5 ) / T )
270         IF( S.LE.T ) THEN
271            S = Z( NN-3 )*( Z( NN-5 ) /
272     $          ( T*( ONE+SQRT( ONE+S / T ) ) ) )
273         ELSE
274            S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
275         END IF
276         T = Z( NN-7 ) + ( S+Z( NN-5 ) )
277         Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
278         Z( NN-7 ) = T
279      END IF
280      Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
281      Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
282      N0 = N0 - 2
283      GO TO 10
284*
285   50 CONTINUE
286      IF( PP.EQ.2 )
287     $   PP = 0
288*
289*     Reverse the qd-array, if warranted.
290*
291      IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
292         IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
293            IPN4 = 4*( I0+N0 )
294            DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
295               TEMP = Z( J4-3 )
296               Z( J4-3 ) = Z( IPN4-J4-3 )
297               Z( IPN4-J4-3 ) = TEMP
298               TEMP = Z( J4-2 )
299               Z( J4-2 ) = Z( IPN4-J4-2 )
300               Z( IPN4-J4-2 ) = TEMP
301               TEMP = Z( J4-1 )
302               Z( J4-1 ) = Z( IPN4-J4-5 )
303               Z( IPN4-J4-5 ) = TEMP
304               TEMP = Z( J4 )
305               Z( J4 ) = Z( IPN4-J4-4 )
306               Z( IPN4-J4-4 ) = TEMP
307   60       CONTINUE
308            IF( N0-I0.LE.4 ) THEN
309               Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
310               Z( 4*N0-PP ) = Z( 4*I0-PP )
311            END IF
312            DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
313            Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
314     $                            Z( 4*I0+PP+3 ) )
315            Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
316     $                          Z( 4*I0-PP+4 ) )
317            QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
318            DMIN = -ZERO
319         END IF
320      END IF
321*
322*     Choose a shift.
323*
324      CALL SLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
325     $             DN2, TAU, TTYPE, G )
326*
327*     Call dqds until DMIN > 0.
328*
329   70 CONTINUE
330*
331      CALL SLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
332     $             DN1, DN2, IEEE, EPS )
333*
334      NDIV = NDIV + ( N0-I0+2 )
335      ITER = ITER + 1
336*
337*     Check status.
338*
339      IF( DMIN.GE.ZERO .AND. DMIN1.GE.ZERO ) THEN
340*
341*        Success.
342*
343         GO TO 90
344*
345      ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND.
346     $         Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
347     $         ABS( DN ).LT.TOL*SIGMA ) THEN
348*
349*        Convergence hidden by negative DN.
350*
351         Z( 4*( N0-1 )-PP+2 ) = ZERO
352         DMIN = ZERO
353         GO TO 90
354      ELSE IF( DMIN.LT.ZERO ) THEN
355*
356*        TAU too big. Select new TAU and try again.
357*
358         NFAIL = NFAIL + 1
359         IF( TTYPE.LT.-22 ) THEN
360*
361*           Failed twice. Play it safe.
362*
363            TAU = ZERO
364         ELSE IF( DMIN1.GT.ZERO ) THEN
365*
366*           Late failure. Gives excellent shift.
367*
368            TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
369            TTYPE = TTYPE - 11
370         ELSE
371*
372*           Early failure. Divide by 4.
373*
374            TAU = QURTR*TAU
375            TTYPE = TTYPE - 12
376         END IF
377         GO TO 70
378      ELSE IF( SISNAN( DMIN ) ) THEN
379*
380*        NaN.
381*
382         IF( TAU.EQ.ZERO ) THEN
383            GO TO 80
384         ELSE
385            TAU = ZERO
386            GO TO 70
387         END IF
388      ELSE
389*
390*        Possible underflow. Play it safe.
391*
392         GO TO 80
393      END IF
394*
395*     Risk of underflow.
396*
397   80 CONTINUE
398      CALL SLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
399      NDIV = NDIV + ( N0-I0+2 )
400      ITER = ITER + 1
401      TAU = ZERO
402*
403   90 CONTINUE
404      IF( TAU.LT.SIGMA ) THEN
405         DESIG = DESIG + TAU
406         T = SIGMA + DESIG
407         DESIG = DESIG - ( T-SIGMA )
408      ELSE
409         T = SIGMA + TAU
410         DESIG = SIGMA - ( T-TAU ) + DESIG
411      END IF
412      SIGMA = T
413*
414      RETURN
415*
416*     End of SLASQ3
417*
418      END
419