LAPACK  3.6.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 *> \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
dlasq5
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:146