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] 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.```
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: