LAPACK 3.12.1
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 140 of file dlasq5.f.

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