LAPACK 3.12.1
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*> Download DLASQ5 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq5.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq5.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq5.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DLASQ5( 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* DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, SIGMA, EPS
26* ..
27* .. Array Arguments ..
28* DOUBLE PRECISION Z( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> DLASQ5 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 DOUBLE PRECISION 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 DOUBLE PRECISION
72*> This is the shift.
73*> \endverbatim
74*>
75*> \param[in] SIGMA
76*> \verbatim
77*> SIGMA is DOUBLE PRECISION
78*> This is the accumulated shift up to this step.
79*> \endverbatim
80*>
81*> \param[out] DMIN
82*> \verbatim
83*> DMIN is DOUBLE PRECISION
84*> Minimum value of d.
85*> \endverbatim
86*>
87*> \param[out] DMIN1
88*> \verbatim
89*> DMIN1 is DOUBLE PRECISION
90*> Minimum value of d, excluding D( N0 ).
91*> \endverbatim
92*>
93*> \param[out] DMIN2
94*> \verbatim
95*> DMIN2 is DOUBLE PRECISION
96*> Minimum value of d, excluding D( N0 ) and D( N0-1 ).
97*> \endverbatim
98*>
99*> \param[out] DN
100*> \verbatim
101*> DN is DOUBLE PRECISION
102*> d(N0), the last value of d.
103*> \endverbatim
104*>
105*> \param[out] DNM1
106*> \verbatim
107*> DNM1 is DOUBLE PRECISION
108*> d(N0-1).
109*> \endverbatim
110*>
111*> \param[out] DNM2
112*> \verbatim
113*> DNM2 is DOUBLE PRECISION
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 DOUBLE PRECISION
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 dlasq5( 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 DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
152 $ sigma, eps
153* ..
154* .. Array Arguments ..
155 DOUBLE PRECISION Z( * )
156* ..
157*
158* =====================================================================
159*
160* .. Parameter ..
161 DOUBLE PRECISION ZERO, HALF
162 PARAMETER ( ZERO = 0.0d0, half = 0.5 )
163* ..
164* .. Local Scalars ..
165 INTEGER J4, J4P2
166 DOUBLE PRECISION 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 ELSE
287* This is the version that sets d's to zero if they are small enough
288 j4 = 4*i0 + pp - 3
289 emin = z( j4+4 )
290 d = z( j4 ) - tau
291 dmin = d
292 dmin1 = -z( j4 )
293 IF( ieee ) THEN
294*
295* Code for IEEE arithmetic.
296*
297 IF( pp.EQ.0 ) THEN
298 DO 50 j4 = 4*i0, 4*( n0-3 ), 4
299 z( j4-2 ) = d + z( j4-1 )
300 temp = z( j4+1 ) / z( j4-2 )
301 d = d*temp - tau
302 IF( d.LT.dthresh ) d = zero
303 dmin = min( dmin, d )
304 z( j4 ) = z( j4-1 )*temp
305 emin = min( z( j4 ), emin )
306 50 CONTINUE
307 ELSE
308 DO 60 j4 = 4*i0, 4*( n0-3 ), 4
309 z( j4-3 ) = d + z( j4 )
310 temp = z( j4+2 ) / z( j4-3 )
311 d = d*temp - tau
312 IF( d.LT.dthresh ) d = zero
313 dmin = min( dmin, d )
314 z( j4-1 ) = z( j4 )*temp
315 emin = min( z( j4-1 ), emin )
316 60 CONTINUE
317 END IF
318*
319* Unroll last two steps.
320*
321 dnm2 = d
322 dmin2 = dmin
323 j4 = 4*( n0-2 ) - pp
324 j4p2 = j4 + 2*pp - 1
325 z( j4-2 ) = dnm2 + z( j4p2 )
326 z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
327 dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
328 dmin = min( dmin, dnm1 )
329*
330 dmin1 = dmin
331 j4 = j4 + 4
332 j4p2 = j4 + 2*pp - 1
333 z( j4-2 ) = dnm1 + z( j4p2 )
334 z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
335 dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
336 dmin = min( dmin, dn )
337*
338 ELSE
339*
340* Code for non IEEE arithmetic.
341*
342 IF( pp.EQ.0 ) THEN
343 DO 70 j4 = 4*i0, 4*( n0-3 ), 4
344 z( j4-2 ) = d + z( j4-1 )
345 IF( d.LT.zero ) THEN
346 RETURN
347 ELSE
348 z( j4 ) = z( j4+1 )*( z( j4-1 ) / z( j4-2 ) )
349 d = z( j4+1 )*( d / z( j4-2 ) ) - tau
350 END IF
351 IF( d.LT.dthresh) d = zero
352 dmin = min( dmin, d )
353 emin = min( emin, z( j4 ) )
354 70 CONTINUE
355 ELSE
356 DO 80 j4 = 4*i0, 4*( n0-3 ), 4
357 z( j4-3 ) = d + z( j4 )
358 IF( d.LT.zero ) THEN
359 RETURN
360 ELSE
361 z( j4-1 ) = z( j4+2 )*( z( j4 ) / z( j4-3 ) )
362 d = z( j4+2 )*( d / z( j4-3 ) ) - tau
363 END IF
364 IF( d.LT.dthresh) d = zero
365 dmin = min( dmin, d )
366 emin = min( emin, z( j4-1 ) )
367 80 CONTINUE
368 END IF
369*
370* Unroll last two steps.
371*
372 dnm2 = d
373 dmin2 = dmin
374 j4 = 4*( n0-2 ) - pp
375 j4p2 = j4 + 2*pp - 1
376 z( j4-2 ) = dnm2 + z( j4p2 )
377 IF( dnm2.LT.zero ) THEN
378 RETURN
379 ELSE
380 z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
381 dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
382 END IF
383 dmin = min( dmin, dnm1 )
384*
385 dmin1 = dmin
386 j4 = j4 + 4
387 j4p2 = j4 + 2*pp - 1
388 z( j4-2 ) = dnm1 + z( j4p2 )
389 IF( dnm1.LT.zero ) THEN
390 RETURN
391 ELSE
392 z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
393 dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
394 END IF
395 dmin = min( dmin, dn )
396*
397 END IF
398 END IF
399*
400 z( j4+2 ) = dn
401 z( 4*n0-pp ) = emin
402 RETURN
403*
404* End of DLASQ5
405*
406 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:143