LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine slasq5 ( integer  I0,
integer  N0,
real, dimension( * )  Z,
integer  PP,
real  TAU,
real  SIGMA,
real  DMIN,
real  DMIN1,
real  DMIN2,
real  DN,
real  DNM1,
real  DNM2,
logical  IEEE,
real  EPS 
)

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

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

Purpose: 
 SLASQ5 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 REAL 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 REAL
        This is the shift.
[in]SIGMASIGMA is REAL This is the accumulated shift up to this step.
[out]DMIN
          DMIN is REAL
        Minimum value of d.
[out]DMIN1
          DMIN1 is REAL
        Minimum value of d, excluding D( N0 ).
[out]DMIN2
          DMIN2 is REAL
        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
[out]DN
          DN is REAL
        d(N0), the last value of d.
[out]DNM1
          DNM1 is REAL
        d(N0-1).
[out]DNM2
          DNM2 is REAL
        d(N0-2).
[in]IEEE
          IEEE is LOGICAL
        Flag for IEEE or non IEEE arithmetic.
[in]EPS
         EPS is REAL
        This is the value of epsilon used.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015

Definition at line 145 of file slasq5.f.

145 *
146 * -- LAPACK computational routine (version 3.6.0) --
147 * -- LAPACK is a software package provided by Univ. of Tennessee, --
148 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
149 * November 2015
150 *
151 * .. Scalar Arguments ..
152  LOGICAL ieee
153  INTEGER i0, n0, pp
154  REAL dmin, dmin1, dmin2, dn, dnm1, dnm2, tau,
155  $ sigma, eps
156 * ..
157 * .. Array Arguments ..
158  REAL z( * )
159 * ..
160 *
161 * =====================================================================
162 *
163 * .. Parameter ..
164  REAL zero, half
165  parameter ( zero = 0.0e0, half = 0.5 )
166 * ..
167 * .. Local Scalars ..
168  INTEGER j4, j4p2
169  REAL d, emin, temp, dthresh
170 * ..
171 * .. Intrinsic Functions ..
172  INTRINSIC min
173 * ..
174 * .. Executable Statements ..
175 *
176  IF( ( n0-i0-1 ).LE.0 )
177  $ RETURN
178 *
179  dthresh = eps*(sigma+tau)
180  IF( tau.LT.dthresh*half ) tau = zero
181  IF( tau.NE.zero ) THEN
182  j4 = 4*i0 + pp - 3
183  emin = z( j4+4 )
184  d = z( j4 ) - tau
185  dmin = d
186  dmin1 = -z( j4 )
187 *
188  IF( ieee ) THEN
189 *
190 * Code for IEEE arithmetic.
191 *
192  IF( pp.EQ.0 ) THEN
193  DO 10 j4 = 4*i0, 4*( n0-3 ), 4
194  z( j4-2 ) = d + z( j4-1 )
195  temp = z( j4+1 ) / z( j4-2 )
196  d = d*temp - tau
197  dmin = min( dmin, d )
198  z( j4 ) = z( j4-1 )*temp
199  emin = min( z( j4 ), emin )
200  10 CONTINUE
201  ELSE
202  DO 20 j4 = 4*i0, 4*( n0-3 ), 4
203  z( j4-3 ) = d + z( j4 )
204  temp = z( j4+2 ) / z( j4-3 )
205  d = d*temp - tau
206  dmin = min( dmin, d )
207  z( j4-1 ) = z( j4 )*temp
208  emin = min( z( j4-1 ), emin )
209  20 CONTINUE
210  END IF
211 *
212 * Unroll last two steps.
213 *
214  dnm2 = d
215  dmin2 = dmin
216  j4 = 4*( n0-2 ) - pp
217  j4p2 = j4 + 2*pp - 1
218  z( j4-2 ) = dnm2 + z( j4p2 )
219  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
220  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
221  dmin = min( dmin, dnm1 )
222 *
223  dmin1 = dmin
224  j4 = j4 + 4
225  j4p2 = j4 + 2*pp - 1
226  z( j4-2 ) = dnm1 + z( j4p2 )
227  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
228  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
229  dmin = min( dmin, dn )
230 *
231  ELSE
232 *
233 * Code for non IEEE arithmetic.
234 *
235  IF( pp.EQ.0 ) THEN
236  DO 30 j4 = 4*i0, 4*( n0-3 ), 4
237  z( j4-2 ) = d + z( j4-1 )
238  IF( d.LT.zero ) THEN
239  RETURN
240  ELSE
241  z( j4 ) = z( j4+1 )*( z( j4-1 ) / z( j4-2 ) )
242  d = z( j4+1 )*( d / z( j4-2 ) ) - tau
243  END IF
244  dmin = min( dmin, d )
245  emin = min( emin, z( j4 ) )
246  30 CONTINUE
247  ELSE
248  DO 40 j4 = 4*i0, 4*( n0-3 ), 4
249  z( j4-3 ) = d + z( j4 )
250  IF( d.LT.zero ) THEN
251  RETURN
252  ELSE
253  z( j4-1 ) = z( j4+2 )*( z( j4 ) / z( j4-3 ) )
254  d = z( j4+2 )*( d / z( j4-3 ) ) - tau
255  END IF
256  dmin = min( dmin, d )
257  emin = min( emin, z( j4-1 ) )
258  40 CONTINUE
259  END IF
260 *
261 * Unroll last two steps.
262 *
263  dnm2 = d
264  dmin2 = dmin
265  j4 = 4*( n0-2 ) - pp
266  j4p2 = j4 + 2*pp - 1
267  z( j4-2 ) = dnm2 + z( j4p2 )
268  IF( dnm2.LT.zero ) THEN
269  RETURN
270  ELSE
271  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
272  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
273  END IF
274  dmin = min( dmin, dnm1 )
275 *
276  dmin1 = dmin
277  j4 = j4 + 4
278  j4p2 = j4 + 2*pp - 1
279  z( j4-2 ) = dnm1 + z( j4p2 )
280  IF( dnm1.LT.zero ) THEN
281  RETURN
282  ELSE
283  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
284  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
285  END IF
286  dmin = min( dmin, dn )
287 *
288  END IF
289 *
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 *
403  END IF
404  z( j4+2 ) = dn
405  z( 4*n0-pp ) = emin
406  RETURN
407 *
408 * End of SLASQ5
409 *

Here is the caller graph for this function: