LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slasq5.f
Go to the documentation of this file.
1 *> \brief \b SLASQ5 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 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, DMIN, DMIN1, DMIN2, DN,
22 * DNM1, DNM2, IEEE )
23 *
24 * .. Scalar Arguments ..
25 * LOGICAL IEEE
26 * INTEGER I0, N0, PP
27 * REAL DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, 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[out] DMIN
78 *> \verbatim
79 *> DMIN is REAL
80 *> Minimum value of d.
81 *> \endverbatim
82 *>
83 *> \param[out] DMIN1
84 *> \verbatim
85 *> DMIN1 is REAL
86 *> Minimum value of d, excluding D( N0 ).
87 *> \endverbatim
88 *>
89 *> \param[out] DMIN2
90 *> \verbatim
91 *> DMIN2 is REAL
92 *> Minimum value of d, excluding D( N0 ) and D( N0-1 ).
93 *> \endverbatim
94 *>
95 *> \param[out] DN
96 *> \verbatim
97 *> DN is REAL
98 *> d(N0), the last value of d.
99 *> \endverbatim
100 *>
101 *> \param[out] DNM1
102 *> \verbatim
103 *> DNM1 is REAL
104 *> d(N0-1).
105 *> \endverbatim
106 *>
107 *> \param[out] DNM2
108 *> \verbatim
109 *> DNM2 is REAL
110 *> d(N0-2).
111 *> \endverbatim
112 *>
113 *> \param[in] IEEE
114 *> \verbatim
115 *> IEEE is LOGICAL
116 *> Flag for IEEE or non IEEE arithmetic.
117 *> \endverbatim
118 *
119 * Authors:
120 * ========
121 *
122 *> \author Univ. of Tennessee
123 *> \author Univ. of California Berkeley
124 *> \author Univ. of Colorado Denver
125 *> \author NAG Ltd.
126 *
127 *> \date September 2012
128 *
129 *> \ingroup auxOTHERcomputational
130 *
131 * =====================================================================
132  SUBROUTINE slasq5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
133  $ dn, dnm1, dnm2, ieee, eps )
134 *
135 * -- LAPACK computational routine (version 3.4.2) --
136 * -- LAPACK is a software package provided by Univ. of Tennessee, --
137 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138 * September 2012
139 *
140 * .. Scalar Arguments ..
141  LOGICAL ieee
142  INTEGER i0, n0, pp
143  REAL dmin, dmin1, dmin2, dn, dnm1, dnm2, tau,
144  $ sigma, eps
145 * ..
146 * .. Array Arguments ..
147  REAL z( * )
148 * ..
149 *
150 * =====================================================================
151 *
152 * .. Parameter ..
153  REAL zero, half
154  parameter( zero = 0.0e0, half = 0.5 )
155 * ..
156 * .. Local Scalars ..
157  INTEGER j4, j4p2
158  REAL d, emin, temp, dthresh
159 * ..
160 * .. Intrinsic Functions ..
161  INTRINSIC min
162 * ..
163 * .. Executable Statements ..
164 *
165  IF( ( n0-i0-1 ).LE.0 )
166  $ return
167 *
168  dthresh = eps*(sigma+tau)
169  IF( tau.LT.dthresh*half ) tau = zero
170  IF( tau.NE.zero ) THEN
171  j4 = 4*i0 + pp - 3
172  emin = z( j4+4 )
173  d = z( j4 ) - tau
174  dmin = d
175  dmin1 = -z( j4 )
176 *
177  IF( ieee ) THEN
178 *
179 * Code for IEEE arithmetic.
180 *
181  IF( pp.EQ.0 ) THEN
182  DO 10 j4 = 4*i0, 4*( n0-3 ), 4
183  z( j4-2 ) = d + z( j4-1 )
184  temp = z( j4+1 ) / z( j4-2 )
185  d = d*temp - tau
186  dmin = min( dmin, d )
187  z( j4 ) = z( j4-1 )*temp
188  emin = min( z( j4 ), emin )
189  10 continue
190  ELSE
191  DO 20 j4 = 4*i0, 4*( n0-3 ), 4
192  z( j4-3 ) = d + z( j4 )
193  temp = z( j4+2 ) / z( j4-3 )
194  d = d*temp - tau
195  dmin = min( dmin, d )
196  z( j4-1 ) = z( j4 )*temp
197  emin = min( z( j4-1 ), emin )
198  20 continue
199  END IF
200 *
201 * Unroll last two steps.
202 *
203  dnm2 = d
204  dmin2 = dmin
205  j4 = 4*( n0-2 ) - pp
206  j4p2 = j4 + 2*pp - 1
207  z( j4-2 ) = dnm2 + z( j4p2 )
208  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
209  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
210  dmin = min( dmin, dnm1 )
211 *
212  dmin1 = dmin
213  j4 = j4 + 4
214  j4p2 = j4 + 2*pp - 1
215  z( j4-2 ) = dnm1 + z( j4p2 )
216  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
217  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
218  dmin = min( dmin, dn )
219 *
220  ELSE
221 *
222 * Code for non IEEE arithmetic.
223 *
224  IF( pp.EQ.0 ) THEN
225  DO 30 j4 = 4*i0, 4*( n0-3 ), 4
226  z( j4-2 ) = d + z( j4-1 )
227  IF( d.LT.zero ) THEN
228  return
229  ELSE
230  z( j4 ) = z( j4+1 )*( z( j4-1 ) / z( j4-2 ) )
231  d = z( j4+1 )*( d / z( j4-2 ) ) - tau
232  END IF
233  dmin = min( dmin, d )
234  emin = min( emin, z( j4 ) )
235  30 continue
236  ELSE
237  DO 40 j4 = 4*i0, 4*( n0-3 ), 4
238  z( j4-3 ) = d + z( j4 )
239  IF( d.LT.zero ) THEN
240  return
241  ELSE
242  z( j4-1 ) = z( j4+2 )*( z( j4 ) / z( j4-3 ) )
243  d = z( j4+2 )*( d / z( j4-3 ) ) - tau
244  END IF
245  dmin = min( dmin, d )
246  emin = min( emin, z( j4-1 ) )
247  40 continue
248  END IF
249 *
250 * Unroll last two steps.
251 *
252  dnm2 = d
253  dmin2 = dmin
254  j4 = 4*( n0-2 ) - pp
255  j4p2 = j4 + 2*pp - 1
256  z( j4-2 ) = dnm2 + z( j4p2 )
257  IF( dnm2.LT.zero ) THEN
258  return
259  ELSE
260  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
261  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
262  END IF
263  dmin = min( dmin, dnm1 )
264 *
265  dmin1 = dmin
266  j4 = j4 + 4
267  j4p2 = j4 + 2*pp - 1
268  z( j4-2 ) = dnm1 + z( j4p2 )
269  IF( dnm1.LT.zero ) THEN
270  return
271  ELSE
272  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
273  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
274  END IF
275  dmin = min( dmin, dn )
276 *
277  END IF
278 *
279  ELSE
280 * This is the version that sets d's to zero if they are small enough
281  j4 = 4*i0 + pp - 3
282  emin = z( j4+4 )
283  d = z( j4 ) - tau
284  dmin = d
285  dmin1 = -z( j4 )
286  IF( ieee ) THEN
287 *
288 * Code for IEEE arithmetic.
289 *
290  IF( pp.EQ.0 ) THEN
291  DO 50 j4 = 4*i0, 4*( n0-3 ), 4
292  z( j4-2 ) = d + z( j4-1 )
293  temp = z( j4+1 ) / z( j4-2 )
294  d = d*temp - tau
295  IF( d.LT.dthresh ) d = zero
296  dmin = min( dmin, d )
297  z( j4 ) = z( j4-1 )*temp
298  emin = min( z( j4 ), emin )
299  50 continue
300  ELSE
301  DO 60 j4 = 4*i0, 4*( n0-3 ), 4
302  z( j4-3 ) = d + z( j4 )
303  temp = z( j4+2 ) / z( j4-3 )
304  d = d*temp - tau
305  IF( d.LT.dthresh ) d = zero
306  dmin = min( dmin, d )
307  z( j4-1 ) = z( j4 )*temp
308  emin = min( z( j4-1 ), emin )
309  60 continue
310  END IF
311 *
312 * Unroll last two steps.
313 *
314  dnm2 = d
315  dmin2 = dmin
316  j4 = 4*( n0-2 ) - pp
317  j4p2 = j4 + 2*pp - 1
318  z( j4-2 ) = dnm2 + z( j4p2 )
319  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
320  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
321  dmin = min( dmin, dnm1 )
322 *
323  dmin1 = dmin
324  j4 = j4 + 4
325  j4p2 = j4 + 2*pp - 1
326  z( j4-2 ) = dnm1 + z( j4p2 )
327  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
328  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
329  dmin = min( dmin, dn )
330 *
331  ELSE
332 *
333 * Code for non IEEE arithmetic.
334 *
335  IF( pp.EQ.0 ) THEN
336  DO 70 j4 = 4*i0, 4*( n0-3 ), 4
337  z( j4-2 ) = d + z( j4-1 )
338  IF( d.LT.zero ) THEN
339  return
340  ELSE
341  z( j4 ) = z( j4+1 )*( z( j4-1 ) / z( j4-2 ) )
342  d = z( j4+1 )*( d / z( j4-2 ) ) - tau
343  END IF
344  IF( d.LT.dthresh ) d = zero
345  dmin = min( dmin, d )
346  emin = min( emin, z( j4 ) )
347  70 continue
348  ELSE
349  DO 80 j4 = 4*i0, 4*( n0-3 ), 4
350  z( j4-3 ) = d + z( j4 )
351  IF( d.LT.zero ) THEN
352  return
353  ELSE
354  z( j4-1 ) = z( j4+2 )*( z( j4 ) / z( j4-3 ) )
355  d = z( j4+2 )*( d / z( j4-3 ) ) - tau
356  END IF
357  IF( d.LT.dthresh ) d = zero
358  dmin = min( dmin, d )
359  emin = min( emin, z( j4-1 ) )
360  80 continue
361  END IF
362 *
363 * Unroll last two steps.
364 *
365  dnm2 = d
366  dmin2 = dmin
367  j4 = 4*( n0-2 ) - pp
368  j4p2 = j4 + 2*pp - 1
369  z( j4-2 ) = dnm2 + z( j4p2 )
370  IF( dnm2.LT.zero ) THEN
371  return
372  ELSE
373  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
374  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
375  END IF
376  dmin = min( dmin, dnm1 )
377 *
378  dmin1 = dmin
379  j4 = j4 + 4
380  j4p2 = j4 + 2*pp - 1
381  z( j4-2 ) = dnm1 + z( j4p2 )
382  IF( dnm1.LT.zero ) THEN
383  return
384  ELSE
385  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
386  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
387  END IF
388  dmin = min( dmin, dn )
389 *
390  END IF
391 *
392  END IF
393  z( j4+2 ) = dn
394  z( 4*n0-pp ) = emin
395  return
396 *
397 * End of SLASQ5
398 *
399  END