LAPACK  3.7.0
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]SIGMA
          SIGMA 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
December 2016

Definition at line 146 of file slasq5.f.

146 *
147 * -- LAPACK computational routine (version 3.7.0) --
148 * -- LAPACK is a software package provided by Univ. of Tennessee, --
149 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
150 * December 2016
151 *
152 * .. Scalar Arguments ..
153  LOGICAL ieee
154  INTEGER i0, n0, pp
155  REAL dmin, dmin1, dmin2, dn, dnm1, dnm2, tau,
156  $ sigma, eps
157 * ..
158 * .. Array Arguments ..
159  REAL z( * )
160 * ..
161 *
162 * =====================================================================
163 *
164 * .. Parameter ..
165  REAL zero, half
166  parameter ( zero = 0.0e0, half = 0.5 )
167 * ..
168 * .. Local Scalars ..
169  INTEGER j4, j4p2
170  REAL 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 *
291  ELSE
292 * This is the version that sets d's to zero if they are small enough
293  j4 = 4*i0 + pp - 3
294  emin = z( j4+4 )
295  d = z( j4 ) - tau
296  dmin = d
297  dmin1 = -z( j4 )
298  IF( ieee ) THEN
299 *
300 * Code for IEEE arithmetic.
301 *
302  IF( pp.EQ.0 ) THEN
303  DO 50 j4 = 4*i0, 4*( n0-3 ), 4
304  z( j4-2 ) = d + z( j4-1 )
305  temp = z( j4+1 ) / z( j4-2 )
306  d = d*temp - tau
307  IF( d.LT.dthresh ) d = zero
308  dmin = min( dmin, d )
309  z( j4 ) = z( j4-1 )*temp
310  emin = min( z( j4 ), emin )
311  50 CONTINUE
312  ELSE
313  DO 60 j4 = 4*i0, 4*( n0-3 ), 4
314  z( j4-3 ) = d + z( j4 )
315  temp = z( j4+2 ) / z( j4-3 )
316  d = d*temp - tau
317  IF( d.LT.dthresh ) d = zero
318  dmin = min( dmin, d )
319  z( j4-1 ) = z( j4 )*temp
320  emin = min( z( j4-1 ), emin )
321  60 CONTINUE
322  END IF
323 *
324 * Unroll last two steps.
325 *
326  dnm2 = d
327  dmin2 = dmin
328  j4 = 4*( n0-2 ) - pp
329  j4p2 = j4 + 2*pp - 1
330  z( j4-2 ) = dnm2 + z( j4p2 )
331  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
332  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
333  dmin = min( dmin, dnm1 )
334 *
335  dmin1 = dmin
336  j4 = j4 + 4
337  j4p2 = j4 + 2*pp - 1
338  z( j4-2 ) = dnm1 + z( j4p2 )
339  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
340  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
341  dmin = min( dmin, dn )
342 *
343  ELSE
344 *
345 * Code for non IEEE arithmetic.
346 *
347  IF( pp.EQ.0 ) THEN
348  DO 70 j4 = 4*i0, 4*( n0-3 ), 4
349  z( j4-2 ) = d + z( j4-1 )
350  IF( d.LT.zero ) THEN
351  RETURN
352  ELSE
353  z( j4 ) = z( j4+1 )*( z( j4-1 ) / z( j4-2 ) )
354  d = z( j4+1 )*( d / z( j4-2 ) ) - tau
355  END IF
356  IF( d.LT.dthresh ) d = zero
357  dmin = min( dmin, d )
358  emin = min( emin, z( j4 ) )
359  70 CONTINUE
360  ELSE
361  DO 80 j4 = 4*i0, 4*( n0-3 ), 4
362  z( j4-3 ) = d + z( j4 )
363  IF( d.LT.zero ) THEN
364  RETURN
365  ELSE
366  z( j4-1 ) = z( j4+2 )*( z( j4 ) / z( j4-3 ) )
367  d = z( j4+2 )*( d / z( j4-3 ) ) - tau
368  END IF
369  IF( d.LT.dthresh ) d = zero
370  dmin = min( dmin, d )
371  emin = min( emin, z( j4-1 ) )
372  80 CONTINUE
373  END IF
374 *
375 * Unroll last two steps.
376 *
377  dnm2 = d
378  dmin2 = dmin
379  j4 = 4*( n0-2 ) - pp
380  j4p2 = j4 + 2*pp - 1
381  z( j4-2 ) = dnm2 + z( j4p2 )
382  IF( dnm2.LT.zero ) THEN
383  RETURN
384  ELSE
385  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
386  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
387  END IF
388  dmin = min( dmin, dnm1 )
389 *
390  dmin1 = dmin
391  j4 = j4 + 4
392  j4p2 = j4 + 2*pp - 1
393  z( j4-2 ) = dnm1 + z( j4p2 )
394  IF( dnm1.LT.zero ) THEN
395  RETURN
396  ELSE
397  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
398  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
399  END IF
400  dmin = min( dmin, dn )
401 *
402  END IF
403 *
404  END IF
405  z( j4+2 ) = dn
406  z( 4*n0-pp ) = emin
407  RETURN
408 *
409 * End of SLASQ5
410 *

Here is the caller graph for this function: