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*> \date September 2012
140*
141*> \ingroup auxOTHERcomputational
142*
143*  =====================================================================
144      SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
145     $                   DN, DNM1, DNM2, IEEE, EPS )
146*
147*  -- LAPACK computational routine (version 3.4.2) --
148*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
149*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
150*     September 2012
151*
152*     .. Scalar Arguments ..
153      LOGICAL            IEEE
154      INTEGER            I0, N0, PP
155      DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
156     $                   SIGMA, EPS
157*     ..
158*     .. Array Arguments ..
159      DOUBLE PRECISION   Z( * )
160*     ..
161*
162*  =====================================================================
163*
164*     .. Parameter ..
165      DOUBLE PRECISION   ZERO, HALF
166      PARAMETER          ( ZERO = 0.0D0, HALF = 0.5 )
167*     ..
168*     .. Local Scalars ..
169      INTEGER            J4, J4P2
170      DOUBLE PRECISION   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      ELSE
291*     This is the version that sets d's to zero if they are small enough
292         J4 = 4*I0 + PP - 3
293         EMIN = Z( J4+4 )
294         D = Z( J4 ) - TAU
295         DMIN = D
296         DMIN1 = -Z( J4 )
297         IF( IEEE ) THEN
298*
299*     Code for IEEE arithmetic.
300*
301            IF( PP.EQ.0 ) THEN
302               DO 50 J4 = 4*I0, 4*( N0-3 ), 4
303                  Z( J4-2 ) = D + Z( J4-1 )
304                  TEMP = Z( J4+1 ) / Z( J4-2 )
305                  D = D*TEMP - TAU
306                  IF( D.LT.DTHRESH ) D = ZERO
307                  DMIN = MIN( DMIN, D )
308                  Z( J4 ) = Z( J4-1 )*TEMP
309                  EMIN = MIN( Z( J4 ), EMIN )
310 50            CONTINUE
311            ELSE
312               DO 60 J4 = 4*I0, 4*( N0-3 ), 4
313                  Z( J4-3 ) = D + Z( J4 )
314                  TEMP = Z( J4+2 ) / Z( J4-3 )
315                  D = D*TEMP - TAU
316                  IF( D.LT.DTHRESH ) D = ZERO
317                  DMIN = MIN( DMIN, D )
318                  Z( J4-1 ) = Z( J4 )*TEMP
319                  EMIN = MIN( Z( J4-1 ), EMIN )
320 60            CONTINUE
321            END IF
322*
323*     Unroll last two steps.
324*
325            DNM2 = D
326            DMIN2 = DMIN
327            J4 = 4*( N0-2 ) - PP
328            J4P2 = J4 + 2*PP - 1
329            Z( J4-2 ) = DNM2 + Z( J4P2 )
330            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
331            DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
332            DMIN = MIN( DMIN, DNM1 )
333*
334            DMIN1 = DMIN
335            J4 = J4 + 4
336            J4P2 = J4 + 2*PP - 1
337            Z( J4-2 ) = DNM1 + Z( J4P2 )
338            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
339            DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
340            DMIN = MIN( DMIN, DN )
341*
342         ELSE
343*
344*     Code for non IEEE arithmetic.
345*
346            IF( PP.EQ.0 ) THEN
347               DO 70 J4 = 4*I0, 4*( N0-3 ), 4
348                  Z( J4-2 ) = D + Z( J4-1 )
349                  IF( D.LT.ZERO ) THEN
350                     RETURN
351                  ELSE
352                     Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
353                     D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
354                  END IF
355                  IF( D.LT.DTHRESH) D = ZERO
356                  DMIN = MIN( DMIN, D )
357                  EMIN = MIN( EMIN, Z( J4 ) )
358 70            CONTINUE
359            ELSE
360               DO 80 J4 = 4*I0, 4*( N0-3 ), 4
361                  Z( J4-3 ) = D + Z( J4 )
362                  IF( D.LT.ZERO ) THEN
363                     RETURN
364                  ELSE
365                     Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
366                     D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
367                  END IF
368                  IF( D.LT.DTHRESH) D = ZERO
369                  DMIN = MIN( DMIN, D )
370                  EMIN = MIN( EMIN, Z( J4-1 ) )
371 80            CONTINUE
372            END IF
373*
374*     Unroll last two steps.
375*
376            DNM2 = D
377            DMIN2 = DMIN
378            J4 = 4*( N0-2 ) - PP
379            J4P2 = J4 + 2*PP - 1
380            Z( J4-2 ) = DNM2 + Z( J4P2 )
381            IF( DNM2.LT.ZERO ) THEN
382               RETURN
383            ELSE
384               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
385               DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
386            END IF
387            DMIN = MIN( DMIN, DNM1 )
388*
389            DMIN1 = DMIN
390            J4 = J4 + 4
391            J4P2 = J4 + 2*PP - 1
392            Z( J4-2 ) = DNM1 + Z( J4P2 )
393            IF( DNM1.LT.ZERO ) THEN
394               RETURN
395            ELSE
396               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
397               DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
398            END IF
399            DMIN = MIN( DMIN, DN )
400*
401         END IF
402      END IF
403*
404      Z( J4+2 ) = DN
405      Z( 4*N0-PP ) = EMIN
406      RETURN
407*
408*     End of DLASQ5
409*
410      END
411c $Id$
412