LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
slasq5.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine slasq5 (I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, IEEE, EPS)
 SLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr. More...
 

Function/Subroutine Documentation

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.
[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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 134 of file slasq5.f.

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

Here is the caller graph for this function: