LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
dlasq5.f
Go to the documentation of this file.
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
subroutine dlasq5(I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, IEEE, EPS)
DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
Definition: dlasq5.f:144