LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlasq5()

subroutine dlasq5 ( integer  i0,
integer  n0,
double precision, dimension( * )  z,
integer  pp,
double precision  tau,
double precision  sigma,
double precision  dmin,
double precision  dmin1,
double precision  dmin2,
double precision  dn,
double precision  dnm1,
double precision  dnm2,
logical  ieee,
double precision  eps 
)

DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.

Download DLASQ5 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DLASQ5 computes one dqds transform in ping-pong form, one
 version for IEEE machines another for non IEEE machines.
Parameters
[in]I0
          I0 is INTEGER
        First index.
[in]N0
          N0 is INTEGER
        Last index.
[in]Z
          Z is DOUBLE PRECISION array, dimension ( 4*N )
        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
        an extra argument.
[in]PP
          PP is INTEGER
        PP=0 for ping, PP=1 for pong.
[in]TAU
          TAU is DOUBLE PRECISION
        This is the shift.
[in]SIGMA
          SIGMA is DOUBLE PRECISION
        This is the accumulated shift up to this step.
[out]DMIN
          DMIN is DOUBLE PRECISION
        Minimum value of d.
[out]DMIN1
          DMIN1 is DOUBLE PRECISION
        Minimum value of d, excluding D( N0 ).
[out]DMIN2
          DMIN2 is DOUBLE PRECISION
        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
[out]DN
          DN is DOUBLE PRECISION
        d(N0), the last value of d.
[out]DNM1
          DNM1 is DOUBLE PRECISION
        d(N0-1).
[out]DNM2
          DNM2 is DOUBLE PRECISION
        d(N0-2).
[in]IEEE
          IEEE is LOGICAL
        Flag for IEEE or non IEEE arithmetic.
[in]EPS
          EPS is DOUBLE PRECISION
        This is the value of epsilon used.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 142 of file dlasq5.f.

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*
Here is the caller graph for this function: