LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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, SIGMA, DMIN, DMIN1, DMIN2, DN,
22 * DNM1, DNM2, IEEE, EPS )
23 *
24 * .. Scalar Arguments ..
25 * LOGICAL IEEE
26 * INTEGER I0, N0, PP
27 * REAL EPS, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, SIGMA, 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[in] SIGMA
78 *> SIGMA is REAL
79 *> This is the accumulated shift up to this step.
80 *> \endverbatim
81 *>
82 *> \param[out] DMIN
83 *> \verbatim
84 *> DMIN is REAL
85 *> Minimum value of d.
86 *> \endverbatim
87 *>
88 *> \param[out] DMIN1
89 *> \verbatim
90 *> DMIN1 is REAL
91 *> Minimum value of d, excluding D( N0 ).
92 *> \endverbatim
93 *>
94 *> \param[out] DMIN2
95 *> \verbatim
96 *> DMIN2 is REAL
97 *> Minimum value of d, excluding D( N0 ) and D( N0-1 ).
98 *> \endverbatim
99 *>
100 *> \param[out] DN
101 *> \verbatim
102 *> DN is REAL
103 *> d(N0), the last value of d.
104 *> \endverbatim
105 *>
106 *> \param[out] DNM1
107 *> \verbatim
108 *> DNM1 is REAL
109 *> d(N0-1).
110 *> \endverbatim
111 *>
112 *> \param[out] DNM2
113 *> \verbatim
114 *> DNM2 is REAL
115 *> d(N0-2).
116 *> \endverbatim
117 *>
118 *> \param[in] IEEE
119 *> \verbatim
120 *> IEEE is LOGICAL
121 *> Flag for IEEE or non IEEE arithmetic.
122 *> \endverbatim
123 *>
124 *> \param[in] EPS
125 *> \verbatim
126 *> EPS is REAL
127 *> This is the value of epsilon used.
128 *> \endverbatim
129 *
130 * Authors:
131 * ========
132 *
133 *> \author Univ. of Tennessee
134 *> \author Univ. of California Berkeley
135 *> \author Univ. of Colorado Denver
136 *> \author NAG Ltd.
137 *
138 *> \date November 2015
139 *
140 *> \ingroup auxOTHERcomputational
141 *
142 * =====================================================================
143  SUBROUTINE slasq5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
144  $ dn, dnm1, dnm2, ieee, eps )
145 *
146 * -- LAPACK computational routine (version 3.6.0) --
147 * -- LAPACK is a software package provided by Univ. of Tennessee, --
148 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
149 * November 2015
150 *
151 * .. Scalar Arguments ..
152  LOGICAL IEEE
153  INTEGER I0, N0, PP
154  REAL DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
155  $ sigma, eps
156 * ..
157 * .. Array Arguments ..
158  REAL Z( * )
159 * ..
160 *
161 * =====================================================================
162 *
163 * .. Parameter ..
164  REAL ZERO, HALF
165  parameter ( zero = 0.0e0, half = 0.5 )
166 * ..
167 * .. Local Scalars ..
168  INTEGER J4, J4P2
169  REAL D, EMIN, TEMP, DTHRESH
170 * ..
171 * .. Intrinsic Functions ..
172  INTRINSIC min
173 * ..
174 * .. Executable Statements ..
175 *
176  IF( ( n0-i0-1 ).LE.0 )
177  $ RETURN
178 *
179  dthresh = eps*(sigma+tau)
180  IF( tau.LT.dthresh*half ) tau = zero
181  IF( tau.NE.zero ) THEN
182  j4 = 4*i0 + pp - 3
183  emin = z( j4+4 )
184  d = z( j4 ) - tau
185  dmin = d
186  dmin1 = -z( j4 )
187 *
188  IF( ieee ) THEN
189 *
190 * Code for IEEE arithmetic.
191 *
192  IF( pp.EQ.0 ) THEN
193  DO 10 j4 = 4*i0, 4*( n0-3 ), 4
194  z( j4-2 ) = d + z( j4-1 )
195  temp = z( j4+1 ) / z( j4-2 )
196  d = d*temp - tau
197  dmin = min( dmin, d )
198  z( j4 ) = z( j4-1 )*temp
199  emin = min( z( j4 ), emin )
200  10 CONTINUE
201  ELSE
202  DO 20 j4 = 4*i0, 4*( n0-3 ), 4
203  z( j4-3 ) = d + z( j4 )
204  temp = z( j4+2 ) / z( j4-3 )
205  d = d*temp - tau
206  dmin = min( dmin, d )
207  z( j4-1 ) = z( j4 )*temp
208  emin = min( z( j4-1 ), emin )
209  20 CONTINUE
210  END IF
211 *
212 * Unroll last two steps.
213 *
214  dnm2 = d
215  dmin2 = dmin
216  j4 = 4*( n0-2 ) - pp
217  j4p2 = j4 + 2*pp - 1
218  z( j4-2 ) = dnm2 + z( j4p2 )
219  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
220  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
221  dmin = min( dmin, dnm1 )
222 *
223  dmin1 = dmin
224  j4 = j4 + 4
225  j4p2 = j4 + 2*pp - 1
226  z( j4-2 ) = dnm1 + z( j4p2 )
227  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
228  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
229  dmin = min( dmin, dn )
230 *
231  ELSE
232 *
233 * Code for non IEEE arithmetic.
234 *
235  IF( pp.EQ.0 ) THEN
236  DO 30 j4 = 4*i0, 4*( n0-3 ), 4
237  z( j4-2 ) = d + z( j4-1 )
238  IF( d.LT.zero ) THEN
239  RETURN
240  ELSE
241  z( j4 ) = z( j4+1 )*( z( j4-1 ) / z( j4-2 ) )
242  d = z( j4+1 )*( d / z( j4-2 ) ) - tau
243  END IF
244  dmin = min( dmin, d )
245  emin = min( emin, z( j4 ) )
246  30 CONTINUE
247  ELSE
248  DO 40 j4 = 4*i0, 4*( n0-3 ), 4
249  z( j4-3 ) = d + z( j4 )
250  IF( d.LT.zero ) THEN
251  RETURN
252  ELSE
253  z( j4-1 ) = z( j4+2 )*( z( j4 ) / z( j4-3 ) )
254  d = z( j4+2 )*( d / z( j4-3 ) ) - tau
255  END IF
256  dmin = min( dmin, d )
257  emin = min( emin, z( j4-1 ) )
258  40 CONTINUE
259  END IF
260 *
261 * Unroll last two steps.
262 *
263  dnm2 = d
264  dmin2 = dmin
265  j4 = 4*( n0-2 ) - pp
266  j4p2 = j4 + 2*pp - 1
267  z( j4-2 ) = dnm2 + z( j4p2 )
268  IF( dnm2.LT.zero ) THEN
269  RETURN
270  ELSE
271  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
272  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
273  END IF
274  dmin = min( dmin, dnm1 )
275 *
276  dmin1 = dmin
277  j4 = j4 + 4
278  j4p2 = j4 + 2*pp - 1
279  z( j4-2 ) = dnm1 + z( j4p2 )
280  IF( dnm1.LT.zero ) THEN
281  RETURN
282  ELSE
283  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
284  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
285  END IF
286  dmin = min( dmin, dn )
287 *
288  END IF
289 *
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 *
403  END IF
404  z( j4+2 ) = dn
405  z( 4*n0-pp ) = emin
406  RETURN
407 *
408 * End of SLASQ5
409 *
410  END
slasq5
subroutine slasq5(I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, IEEE, EPS)
SLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
Definition: slasq5.f:145