LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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. </b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download SLASQ5 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq5.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq5.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq5.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
20* DNM1, DNM2, IEEE, EPS )
21*
22* .. Scalar Arguments ..
23* LOGICAL IEEE
24* INTEGER I0, N0, PP
25* REAL EPS, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, SIGMA, TAU
26* ..
27* .. Array Arguments ..
28* REAL Z( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> SLASQ5 computes one dqds transform in ping-pong form, one
38*> version for IEEE machines another for non IEEE machines.
39*> \endverbatim
40*
41* Arguments:
42* ==========
43*
44*> \param[in] I0
45*> \verbatim
46*> I0 is INTEGER
47*> First index.
48*> \endverbatim
49*>
50*> \param[in] N0
51*> \verbatim
52*> N0 is INTEGER
53*> Last index.
54*> \endverbatim
55*>
56*> \param[in] Z
57*> \verbatim
58*> Z is REAL array, dimension ( 4*N )
59*> Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
60*> an extra argument.
61*> \endverbatim
62*>
63*> \param[in] PP
64*> \verbatim
65*> PP is INTEGER
66*> PP=0 for ping, PP=1 for pong.
67*> \endverbatim
68*>
69*> \param[in] TAU
70*> \verbatim
71*> TAU is REAL
72*> This is the shift.
73*> \endverbatim
74*>
75*> \param[in] SIGMA
76*> \verbatim
77*> SIGMA is REAL
78*> This is the accumulated shift up to this step.
79*> \endverbatim
80*>
81*> \param[out] DMIN
82*> \verbatim
83*> DMIN is REAL
84*> Minimum value of d.
85*> \endverbatim
86*>
87*> \param[out] DMIN1
88*> \verbatim
89*> DMIN1 is REAL
90*> Minimum value of d, excluding D( N0 ).
91*> \endverbatim
92*>
93*> \param[out] DMIN2
94*> \verbatim
95*> DMIN2 is REAL
96*> Minimum value of d, excluding D( N0 ) and D( N0-1 ).
97*> \endverbatim
98*>
99*> \param[out] DN
100*> \verbatim
101*> DN is REAL
102*> d(N0), the last value of d.
103*> \endverbatim
104*>
105*> \param[out] DNM1
106*> \verbatim
107*> DNM1 is REAL
108*> d(N0-1).
109*> \endverbatim
110*>
111*> \param[out] DNM2
112*> \verbatim
113*> DNM2 is REAL
114*> d(N0-2).
115*> \endverbatim
116*>
117*> \param[in] IEEE
118*> \verbatim
119*> IEEE is LOGICAL
120*> Flag for IEEE or non IEEE arithmetic.
121*> \endverbatim
122*>
123*> \param[in] EPS
124*> \verbatim
125*> EPS is REAL
126*> This is the value of epsilon used.
127*> \endverbatim
128*
129* Authors:
130* ========
131*
132*> \author Univ. of Tennessee
133*> \author Univ. of California Berkeley
134*> \author Univ. of Colorado Denver
135*> \author NAG Ltd.
136*
137*> \ingroup lasq5
138*
139* =====================================================================
140 SUBROUTINE slasq5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1,
141 $ DMIN2,
142 $ DN, DNM1, DNM2, IEEE, EPS )
143*
144* -- LAPACK computational routine --
145* -- LAPACK is a software package provided by Univ. of Tennessee, --
146* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
147*
148* .. Scalar Arguments ..
149 LOGICAL IEEE
150 INTEGER I0, N0, PP
151 REAL DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
152 $ sigma, eps
153* ..
154* .. Array Arguments ..
155 REAL Z( * )
156* ..
157*
158* =====================================================================
159*
160* .. Parameter ..
161 REAL ZERO, HALF
162 PARAMETER ( ZERO = 0.0e0, half = 0.5 )
163* ..
164* .. Local Scalars ..
165 INTEGER J4, J4P2
166 REAL D, EMIN, TEMP, DTHRESH
167* ..
168* .. Intrinsic Functions ..
169 INTRINSIC min
170* ..
171* .. Executable Statements ..
172*
173 IF( ( n0-i0-1 ).LE.0 )
174 $ RETURN
175*
176 dthresh = eps*(sigma+tau)
177 IF( tau.LT.dthresh*half ) tau = zero
178 IF( tau.NE.zero ) THEN
179 j4 = 4*i0 + pp - 3
180 emin = z( j4+4 )
181 d = z( j4 ) - tau
182 dmin = d
183 dmin1 = -z( j4 )
184*
185 IF( ieee ) THEN
186*
187* Code for IEEE arithmetic.
188*
189 IF( pp.EQ.0 ) THEN
190 DO 10 j4 = 4*i0, 4*( n0-3 ), 4
191 z( j4-2 ) = d + z( j4-1 )
192 temp = z( j4+1 ) / z( j4-2 )
193 d = d*temp - tau
194 dmin = min( dmin, d )
195 z( j4 ) = z( j4-1 )*temp
196 emin = min( z( j4 ), emin )
197 10 CONTINUE
198 ELSE
199 DO 20 j4 = 4*i0, 4*( n0-3 ), 4
200 z( j4-3 ) = d + z( j4 )
201 temp = z( j4+2 ) / z( j4-3 )
202 d = d*temp - tau
203 dmin = min( dmin, d )
204 z( j4-1 ) = z( j4 )*temp
205 emin = min( z( j4-1 ), emin )
206 20 CONTINUE
207 END IF
208*
209* Unroll last two steps.
210*
211 dnm2 = d
212 dmin2 = dmin
213 j4 = 4*( n0-2 ) - pp
214 j4p2 = j4 + 2*pp - 1
215 z( j4-2 ) = dnm2 + z( j4p2 )
216 z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
217 dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
218 dmin = min( dmin, dnm1 )
219*
220 dmin1 = dmin
221 j4 = j4 + 4
222 j4p2 = j4 + 2*pp - 1
223 z( j4-2 ) = dnm1 + z( j4p2 )
224 z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
225 dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
226 dmin = min( dmin, dn )
227*
228 ELSE
229*
230* Code for non IEEE arithmetic.
231*
232 IF( pp.EQ.0 ) THEN
233 DO 30 j4 = 4*i0, 4*( n0-3 ), 4
234 z( j4-2 ) = d + z( j4-1 )
235 IF( d.LT.zero ) THEN
236 RETURN
237 ELSE
238 z( j4 ) = z( j4+1 )*( z( j4-1 ) / z( j4-2 ) )
239 d = z( j4+1 )*( d / z( j4-2 ) ) - tau
240 END IF
241 dmin = min( dmin, d )
242 emin = min( emin, z( j4 ) )
243 30 CONTINUE
244 ELSE
245 DO 40 j4 = 4*i0, 4*( n0-3 ), 4
246 z( j4-3 ) = d + z( j4 )
247 IF( d.LT.zero ) THEN
248 RETURN
249 ELSE
250 z( j4-1 ) = z( j4+2 )*( z( j4 ) / z( j4-3 ) )
251 d = z( j4+2 )*( d / z( j4-3 ) ) - tau
252 END IF
253 dmin = min( dmin, d )
254 emin = min( emin, z( j4-1 ) )
255 40 CONTINUE
256 END IF
257*
258* Unroll last two steps.
259*
260 dnm2 = d
261 dmin2 = dmin
262 j4 = 4*( n0-2 ) - pp
263 j4p2 = j4 + 2*pp - 1
264 z( j4-2 ) = dnm2 + z( j4p2 )
265 IF( dnm2.LT.zero ) THEN
266 RETURN
267 ELSE
268 z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
269 dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
270 END IF
271 dmin = min( dmin, dnm1 )
272*
273 dmin1 = dmin
274 j4 = j4 + 4
275 j4p2 = j4 + 2*pp - 1
276 z( j4-2 ) = dnm1 + z( j4p2 )
277 IF( dnm1.LT.zero ) THEN
278 RETURN
279 ELSE
280 z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
281 dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
282 END IF
283 dmin = min( dmin, dn )
284*
285 END IF
286*
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*
400 END IF
401 z( j4+2 ) = dn
402 z( 4*n0-pp ) = emin
403 RETURN
404*
405* End of SLASQ5
406*
407 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:143