133 SUBROUTINE zsteqr( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
145 DOUBLE PRECISION d( * ), e( * ), work( * )
146 COMPLEX*16 z( ldz, * )
152 DOUBLE PRECISION zero, one, two, three
153 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
155 COMPLEX*16 czero, cone
156 parameter( czero = ( 0.0d0, 0.0d0 ),
157 $ cone = ( 1.0d0, 0.0d0 ) )
159 parameter( maxit = 30 )
162 INTEGER i, icompz, ii, iscale,
j, jtot, k, l, l1, lend,
163 $ lendm1, lendp1, lendsv, lm1, lsv, m, mm, mm1,
165 DOUBLE PRECISION anorm,
b, c, eps, eps2, f, g, p, r, rt1, rt2,
166 $ s, safmax, safmin, ssfmax, ssfmin, tst
178 INTRINSIC abs, max, sign, sqrt
186 IF(
lsame( compz,
'N' ) )
THEN
188 ELSE IF(
lsame( compz,
'V' ) )
THEN
190 ELSE IF(
lsame( compz,
'I' ) )
THEN
195 IF( icompz.LT.0 )
THEN
197 ELSE IF( n.LT.0 )
THEN
199 ELSE IF( ( ldz.LT.1 ) .OR. ( icompz.GT.0 .AND. ldz.LT.max( 1,
204 CALL
xerbla(
'ZSTEQR', -info )
224 safmax = one / safmin
225 ssfmax = sqrt( safmax ) / three
226 ssfmin = sqrt( safmin ) / eps2
232 $ CALL
zlaset(
'Full', n, n, czero, cone, z, ldz )
254 IF( tst.LE.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
255 $ 1 ) ) ) )*eps )
THEN
274 anorm =
dlanst(
'I', lend-l+1, d( l ), e( l ) )
278 IF( anorm.GT.ssfmax )
THEN
280 CALL
dlascl(
'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
282 CALL
dlascl(
'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
284 ELSE IF( anorm.LT.ssfmin )
THEN
286 CALL
dlascl(
'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
288 CALL
dlascl(
'G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
294 IF( abs( d( lend ) ).LT.abs( d( l ) ) )
THEN
309 tst = abs( e( m ) )**2
310 IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m+1 ) )+
328 IF( icompz.GT.0 )
THEN
329 CALL
dlaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
332 CALL
zlasr(
'R',
'V',
'B', n, 2, work( l ),
333 $ work( n-1+l ), z( 1, l ), ldz )
335 CALL
dlae2( d( l ), e( l ), d( l+1 ), rt1, rt2 )
352 g = ( d( l+1 )-p ) / ( two*e( l ) )
354 g = d( m ) - p + ( e( l ) / ( g+sign( r, g ) ) )
366 CALL
dlartg( g, f, c, s, r )
370 r = ( d( i )-g )*s + two*c*
b
377 IF( icompz.GT.0 )
THEN
386 IF( icompz.GT.0 )
THEN
388 CALL
zlasr(
'R',
'V',
'B', n, mm, work( l ), work( n-1+l ),
415 DO 100 m = l, lendp1, -1
416 tst = abs( e( m-1 ) )**2
417 IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m-1 ) )+
435 IF( icompz.GT.0 )
THEN
436 CALL
dlaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )
439 CALL
zlasr(
'R',
'V',
'F', n, 2, work( m ),
440 $ work( n-1+m ), z( 1, l-1 ), ldz )
442 CALL
dlae2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2 )
459 g = ( d( l-1 )-p ) / ( two*e( l-1 ) )
461 g = d( m ) - p + ( e( l-1 ) / ( g+sign( r, g ) ) )
473 CALL
dlartg( g, f, c, s, r )
477 r = ( d( i+1 )-g )*s + two*c*
b
484 IF( icompz.GT.0 )
THEN
493 IF( icompz.GT.0 )
THEN
495 CALL
zlasr(
'R',
'V',
'F', n, mm, work( m ), work( n-1+m ),
518 IF( iscale.EQ.1 )
THEN
519 CALL
dlascl(
'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
520 $ d( lsv ), n, info )
521 CALL
dlascl(
'G', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
523 ELSE IF( iscale.EQ.2 )
THEN
524 CALL
dlascl(
'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
525 $ d( lsv ), n, info )
526 CALL
dlascl(
'G', 0, 0, ssfmin, anorm, lendsv-lsv, 1, e( lsv ),
533 IF( jtot.EQ.nmaxit )
THEN
545 IF( icompz.EQ.0 )
THEN
549 CALL
dlasrt(
'I', n, d, info )
560 IF( d(
j ).LT.p )
THEN
568 CALL
zswap( n, z( 1, i ), 1, z( 1, k ), 1 )
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine zlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
ZLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine dlae2(A, B, C, RT1, RT2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
logical function lsame(CA, CB)
LSAME
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dlaev2(A, B, C, RT1, RT2, CS1, SN1)
DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
double precision function dlanst(NORM, N, D, E)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j