LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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: