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*> \ingroup auxOTHERcomputational
140*
141*  =====================================================================
142      SUBROUTINE SLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
143     $                   DN, DNM1, DNM2, IEEE, EPS )
144*
145*  -- LAPACK computational routine --
146*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
147*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148*
149*     .. Scalar Arguments ..
150      LOGICAL            IEEE
151      INTEGER            I0, N0, PP
152      REAL               DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
153     $                   SIGMA, EPS
154*     ..
155*     .. Array Arguments ..
156      REAL               Z( * )
157*     ..
158*
159*  =====================================================================
160*
161*     .. Parameter ..
162      REAL               ZERO, HALF
163      PARAMETER          ( ZERO = 0.0E0, HALF = 0.5 )
164*     ..
165*     .. Local Scalars ..
166      INTEGER            J4, J4P2
167      REAL               D, EMIN, TEMP, DTHRESH
168*     ..
169*     .. Intrinsic Functions ..
170      INTRINSIC          MIN
171*     ..
172*     .. Executable Statements ..
173*
174      IF( ( N0-I0-1 ).LE.0 )
175     $   RETURN
176*
177      DTHRESH = EPS*(SIGMA+TAU)
178      IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
179      IF( TAU.NE.ZERO ) THEN
180         J4 = 4*I0 + PP - 3
181         EMIN = Z( J4+4 )
182         D = Z( J4 ) - TAU
183         DMIN = D
184         DMIN1 = -Z( J4 )
185*
186         IF( IEEE ) THEN
187*
188*     Code for IEEE arithmetic.
189*
190            IF( PP.EQ.0 ) THEN
191               DO 10 J4 = 4*I0, 4*( N0-3 ), 4
192                  Z( J4-2 ) = D + Z( J4-1 )
193                  TEMP = Z( J4+1 ) / Z( J4-2 )
194                  D = D*TEMP - TAU
195                  DMIN = MIN( DMIN, D )
196                  Z( J4 ) = Z( J4-1 )*TEMP
197                  EMIN = MIN( Z( J4 ), EMIN )
198 10            CONTINUE
199            ELSE
200               DO 20 J4 = 4*I0, 4*( N0-3 ), 4
201                  Z( J4-3 ) = D + Z( J4 )
202                  TEMP = Z( J4+2 ) / Z( J4-3 )
203                  D = D*TEMP - TAU
204                  DMIN = MIN( DMIN, D )
205                  Z( J4-1 ) = Z( J4 )*TEMP
206                  EMIN = MIN( Z( J4-1 ), EMIN )
207 20            CONTINUE
208            END IF
209*
210*     Unroll last two steps.
211*
212            DNM2 = D
213            DMIN2 = DMIN
214            J4 = 4*( N0-2 ) - PP
215            J4P2 = J4 + 2*PP - 1
216            Z( J4-2 ) = DNM2 + Z( J4P2 )
217            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
218            DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
219            DMIN = MIN( DMIN, DNM1 )
220*
221            DMIN1 = DMIN
222            J4 = J4 + 4
223            J4P2 = J4 + 2*PP - 1
224            Z( J4-2 ) = DNM1 + Z( J4P2 )
225            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
226            DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
227            DMIN = MIN( DMIN, DN )
228*
229         ELSE
230*
231*     Code for non IEEE arithmetic.
232*
233            IF( PP.EQ.0 ) THEN
234               DO 30 J4 = 4*I0, 4*( N0-3 ), 4
235                  Z( J4-2 ) = D + Z( J4-1 )
236                  IF( D.LT.ZERO ) THEN
237                     RETURN
238                  ELSE
239                     Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
240                     D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
241                  END IF
242                  DMIN = MIN( DMIN, D )
243                  EMIN = MIN( EMIN, Z( J4 ) )
244 30            CONTINUE
245            ELSE
246               DO 40 J4 = 4*I0, 4*( N0-3 ), 4
247                  Z( J4-3 ) = D + Z( J4 )
248                  IF( D.LT.ZERO ) THEN
249                     RETURN
250                  ELSE
251                     Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
252                     D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
253                  END IF
254                  DMIN = MIN( DMIN, D )
255                  EMIN = MIN( EMIN, Z( J4-1 ) )
256 40            CONTINUE
257            END IF
258*
259*     Unroll last two steps.
260*
261            DNM2 = D
262            DMIN2 = DMIN
263            J4 = 4*( N0-2 ) - PP
264            J4P2 = J4 + 2*PP - 1
265            Z( J4-2 ) = DNM2 + Z( J4P2 )
266            IF( DNM2.LT.ZERO ) THEN
267               RETURN
268            ELSE
269               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
270               DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
271            END IF
272            DMIN = MIN( DMIN, DNM1 )
273*
274            DMIN1 = DMIN
275            J4 = J4 + 4
276            J4P2 = J4 + 2*PP - 1
277            Z( J4-2 ) = DNM1 + Z( J4P2 )
278            IF( DNM1.LT.ZERO ) THEN
279               RETURN
280            ELSE
281               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
282               DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
283            END IF
284            DMIN = MIN( DMIN, DN )
285*
286         END IF
287*
288      ELSE
289*     This is the version that sets d's to zero if they are small enough
290         J4 = 4*I0 + PP - 3
291         EMIN = Z( J4+4 )
292         D = Z( J4 ) - TAU
293         DMIN = D
294         DMIN1 = -Z( J4 )
295         IF( IEEE ) THEN
296*
297*     Code for IEEE arithmetic.
298*
299            IF( PP.EQ.0 ) THEN
300               DO 50 J4 = 4*I0, 4*( N0-3 ), 4
301                  Z( J4-2 ) = D + Z( J4-1 )
302                  TEMP = Z( J4+1 ) / Z( J4-2 )
303                  D = D*TEMP - TAU
304                  IF( D.LT.DTHRESH ) D = ZERO
305                  DMIN = MIN( DMIN, D )
306                  Z( J4 ) = Z( J4-1 )*TEMP
307                  EMIN = MIN( Z( J4 ), EMIN )
308 50            CONTINUE
309            ELSE
310               DO 60 J4 = 4*I0, 4*( N0-3 ), 4
311                  Z( J4-3 ) = D + Z( J4 )
312                  TEMP = Z( J4+2 ) / Z( J4-3 )
313                  D = D*TEMP - TAU
314                  IF( D.LT.DTHRESH ) D = ZERO
315                  DMIN = MIN( DMIN, D )
316                  Z( J4-1 ) = Z( J4 )*TEMP
317                  EMIN = MIN( Z( J4-1 ), EMIN )
318 60            CONTINUE
319            END IF
320*
321*     Unroll last two steps.
322*
323            DNM2 = D
324            DMIN2 = DMIN
325            J4 = 4*( N0-2 ) - PP
326            J4P2 = J4 + 2*PP - 1
327            Z( J4-2 ) = DNM2 + Z( J4P2 )
328            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
329            DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
330            DMIN = MIN( DMIN, DNM1 )
331*
332            DMIN1 = DMIN
333            J4 = J4 + 4
334            J4P2 = J4 + 2*PP - 1
335            Z( J4-2 ) = DNM1 + Z( J4P2 )
336            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
337            DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
338            DMIN = MIN( DMIN, DN )
339*
340         ELSE
341*
342*     Code for non IEEE arithmetic.
343*
344            IF( PP.EQ.0 ) THEN
345               DO 70 J4 = 4*I0, 4*( N0-3 ), 4
346                  Z( J4-2 ) = D + Z( J4-1 )
347                  IF( D.LT.ZERO ) THEN
348                     RETURN
349                  ELSE
350                     Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
351                     D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
352                  END IF
353                  IF( D.LT.DTHRESH ) D = ZERO
354                  DMIN = MIN( DMIN, D )
355                  EMIN = MIN( EMIN, Z( J4 ) )
356 70            CONTINUE
357            ELSE
358               DO 80 J4 = 4*I0, 4*( N0-3 ), 4
359                  Z( J4-3 ) = D + Z( J4 )
360                  IF( D.LT.ZERO ) THEN
361                     RETURN
362                  ELSE
363                     Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
364                     D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
365                  END IF
366                  IF( D.LT.DTHRESH ) D = ZERO
367                  DMIN = MIN( DMIN, D )
368                  EMIN = MIN( EMIN, Z( J4-1 ) )
369 80            CONTINUE
370            END IF
371*
372*     Unroll last two steps.
373*
374            DNM2 = D
375            DMIN2 = DMIN
376            J4 = 4*( N0-2 ) - PP
377            J4P2 = J4 + 2*PP - 1
378            Z( J4-2 ) = DNM2 + Z( J4P2 )
379            IF( DNM2.LT.ZERO ) THEN
380               RETURN
381            ELSE
382               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
383               DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
384            END IF
385            DMIN = MIN( DMIN, DNM1 )
386*
387            DMIN1 = DMIN
388            J4 = J4 + 4
389            J4P2 = J4 + 2*PP - 1
390            Z( J4-2 ) = DNM1 + Z( J4P2 )
391            IF( DNM1.LT.ZERO ) THEN
392               RETURN
393            ELSE
394               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
395               DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
396            END IF
397            DMIN = MIN( DMIN, DN )
398*
399         END IF
400*
401      END IF
402      Z( J4+2 ) = DN
403      Z( 4*N0-PP ) = EMIN
404      RETURN
405*
406*     End of SLASQ5
407*
408      END
409