LAPACK 3.12.1
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 140 of file slasq5.f.

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