131 parameter ( cbias = 1.50e0 )
132 REAL zero, half, one, two, four, hundrd
133 parameter ( zero = 0.0e0, half = 0.5e0, one = 1.0e0,
134 $ two = 2.0e0, four = 4.0e0, hundrd = 100.0e0 )
138 INTEGER i0, i4, iinfo, ipn4, iter, iwhila, iwhilb, k,
139 $ kmin, n0, nbig, ndiv, nfail, pp, splt, ttype,
141 REAL d, dee, deemin, desig, dmin, dmin1, dmin2, dn,
142 $ dn1, dn2, e, emax, emin, eps, g, oldemn, qmax,
143 $ qmin, s, safmin, sigma, t, tau, temp, tol,
144 $ tol2, trace, zmax, tempe, tempq
154 INTRINSIC abs, max, min,
REAL, sqrt
162 eps =
slamch(
'Precision' )
163 safmin =
slamch(
'Safe minimum' )
169 CALL xerbla(
'SLASQ2', 1 )
171 ELSE IF( n.EQ.0 )
THEN
173 ELSE IF( n.EQ.1 )
THEN
177 IF( z( 1 ).LT.zero )
THEN
179 CALL xerbla(
'SLASQ2', 2 )
182 ELSE IF( n.EQ.2 )
THEN
186 IF( z( 2 ).LT.zero .OR. z( 3 ).LT.zero )
THEN
188 CALL xerbla(
'SLASQ2', 2 )
190 ELSE IF( z( 3 ).GT.z( 1 ) )
THEN
195 z( 5 ) = z( 1 ) + z( 2 ) + z( 3 )
196 IF( z( 2 ).GT.z( 3 )*tol2 )
THEN
197 t = half*( ( z( 1 )-z( 3 ) )+z( 2 ) )
198 s = z( 3 )*( z( 2 ) / t )
200 s = z( 3 )*( z( 2 ) / ( t*( one+sqrt( one+s / t ) ) ) )
202 s = z( 3 )*( z( 2 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
204 t = z( 1 ) + ( s+z( 2 ) )
205 z( 3 ) = z( 3 )*( z( 1 ) / t )
209 z( 6 ) = z( 2 ) + z( 1 )
222 DO 10 k = 1, 2*( n-1 ), 2
223 IF( z( k ).LT.zero )
THEN
225 CALL xerbla(
'SLASQ2', 2 )
227 ELSE IF( z( k+1 ).LT.zero )
THEN
229 CALL xerbla(
'SLASQ2', 2 )
234 qmax = max( qmax, z( k ) )
235 emin = min( emin, z( k+1 ) )
236 zmax = max( qmax, zmax, z( k+1 ) )
238 IF( z( 2*n-1 ).LT.zero )
THEN
239 info = -( 200+2*n-1 )
240 CALL xerbla(
'SLASQ2', 2 )
244 qmax = max( qmax, z( 2*n-1 ) )
245 zmax = max( qmax, zmax )
253 CALL slasrt(
'D', n, z, iinfo )
262 IF( trace.EQ.zero )
THEN
283 z( 2*k-3 ) = z( k-1 )
291 IF( cbias*z( 4*i0-3 ).LT.z( 4*n0-3 ) )
THEN
293 DO 40 i4 = 4*i0, 2*( i0+n0-1 ), 4
295 z( i4-3 ) = z( ipn4-i4-3 )
296 z( ipn4-i4-3 ) = temp
298 z( i4-1 ) = z( ipn4-i4-5 )
299 z( ipn4-i4-5 ) = temp
310 DO 50 i4 = 4*( n0-1 ) + pp, 4*i0 + pp, -4
311 IF( z( i4-1 ).LE.tol2*d )
THEN
315 d = z( i4-3 )*( d / ( d+z( i4-1 ) ) )
321 emin = z( 4*i0+pp+1 )
323 DO 60 i4 = 4*i0 + pp, 4*( n0-1 ) + pp, 4
324 z( i4-2*pp-2 ) = d + z( i4-1 )
325 IF( z( i4-1 ).LE.tol2*d )
THEN
330 ELSE IF( safmin*z( i4+1 ).LT.z( i4-2*pp-2 ) .AND.
331 $ safmin*z( i4-2*pp-2 ).LT.z( i4+1 ) )
THEN
332 temp = z( i4+1 ) / z( i4-2*pp-2 )
333 z( i4-2*pp ) = z( i4-1 )*temp
336 z( i4-2*pp ) = z( i4+1 )*( z( i4-1 ) / z( i4-2*pp-2 ) )
337 d = z( i4+1 )*( d / z( i4-2*pp-2 ) )
339 emin = min( emin, z( i4-2*pp ) )
345 qmax = z( 4*i0-pp-2 )
346 DO 70 i4 = 4*i0 - pp + 2, 4*n0 - pp - 2, 4
347 qmax = max( qmax, z( i4 ) )
370 DO 160 iwhila = 1, n + 1
385 IF( sigma.LT.zero )
THEN
395 emin = abs( z( 4*n0-5 ) )
401 DO 90 i4 = 4*n0, 8, -4
402 IF( z( i4-5 ).LE.zero )
404 IF( qmin.GE.four*emax )
THEN
405 qmin = min( qmin, z( i4-3 ) )
406 emax = max( emax, z( i4-5 ) )
408 qmax = max( qmax, z( i4-7 )+z( i4-5 ) )
409 emin = min( emin, z( i4-5 ) )
417 IF( n0-i0.GT.1 )
THEN
421 DO 110 i4 = 4*i0+1, 4*n0-3, 4
422 dee = z( i4 )*( dee /( dee+z( i4-2 ) ) )
423 IF( dee.LE.deemin )
THEN
428 IF( (kmin-i0)*2.LT.n0-kmin .AND.
429 $ deemin.LE.half*z(4*n0-3) )
THEN
432 DO 120 i4 = 4*i0, 2*( i0+n0-1 ), 4
434 z( i4-3 ) = z( ipn4-i4-3 )
435 z( ipn4-i4-3 ) = temp
437 z( i4-2 ) = z( ipn4-i4-2 )
438 z( ipn4-i4-2 ) = temp
440 z( i4-1 ) = z( ipn4-i4-5 )
441 z( ipn4-i4-5 ) = temp
443 z( i4 ) = z( ipn4-i4-4 )
444 z( ipn4-i4-4 ) = temp
451 dmin = -max( zero, qmin-two*sqrt( qmin )*sqrt( emax ) )
459 nbig = 100*( n0-i0+1 )
460 DO 140 iwhilb = 1, nbig
466 CALL slasq3( i0, n0, z, pp, dmin, sigma, desig, qmax, nfail,
467 $ iter, ndiv, ieee, ttype, dmin1, dmin2, dn, dn1,
474 IF( pp.EQ.0 .AND. n0-i0.GE.3 )
THEN
475 IF( z( 4*n0 ).LE.tol2*qmax .OR.
476 $ z( 4*n0-1 ).LE.tol2*sigma )
THEN
481 DO 130 i4 = 4*i0, 4*( n0-3 ), 4
482 IF( z( i4 ).LE.tol2*z( i4-3 ) .OR.
483 $ z( i4-1 ).LE.tol2*sigma )
THEN
490 qmax = max( qmax, z( i4+1 ) )
491 emin = min( emin, z( i4-1 ) )
492 oldemn = min( oldemn, z( i4 ) )
513 z( 4*i0-3 ) = z( 4*i0-3 ) + sigma
516 z( 4*k-5 ) = z( 4*k-5 ) * (tempq / z( 4*k-7 ))
518 z( 4*k-3 ) = z( 4*k-3 ) + sigma + tempe - z( 4*k-5 )
525 DO WHILE( ( i1.GE.2 ) .AND. ( z(4*i1-5).GE.zero ) )
535 z( 2*k-1 ) = z( 4*k-3 )
542 z( 2*k ) = z( 4*k-1 )
570 CALL slasrt(
'D', n, z, iinfo )
581 z( 2*n+3 ) =
REAL( iter )
582 z( 2*n+4 ) =
REAL( NDIV ) /
REAL( n**2 )
583 z( 2*n+5 ) = hundrd*nfail /
REAL( iter )
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
subroutine xerbla(SRNAME, INFO)
XERBLA
real function slamch(CMACH)
SLAMCH
subroutine slasq3(I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, DN2, G, TAU)
SLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.