LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 lasq5
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
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:144