1*> \brief <b> SLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr. </b>
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SLASQ5 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq5.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq5.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq5.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
22*                          DNM1, DNM2, IEEE, EPS )
23*
24*       .. Scalar Arguments ..
25*       LOGICAL            IEEE
26*       INTEGER            I0, N0, PP
27*       REAL               EPS, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, SIGMA, TAU
28*       ..
29*       .. Array Arguments ..
30*       REAL               Z( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> SLASQ5 computes one dqds transform in ping-pong form, one
40*> version for IEEE machines another for non IEEE machines.
41*> \endverbatim
42*
43*  Arguments:
44*  ==========
45*
46*> \param[in] I0
47*> \verbatim
48*>          I0 is INTEGER
49*>        First index.
50*> \endverbatim
51*>
52*> \param[in] N0
53*> \verbatim
54*>          N0 is INTEGER
55*>        Last index.
56*> \endverbatim
57*>
58*> \param[in] Z
59*> \verbatim
60*>          Z is REAL array, dimension ( 4*N )
61*>        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
62*>        an extra argument.
63*> \endverbatim
64*>
65*> \param[in] PP
66*> \verbatim
67*>          PP is INTEGER
68*>        PP=0 for ping, PP=1 for pong.
69*> \endverbatim
70*>
71*> \param[in] TAU
72*> \verbatim
73*>          TAU is REAL
74*>        This is the shift.
75*> \endverbatim
76*>
77*> \param[in] SIGMA
78*> \verbatim
79*>          SIGMA is REAL
80*>        This is the accumulated shift up to this step.
81*> \endverbatim
82*>
83*> \param[out] DMIN
84*> \verbatim
85*>          DMIN is REAL
86*>        Minimum value of d.
87*> \endverbatim
88*>
89*> \param[out] DMIN1
90*> \verbatim
91*>          DMIN1 is REAL
92*>        Minimum value of d, excluding D( N0 ).
93*> \endverbatim
94*>
95*> \param[out] DMIN2
96*> \verbatim
97*>          DMIN2 is REAL
98*>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
99*> \endverbatim
100*>
101*> \param[out] DN
102*> \verbatim
103*>          DN is REAL
104*>        d(N0), the last value of d.
105*> \endverbatim
106*>
107*> \param[out] DNM1
108*> \verbatim
109*>          DNM1 is REAL
110*>        d(N0-1).
111*> \endverbatim
112*>
113*> \param[out] DNM2
114*> \verbatim
115*>          DNM2 is REAL
116*>        d(N0-2).
117*> \endverbatim
118*>
119*> \param[in] IEEE
120*> \verbatim
121*>          IEEE is LOGICAL
122*>        Flag for IEEE or non IEEE arithmetic.
123*> \endverbatim
124*>
125*> \param[in] EPS
126*> \verbatim
127*>         EPS is REAL
128*>        This is the value of epsilon used.
129*> \endverbatim
130*
131*  Authors:
132*  ========
133*
134*> \author Univ. of Tennessee
135*> \author Univ. of California Berkeley
136*> \author Univ. of Colorado Denver
137*> \author NAG Ltd.
138*
139*> \date December 2016
140*
141*> \ingroup auxOTHERcomputational
142*
143*  =====================================================================
144      SUBROUTINE SLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
145     $                   DN, DNM1, DNM2, IEEE, EPS )
146*
147*  -- LAPACK computational routine (version 3.7.0) --
148*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
149*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
150*     December 2016
151*
152*     .. Scalar Arguments ..
153      LOGICAL            IEEE
154      INTEGER            I0, N0, PP
155      REAL               DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
156     $                   SIGMA, EPS
157*     ..
158*     .. Array Arguments ..
159      REAL               Z( * )
160*     ..
161*
162*  =====================================================================
163*
164*     .. Parameter ..
165      REAL               ZERO, HALF
166      PARAMETER          ( ZERO = 0.0E0, HALF = 0.5 )
167*     ..
168*     .. Local Scalars ..
169      INTEGER            J4, J4P2
170      REAL               D, EMIN, TEMP, DTHRESH
171*     ..
172*     .. Intrinsic Functions ..
173      INTRINSIC          MIN
174*     ..
175*     .. Executable Statements ..
176*
177      IF( ( N0-I0-1 ).LE.0 )
178     $   RETURN
179*
180      DTHRESH = EPS*(SIGMA+TAU)
181      IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
182      IF( TAU.NE.ZERO ) THEN
183         J4 = 4*I0 + PP - 3
184         EMIN = Z( J4+4 )
185         D = Z( J4 ) - TAU
186         DMIN = D
187         DMIN1 = -Z( J4 )
188*
189         IF( IEEE ) THEN
190*
191*     Code for IEEE arithmetic.
192*
193            IF( PP.EQ.0 ) THEN
194               DO 10 J4 = 4*I0, 4*( N0-3 ), 4
195                  Z( J4-2 ) = D + Z( J4-1 )
196                  TEMP = Z( J4+1 ) / Z( J4-2 )
197                  D = D*TEMP - TAU
198                  DMIN = MIN( DMIN, D )
199                  Z( J4 ) = Z( J4-1 )*TEMP
200                  EMIN = MIN( Z( J4 ), EMIN )
201 10            CONTINUE
202            ELSE
203               DO 20 J4 = 4*I0, 4*( N0-3 ), 4
204                  Z( J4-3 ) = D + Z( J4 )
205                  TEMP = Z( J4+2 ) / Z( J4-3 )
206                  D = D*TEMP - TAU
207                  DMIN = MIN( DMIN, D )
208                  Z( J4-1 ) = Z( J4 )*TEMP
209                  EMIN = MIN( Z( J4-1 ), EMIN )
210 20            CONTINUE
211            END IF
212*
213*     Unroll last two steps.
214*
215            DNM2 = D
216            DMIN2 = DMIN
217            J4 = 4*( N0-2 ) - PP
218            J4P2 = J4 + 2*PP - 1
219            Z( J4-2 ) = DNM2 + Z( J4P2 )
220            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
221            DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
222            DMIN = MIN( DMIN, DNM1 )
223*
224            DMIN1 = DMIN
225            J4 = J4 + 4
226            J4P2 = J4 + 2*PP - 1
227            Z( J4-2 ) = DNM1 + Z( J4P2 )
228            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
229            DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
230            DMIN = MIN( DMIN, DN )
231*
232         ELSE
233*
234*     Code for non IEEE arithmetic.
235*
236            IF( PP.EQ.0 ) THEN
237               DO 30 J4 = 4*I0, 4*( N0-3 ), 4
238                  Z( J4-2 ) = D + Z( J4-1 )
239                  IF( D.LT.ZERO ) THEN
240                     RETURN
241                  ELSE
242                     Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
243                     D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
244                  END IF
245                  DMIN = MIN( DMIN, D )
246                  EMIN = MIN( EMIN, Z( J4 ) )
247 30            CONTINUE
248            ELSE
249               DO 40 J4 = 4*I0, 4*( N0-3 ), 4
250                  Z( J4-3 ) = D + Z( J4 )
251                  IF( D.LT.ZERO ) THEN
252                     RETURN
253                  ELSE
254                     Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
255                     D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
256                  END IF
257                  DMIN = MIN( DMIN, D )
258                  EMIN = MIN( EMIN, Z( J4-1 ) )
259 40            CONTINUE
260            END IF
261*
262*     Unroll last two steps.
263*
264            DNM2 = D
265            DMIN2 = DMIN
266            J4 = 4*( N0-2 ) - PP
267            J4P2 = J4 + 2*PP - 1
268            Z( J4-2 ) = DNM2 + Z( J4P2 )
269            IF( DNM2.LT.ZERO ) THEN
270               RETURN
271            ELSE
272               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
273               DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
274            END IF
275            DMIN = MIN( DMIN, DNM1 )
276*
277            DMIN1 = DMIN
278            J4 = J4 + 4
279            J4P2 = J4 + 2*PP - 1
280            Z( J4-2 ) = DNM1 + Z( J4P2 )
281            IF( DNM1.LT.ZERO ) THEN
282               RETURN
283            ELSE
284               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
285               DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
286            END IF
287            DMIN = MIN( DMIN, DN )
288*
289         END IF
290*
291      ELSE
292*     This is the version that sets d's to zero if they are small enough
293         J4 = 4*I0 + PP - 3
294         EMIN = Z( J4+4 )
295         D = Z( J4 ) - TAU
296         DMIN = D
297         DMIN1 = -Z( J4 )
298         IF( IEEE ) THEN
299*
300*     Code for IEEE arithmetic.
301*
302            IF( PP.EQ.0 ) THEN
303               DO 50 J4 = 4*I0, 4*( N0-3 ), 4
304                  Z( J4-2 ) = D + Z( J4-1 )
305                  TEMP = Z( J4+1 ) / Z( J4-2 )
306                  D = D*TEMP - TAU
307                  IF( D.LT.DTHRESH ) D = ZERO
308                  DMIN = MIN( DMIN, D )
309                  Z( J4 ) = Z( J4-1 )*TEMP
310                  EMIN = MIN( Z( J4 ), EMIN )
311 50            CONTINUE
312            ELSE
313               DO 60 J4 = 4*I0, 4*( N0-3 ), 4
314                  Z( J4-3 ) = D + Z( J4 )
315                  TEMP = Z( J4+2 ) / Z( J4-3 )
316                  D = D*TEMP - TAU
317                  IF( D.LT.DTHRESH ) D = ZERO
318                  DMIN = MIN( DMIN, D )
319                  Z( J4-1 ) = Z( J4 )*TEMP
320                  EMIN = MIN( Z( J4-1 ), EMIN )
321 60            CONTINUE
322            END IF
323*
324*     Unroll last two steps.
325*
326            DNM2 = D
327            DMIN2 = DMIN
328            J4 = 4*( N0-2 ) - PP
329            J4P2 = J4 + 2*PP - 1
330            Z( J4-2 ) = DNM2 + Z( J4P2 )
331            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
332            DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
333            DMIN = MIN( DMIN, DNM1 )
334*
335            DMIN1 = DMIN
336            J4 = J4 + 4
337            J4P2 = J4 + 2*PP - 1
338            Z( J4-2 ) = DNM1 + Z( J4P2 )
339            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
340            DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
341            DMIN = MIN( DMIN, DN )
342*
343         ELSE
344*
345*     Code for non IEEE arithmetic.
346*
347            IF( PP.EQ.0 ) THEN
348               DO 70 J4 = 4*I0, 4*( N0-3 ), 4
349                  Z( J4-2 ) = D + Z( J4-1 )
350                  IF( D.LT.ZERO ) THEN
351                     RETURN
352                  ELSE
353                     Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
354                     D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
355                  END IF
356                  IF( D.LT.DTHRESH ) D = ZERO
357                  DMIN = MIN( DMIN, D )
358                  EMIN = MIN( EMIN, Z( J4 ) )
359 70            CONTINUE
360            ELSE
361               DO 80 J4 = 4*I0, 4*( N0-3 ), 4
362                  Z( J4-3 ) = D + Z( J4 )
363                  IF( D.LT.ZERO ) THEN
364                     RETURN
365                  ELSE
366                     Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
367                     D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
368                  END IF
369                  IF( D.LT.DTHRESH ) D = ZERO
370                  DMIN = MIN( DMIN, D )
371                  EMIN = MIN( EMIN, Z( J4-1 ) )
372 80            CONTINUE
373            END IF
374*
375*     Unroll last two steps.
376*
377            DNM2 = D
378            DMIN2 = DMIN
379            J4 = 4*( N0-2 ) - PP
380            J4P2 = J4 + 2*PP - 1
381            Z( J4-2 ) = DNM2 + Z( J4P2 )
382            IF( DNM2.LT.ZERO ) THEN
383               RETURN
384            ELSE
385               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
386               DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
387            END IF
388            DMIN = MIN( DMIN, DNM1 )
389*
390            DMIN1 = DMIN
391            J4 = J4 + 4
392            J4P2 = J4 + 2*PP - 1
393            Z( J4-2 ) = DNM1 + Z( J4P2 )
394            IF( DNM1.LT.ZERO ) THEN
395               RETURN
396            ELSE
397               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
398               DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
399            END IF
400            DMIN = MIN( DMIN, DN )
401*
402         END IF
403*
404      END IF
405      Z( J4+2 ) = DN
406      Z( 4*N0-PP ) = EMIN
407      RETURN
408*
409*     End of SLASQ5
410*
411      END
412