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
155 INTRINSIC abs, max, min,
REAL, sqrt
163 eps =
slamch(
'Precision' )
164 safmin =
slamch(
'Safe minimum' )
170 CALL
xerbla(
'SLASQ2', 1 )
172 ELSE IF( n.EQ.0 )
THEN
174 ELSE IF( n.EQ.1 )
THEN
178 IF( z( 1 ).LT.zero )
THEN
180 CALL
xerbla(
'SLASQ2', 2 )
183 ELSE IF( n.EQ.2 )
THEN
187 IF( z( 2 ).LT.zero .OR. z( 3 ).LT.zero )
THEN
189 CALL
xerbla(
'SLASQ2', 2 )
191 ELSE IF( z( 3 ).GT.z( 1 ) )
THEN
196 z( 5 ) = z( 1 ) + z( 2 ) + z( 3 )
197 IF( z( 2 ).GT.z( 3 )*tol2 )
THEN
198 t = half*( ( z( 1 )-z( 3 ) )+z( 2 ) )
199 s = z( 3 )*( z( 2 ) / t )
201 s = z( 3 )*( z( 2 ) / ( t*( one+sqrt( one+s / t ) ) ) )
203 s = z( 3 )*( z( 2 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
205 t = z( 1 ) + ( s+z( 2 ) )
206 z( 3 ) = z( 3 )*( z( 1 ) / t )
210 z( 6 ) = z( 2 ) + z( 1 )
223 DO 10 k = 1, 2*( n-1 ), 2
224 IF( z( k ).LT.zero )
THEN
226 CALL
xerbla(
'SLASQ2', 2 )
228 ELSE IF( z( k+1 ).LT.zero )
THEN
230 CALL
xerbla(
'SLASQ2', 2 )
235 qmax = max( qmax, z( k ) )
236 emin = min( emin, z( k+1 ) )
237 zmax = max( qmax, zmax, z( k+1 ) )
239 IF( z( 2*n-1 ).LT.zero )
THEN
240 info = -( 200+2*n-1 )
241 CALL
xerbla(
'SLASQ2', 2 )
245 qmax = max( qmax, z( 2*n-1 ) )
246 zmax = max( qmax, zmax )
254 CALL
slasrt(
'D', n, z, iinfo )
263 IF( trace.EQ.zero )
THEN
284 z( 2*k-3 ) = z( k-1 )
292 IF( cbias*z( 4*i0-3 ).LT.z( 4*n0-3 ) )
THEN
294 DO 40 i4 = 4*i0, 2*( i0+n0-1 ), 4
296 z( i4-3 ) = z( ipn4-i4-3 )
297 z( ipn4-i4-3 ) = temp
299 z( i4-1 ) = z( ipn4-i4-5 )
300 z( ipn4-i4-5 ) = temp
311 DO 50 i4 = 4*( n0-1 ) + pp, 4*i0 + pp, -4
312 IF( z( i4-1 ).LE.tol2*d )
THEN
316 d = z( i4-3 )*( d / ( d+z( i4-1 ) ) )
322 emin = z( 4*i0+pp+1 )
324 DO 60 i4 = 4*i0 + pp, 4*( n0-1 ) + pp, 4
325 z( i4-2*pp-2 ) = d + z( i4-1 )
326 IF( z( i4-1 ).LE.tol2*d )
THEN
331 ELSE IF( safmin*z( i4+1 ).LT.z( i4-2*pp-2 ) .AND.
332 $ safmin*z( i4-2*pp-2 ).LT.z( i4+1 ) )
THEN
333 temp = z( i4+1 ) / z( i4-2*pp-2 )
334 z( i4-2*pp ) = z( i4-1 )*temp
337 z( i4-2*pp ) = z( i4+1 )*( z( i4-1 ) / z( i4-2*pp-2 ) )
338 d = z( i4+1 )*( d / z( i4-2*pp-2 ) )
340 emin = min( emin, z( i4-2*pp ) )
346 qmax = z( 4*i0-pp-2 )
347 DO 70 i4 = 4*i0 - pp + 2, 4*n0 - pp - 2, 4
348 qmax = max( qmax, z( i4 ) )
371 DO 160 iwhila = 1, n + 1
386 IF( sigma.LT.zero )
THEN
396 emin = abs( z( 4*n0-5 ) )
402 DO 90 i4 = 4*n0, 8, -4
403 IF( z( i4-5 ).LE.zero )
405 IF( qmin.GE.four*emax )
THEN
406 qmin = min( qmin, z( i4-3 ) )
407 emax = max( emax, z( i4-5 ) )
409 qmax = max( qmax, z( i4-7 )+z( i4-5 ) )
410 emin = min( emin, z( i4-5 ) )
418 IF( n0-i0.GT.1 )
THEN
422 DO 110 i4 = 4*i0+1, 4*n0-3, 4
423 dee = z( i4 )*( dee /( dee+z( i4-2 ) ) )
424 IF( dee.LE.deemin )
THEN
429 IF( (kmin-i0)*2.LT.n0-kmin .AND.
430 $ deemin.LE.half*z(4*n0-3) )
THEN
433 DO 120 i4 = 4*i0, 2*( i0+n0-1 ), 4
435 z( i4-3 ) = z( ipn4-i4-3 )
436 z( ipn4-i4-3 ) = temp
438 z( i4-2 ) = z( ipn4-i4-2 )
439 z( ipn4-i4-2 ) = temp
441 z( i4-1 ) = z( ipn4-i4-5 )
442 z( ipn4-i4-5 ) = temp
444 z( i4 ) = z( ipn4-i4-4 )
445 z( ipn4-i4-4 ) = temp
452 dmin = -max( zero, qmin-two*sqrt( qmin )*sqrt( emax ) )
460 nbig = 100*( n0-i0+1 )
461 DO 140 iwhilb = 1, nbig
467 CALL
slasq3( i0, n0, z, pp, dmin, sigma, desig, qmax, nfail,
468 $ iter, ndiv, ieee, ttype, dmin1, dmin2, dn, dn1,
475 IF( pp.EQ.0 .AND. n0-i0.GE.3 )
THEN
476 IF( z( 4*n0 ).LE.tol2*qmax .OR.
477 $ z( 4*n0-1 ).LE.tol2*sigma )
THEN
482 DO 130 i4 = 4*i0, 4*( n0-3 ), 4
483 IF( z( i4 ).LE.tol2*z( i4-3 ) .OR.
484 $ z( i4-1 ).LE.tol2*sigma )
THEN
491 qmax = max( qmax, z( i4+1 ) )
492 emin = min( emin, z( i4-1 ) )
493 oldemn = min( oldemn, z( i4 ) )
514 z( 4*i0-3 ) = z( 4*i0-3 ) + sigma
517 z( 4*k-5 ) = z( 4*k-5 ) * (tempq / z( 4*k-7 ))
519 z( 4*k-3 ) = z( 4*k-3 ) + sigma + tempe - z( 4*k-5 )
526 DO WHILE( ( i1.GE.2 ) .AND. ( z(4*i1-5).GE.zero ) )
536 z( 2*k-1 ) = z( 4*k-3 )
543 z( 2*k ) = z( 4*k-1 )
571 CALL
slasrt(
'D', n, z, iinfo )
582 z( 2*n+3 ) =
REAL( iter )
583 z( 2*n+4 ) =
REAL( NDIV ) /
REAL( n**2 )
584 z( 2*n+5 ) = hundrd*nfail /
REAL( iter )
subroutine slasq2(N, Z, INFO)
SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
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.
REAL function slamch(CMACH)
SLAMCH
subroutine xerbla(SRNAME, INFO)
XERBLA
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.