LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slasq5()

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.

Definition at line 142 of file slasq5.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  REAL DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
153  $ SIGMA, EPS
154 * ..
155 * .. Array Arguments ..
156  REAL Z( * )
157 * ..
158 *
159 * =====================================================================
160 *
161 * .. Parameter ..
162  REAL ZERO, HALF
163  parameter( zero = 0.0e0, half = 0.5 )
164 * ..
165 * .. Local Scalars ..
166  INTEGER J4, J4P2
167  REAL 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 *
288  ELSE
289 * This is the version that sets d's to zero if they are small enough
290  j4 = 4*i0 + pp - 3
291  emin = z( j4+4 )
292  d = z( j4 ) - tau
293  dmin = d
294  dmin1 = -z( j4 )
295  IF( ieee ) THEN
296 *
297 * Code for IEEE arithmetic.
298 *
299  IF( pp.EQ.0 ) THEN
300  DO 50 j4 = 4*i0, 4*( n0-3 ), 4
301  z( j4-2 ) = d + z( j4-1 )
302  temp = z( j4+1 ) / z( j4-2 )
303  d = d*temp - tau
304  IF( d.LT.dthresh ) d = zero
305  dmin = min( dmin, d )
306  z( j4 ) = z( j4-1 )*temp
307  emin = min( z( j4 ), emin )
308  50 CONTINUE
309  ELSE
310  DO 60 j4 = 4*i0, 4*( n0-3 ), 4
311  z( j4-3 ) = d + z( j4 )
312  temp = z( j4+2 ) / z( j4-3 )
313  d = d*temp - tau
314  IF( d.LT.dthresh ) d = zero
315  dmin = min( dmin, d )
316  z( j4-1 ) = z( j4 )*temp
317  emin = min( z( j4-1 ), emin )
318  60 CONTINUE
319  END IF
320 *
321 * Unroll last two steps.
322 *
323  dnm2 = d
324  dmin2 = dmin
325  j4 = 4*( n0-2 ) - pp
326  j4p2 = j4 + 2*pp - 1
327  z( j4-2 ) = dnm2 + z( j4p2 )
328  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
329  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
330  dmin = min( dmin, dnm1 )
331 *
332  dmin1 = dmin
333  j4 = j4 + 4
334  j4p2 = j4 + 2*pp - 1
335  z( j4-2 ) = dnm1 + z( j4p2 )
336  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
337  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
338  dmin = min( dmin, dn )
339 *
340  ELSE
341 *
342 * Code for non IEEE arithmetic.
343 *
344  IF( pp.EQ.0 ) THEN
345  DO 70 j4 = 4*i0, 4*( n0-3 ), 4
346  z( j4-2 ) = d + z( j4-1 )
347  IF( d.LT.zero ) THEN
348  RETURN
349  ELSE
350  z( j4 ) = z( j4+1 )*( z( j4-1 ) / z( j4-2 ) )
351  d = z( j4+1 )*( d / z( j4-2 ) ) - tau
352  END IF
353  IF( d.LT.dthresh ) d = zero
354  dmin = min( dmin, d )
355  emin = min( emin, z( j4 ) )
356  70 CONTINUE
357  ELSE
358  DO 80 j4 = 4*i0, 4*( n0-3 ), 4
359  z( j4-3 ) = d + z( j4 )
360  IF( d.LT.zero ) THEN
361  RETURN
362  ELSE
363  z( j4-1 ) = z( j4+2 )*( z( j4 ) / z( j4-3 ) )
364  d = z( j4+2 )*( d / z( j4-3 ) ) - tau
365  END IF
366  IF( d.LT.dthresh ) d = zero
367  dmin = min( dmin, d )
368  emin = min( emin, z( j4-1 ) )
369  80 CONTINUE
370  END IF
371 *
372 * Unroll last two steps.
373 *
374  dnm2 = d
375  dmin2 = dmin
376  j4 = 4*( n0-2 ) - pp
377  j4p2 = j4 + 2*pp - 1
378  z( j4-2 ) = dnm2 + z( j4p2 )
379  IF( dnm2.LT.zero ) THEN
380  RETURN
381  ELSE
382  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
383  dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) ) - tau
384  END IF
385  dmin = min( dmin, dnm1 )
386 *
387  dmin1 = dmin
388  j4 = j4 + 4
389  j4p2 = j4 + 2*pp - 1
390  z( j4-2 ) = dnm1 + z( j4p2 )
391  IF( dnm1.LT.zero ) THEN
392  RETURN
393  ELSE
394  z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
395  dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) ) - tau
396  END IF
397  dmin = min( dmin, dn )
398 *
399  END IF
400 *
401  END IF
402  z( j4+2 ) = dn
403  z( 4*n0-pp ) = emin
404  RETURN
405 *
406 * End of SLASQ5
407 *
Here is the caller graph for this function: