LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ 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: