1*> \brief \b DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLASQ5 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq5.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq5.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq5.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLASQ5( 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*       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, SIGMA, EPS
28*       ..
29*       .. Array Arguments ..
30*       DOUBLE PRECISION   Z( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> DLASQ5 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 DOUBLE PRECISION 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 DOUBLE PRECISION
74*>        This is the shift.
75*> \endverbatim
76*>
77*> \param[in] SIGMA
78*> \verbatim
79*>          SIGMA is DOUBLE PRECISION
80*>        This is the accumulated shift up to this step.
81*> \endverbatim
82*>
83*> \param[out] DMIN
84*> \verbatim
85*>          DMIN is DOUBLE PRECISION
86*>        Minimum value of d.
87*> \endverbatim
88*>
89*> \param[out] DMIN1
90*> \verbatim
91*>          DMIN1 is DOUBLE PRECISION
92*>        Minimum value of d, excluding D( N0 ).
93*> \endverbatim
94*>
95*> \param[out] DMIN2
96*> \verbatim
97*>          DMIN2 is DOUBLE PRECISION
98*>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
99*> \endverbatim
100*>
101*> \param[out] DN
102*> \verbatim
103*>          DN is DOUBLE PRECISION
104*>        d(N0), the last value of d.
105*> \endverbatim
106*>
107*> \param[out] DNM1
108*> \verbatim
109*>          DNM1 is DOUBLE PRECISION
110*>        d(N0-1).
111*> \endverbatim
112*>
113*> \param[out] DNM2
114*> \verbatim
115*>          DNM2 is DOUBLE PRECISION
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 DOUBLE PRECISION
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 DLASQ5( 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      DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
153     $                   SIGMA, EPS
154*     ..
155*     .. Array Arguments ..
156      DOUBLE PRECISION   Z( * )
157*     ..
158*
159*  =====================================================================
160*
161*     .. Parameter ..
162      DOUBLE PRECISION   ZERO, HALF
163      PARAMETER          ( ZERO = 0.0D0, HALF = 0.5 )
164*     ..
165*     .. Local Scalars ..
166      INTEGER            J4, J4P2
167      DOUBLE PRECISION   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      ELSE
288*     This is the version that sets d's to zero if they are small enough
289         J4 = 4*I0 + PP - 3
290         EMIN = Z( J4+4 )
291         D = Z( J4 ) - TAU
292         DMIN = D
293         DMIN1 = -Z( J4 )
294         IF( IEEE ) THEN
295*
296*     Code for IEEE arithmetic.
297*
298            IF( PP.EQ.0 ) THEN
299               DO 50 J4 = 4*I0, 4*( N0-3 ), 4
300                  Z( J4-2 ) = D + Z( J4-1 )
301                  TEMP = Z( J4+1 ) / Z( J4-2 )
302                  D = D*TEMP - TAU
303                  IF( D.LT.DTHRESH ) D = ZERO
304                  DMIN = MIN( DMIN, D )
305                  Z( J4 ) = Z( J4-1 )*TEMP
306                  EMIN = MIN( Z( J4 ), EMIN )
307 50            CONTINUE
308            ELSE
309               DO 60 J4 = 4*I0, 4*( N0-3 ), 4
310                  Z( J4-3 ) = D + Z( J4 )
311                  TEMP = Z( J4+2 ) / Z( J4-3 )
312                  D = D*TEMP - TAU
313                  IF( D.LT.DTHRESH ) D = ZERO
314                  DMIN = MIN( DMIN, D )
315                  Z( J4-1 ) = Z( J4 )*TEMP
316                  EMIN = MIN( Z( J4-1 ), EMIN )
317 60            CONTINUE
318            END IF
319*
320*     Unroll last two steps.
321*
322            DNM2 = D
323            DMIN2 = DMIN
324            J4 = 4*( N0-2 ) - PP
325            J4P2 = J4 + 2*PP - 1
326            Z( J4-2 ) = DNM2 + Z( J4P2 )
327            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
328            DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
329            DMIN = MIN( DMIN, DNM1 )
330*
331            DMIN1 = DMIN
332            J4 = J4 + 4
333            J4P2 = J4 + 2*PP - 1
334            Z( J4-2 ) = DNM1 + Z( J4P2 )
335            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
336            DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
337            DMIN = MIN( DMIN, DN )
338*
339         ELSE
340*
341*     Code for non IEEE arithmetic.
342*
343            IF( PP.EQ.0 ) THEN
344               DO 70 J4 = 4*I0, 4*( N0-3 ), 4
345                  Z( J4-2 ) = D + Z( J4-1 )
346                  IF( D.LT.ZERO ) THEN
347                     RETURN
348                  ELSE
349                     Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
350                     D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
351                  END IF
352                  IF( D.LT.DTHRESH) D = ZERO
353                  DMIN = MIN( DMIN, D )
354                  EMIN = MIN( EMIN, Z( J4 ) )
355 70            CONTINUE
356            ELSE
357               DO 80 J4 = 4*I0, 4*( N0-3 ), 4
358                  Z( J4-3 ) = D + Z( J4 )
359                  IF( D.LT.ZERO ) THEN
360                     RETURN
361                  ELSE
362                     Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
363                     D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
364                  END IF
365                  IF( D.LT.DTHRESH) D = ZERO
366                  DMIN = MIN( DMIN, D )
367                  EMIN = MIN( EMIN, Z( J4-1 ) )
368 80            CONTINUE
369            END IF
370*
371*     Unroll last two steps.
372*
373            DNM2 = D
374            DMIN2 = DMIN
375            J4 = 4*( N0-2 ) - PP
376            J4P2 = J4 + 2*PP - 1
377            Z( J4-2 ) = DNM2 + Z( J4P2 )
378            IF( DNM2.LT.ZERO ) THEN
379               RETURN
380            ELSE
381               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
382               DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
383            END IF
384            DMIN = MIN( DMIN, DNM1 )
385*
386            DMIN1 = DMIN
387            J4 = J4 + 4
388            J4P2 = J4 + 2*PP - 1
389            Z( J4-2 ) = DNM1 + Z( J4P2 )
390            IF( DNM1.LT.ZERO ) THEN
391               RETURN
392            ELSE
393               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
394               DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
395            END IF
396            DMIN = MIN( DMIN, DN )
397*
398         END IF
399      END IF
400*
401      Z( J4+2 ) = DN
402      Z( 4*N0-PP ) = EMIN
403      RETURN
404*
405*     End of DLASQ5
406*
407      END
408