LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
September 2012

Definition at line 146 of file dlasq5.f.

146 *
147 * -- LAPACK computational routine (version 3.4.2) --
148 * -- LAPACK is a software package provided by Univ. of Tennessee, --
149 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
150 * September 2012
151 *
152 * .. Scalar Arguments ..
153  LOGICAL ieee
154  INTEGER i0, n0, pp
155  DOUBLE PRECISION dmin, dmin1, dmin2, dn, dnm1, dnm2, tau,
156  $ sigma, eps
157 * ..
158 * .. Array Arguments ..
159  DOUBLE PRECISION z( * )
160 * ..
161 *
162 * =====================================================================
163 *
164 * .. Parameter ..
165  DOUBLE PRECISION zero, half
166  parameter ( zero = 0.0d0, half = 0.5 )
167 * ..
168 * .. Local Scalars ..
169  INTEGER j4, j4p2
170  DOUBLE PRECISION d, emin, temp, dthresh
171 * ..
172 * .. Intrinsic Functions ..
173  INTRINSIC min
174 * ..
175 * .. Executable Statements ..
176 *
177  IF( ( n0-i0-1 ).LE.0 )
178  $ RETURN
179 *
180  dthresh = eps*(sigma+tau)
181  IF( tau.LT.dthresh*half ) tau = zero
182  IF( tau.NE.zero ) THEN
183  j4 = 4*i0 + pp - 3
184  emin = z( j4+4 )
185  d = z( j4 ) - tau
186  dmin = d
187  dmin1 = -z( j4 )
188 *
189  IF( ieee ) THEN
190 *
191 * Code for IEEE arithmetic.
192 *
193  IF( pp.EQ.0 ) THEN
194  DO 10 j4 = 4*i0, 4*( n0-3 ), 4
195  z( j4-2 ) = d + z( j4-1 )
196  temp = z( j4+1 ) / z( j4-2 )
197  d = d*temp - tau
198  dmin = min( dmin, d )
199  z( j4 ) = z( j4-1 )*temp
200  emin = min( z( j4 ), emin )
201  10 CONTINUE
202  ELSE
203  DO 20 j4 = 4*i0, 4*( n0-3 ), 4
204  z( j4-3 ) = d + z( j4 )
205  temp = z( j4+2 ) / z( j4-3 )
206  d = d*temp - tau
207  dmin = min( dmin, d )
208  z( j4-1 ) = z( j4 )*temp
209  emin = min( z( j4-1 ), emin )
210  20 CONTINUE
211  END IF
212 *
213 * Unroll last two steps.
214 *
215  dnm2 = d
216  dmin2 = dmin
217  j4 = 4*( n0-2 ) - pp
218  j4p2 = j4 + 2*pp - 1
219  z( j4-2 ) = dnm2 + z( j4p2 )
220  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
221  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
222  dmin = min( dmin, dnm1 )
223 *
224  dmin1 = dmin
225  j4 = j4 + 4
226  j4p2 = j4 + 2*pp - 1
227  z( j4-2 ) = dnm1 + z( j4p2 )
228  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
229  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
230  dmin = min( dmin, dn )
231 *
232  ELSE
233 *
234 * Code for non IEEE arithmetic.
235 *
236  IF( pp.EQ.0 ) THEN
237  DO 30 j4 = 4*i0, 4*( n0-3 ), 4
238  z( j4-2 ) = d + z( j4-1 )
239  IF( d.LT.zero ) THEN
240  RETURN
241  ELSE
242  z( j4 ) = z( j4+1 )*( z( j4-1 ) / z( j4-2 ) )
243  d = z( j4+1 )*( d / z( j4-2 ) ) - tau
244  END IF
245  dmin = min( dmin, d )
246  emin = min( emin, z( j4 ) )
247  30 CONTINUE
248  ELSE
249  DO 40 j4 = 4*i0, 4*( n0-3 ), 4
250  z( j4-3 ) = d + z( j4 )
251  IF( d.LT.zero ) THEN
252  RETURN
253  ELSE
254  z( j4-1 ) = z( j4+2 )*( z( j4 ) / z( j4-3 ) )
255  d = z( j4+2 )*( d / z( j4-3 ) ) - tau
256  END IF
257  dmin = min( dmin, d )
258  emin = min( emin, z( j4-1 ) )
259  40 CONTINUE
260  END IF
261 *
262 * Unroll last two steps.
263 *
264  dnm2 = d
265  dmin2 = dmin
266  j4 = 4*( n0-2 ) - pp
267  j4p2 = j4 + 2*pp - 1
268  z( j4-2 ) = dnm2 + z( j4p2 )
269  IF( dnm2.LT.zero ) THEN
270  RETURN
271  ELSE
272  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
273  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
274  END IF
275  dmin = min( dmin, dnm1 )
276 *
277  dmin1 = dmin
278  j4 = j4 + 4
279  j4p2 = j4 + 2*pp - 1
280  z( j4-2 ) = dnm1 + z( j4p2 )
281  IF( dnm1.LT.zero ) THEN
282  RETURN
283  ELSE
284  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
285  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
286  END IF
287  dmin = min( dmin, dn )
288 *
289  END IF
290  ELSE
291 * This is the version that sets d's to zero if they are small enough
292  j4 = 4*i0 + pp - 3
293  emin = z( j4+4 )
294  d = z( j4 ) - tau
295  dmin = d
296  dmin1 = -z( j4 )
297  IF( ieee ) THEN
298 *
299 * Code for IEEE arithmetic.
300 *
301  IF( pp.EQ.0 ) THEN
302  DO 50 j4 = 4*i0, 4*( n0-3 ), 4
303  z( j4-2 ) = d + z( j4-1 )
304  temp = z( j4+1 ) / z( j4-2 )
305  d = d*temp - tau
306  IF( d.LT.dthresh ) d = zero
307  dmin = min( dmin, d )
308  z( j4 ) = z( j4-1 )*temp
309  emin = min( z( j4 ), emin )
310  50 CONTINUE
311  ELSE
312  DO 60 j4 = 4*i0, 4*( n0-3 ), 4
313  z( j4-3 ) = d + z( j4 )
314  temp = z( j4+2 ) / z( j4-3 )
315  d = d*temp - tau
316  IF( d.LT.dthresh ) d = zero
317  dmin = min( dmin, d )
318  z( j4-1 ) = z( j4 )*temp
319  emin = min( z( j4-1 ), emin )
320  60 CONTINUE
321  END IF
322 *
323 * Unroll last two steps.
324 *
325  dnm2 = d
326  dmin2 = dmin
327  j4 = 4*( n0-2 ) - pp
328  j4p2 = j4 + 2*pp - 1
329  z( j4-2 ) = dnm2 + z( j4p2 )
330  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
331  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
332  dmin = min( dmin, dnm1 )
333 *
334  dmin1 = dmin
335  j4 = j4 + 4
336  j4p2 = j4 + 2*pp - 1
337  z( j4-2 ) = dnm1 + z( j4p2 )
338  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
339  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
340  dmin = min( dmin, dn )
341 *
342  ELSE
343 *
344 * Code for non IEEE arithmetic.
345 *
346  IF( pp.EQ.0 ) THEN
347  DO 70 j4 = 4*i0, 4*( n0-3 ), 4
348  z( j4-2 ) = d + z( j4-1 )
349  IF( d.LT.zero ) THEN
350  RETURN
351  ELSE
352  z( j4 ) = z( j4+1 )*( z( j4-1 ) / z( j4-2 ) )
353  d = z( j4+1 )*( d / z( j4-2 ) ) - tau
354  END IF
355  IF( d.LT.dthresh) d = zero
356  dmin = min( dmin, d )
357  emin = min( emin, z( j4 ) )
358  70 CONTINUE
359  ELSE
360  DO 80 j4 = 4*i0, 4*( n0-3 ), 4
361  z( j4-3 ) = d + z( j4 )
362  IF( d.LT.zero ) THEN
363  RETURN
364  ELSE
365  z( j4-1 ) = z( j4+2 )*( z( j4 ) / z( j4-3 ) )
366  d = z( j4+2 )*( d / z( j4-3 ) ) - tau
367  END IF
368  IF( d.LT.dthresh) d = zero
369  dmin = min( dmin, d )
370  emin = min( emin, z( j4-1 ) )
371  80 CONTINUE
372  END IF
373 *
374 * Unroll last two steps.
375 *
376  dnm2 = d
377  dmin2 = dmin
378  j4 = 4*( n0-2 ) - pp
379  j4p2 = j4 + 2*pp - 1
380  z( j4-2 ) = dnm2 + z( j4p2 )
381  IF( dnm2.LT.zero ) THEN
382  RETURN
383  ELSE
384  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
385  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
386  END IF
387  dmin = min( dmin, dnm1 )
388 *
389  dmin1 = dmin
390  j4 = j4 + 4
391  j4p2 = j4 + 2*pp - 1
392  z( j4-2 ) = dnm1 + z( j4p2 )
393  IF( dnm1.LT.zero ) THEN
394  RETURN
395  ELSE
396  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
397  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
398  END IF
399  dmin = min( dmin, dn )
400 *
401  END IF
402  END IF
403 *
404  z( j4+2 ) = dn
405  z( 4*n0-pp ) = emin
406  RETURN
407 *
408 * End of DLASQ5
409 *

Here is the caller graph for this function: