189 SUBROUTINE dstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
199 INTEGER info, ldz, liwork, lwork, n
203 DOUBLE PRECISION d( * ), e( * ), work( * ), z( ldz, * )
209 DOUBLE PRECISION zero, one, two
210 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
214 INTEGER finish, i, icompz, ii,
j, k, lgn, liwmin,
215 $ lwmin, m, smlsiz, start, storez, strtrw
216 DOUBLE PRECISION eps, orgnrm, p, tiny
229 INTRINSIC abs, dble, int, log, max, mod, sqrt
236 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
238 IF(
lsame( compz,
'N' ) )
THEN
240 ELSE IF(
lsame( compz,
'V' ) )
THEN
242 ELSE IF(
lsame( compz,
'I' ) )
THEN
247 IF( icompz.LT.0 )
THEN
249 ELSE IF( n.LT.0 )
THEN
251 ELSE IF( ( ldz.LT.1 ) .OR.
252 $ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) )
THEN
260 smlsiz =
ilaenv( 9,
'DSTEDC',
' ', 0, 0, 0, 0 )
261 IF( n.LE.1 .OR. icompz.EQ.0 )
THEN
264 ELSE IF( n.LE.smlsiz )
THEN
268 lgn = int( log( dble( n ) )/log( two ) )
273 IF( icompz.EQ.1 )
THEN
274 lwmin = 1 + 3*n + 2*n*lgn + 4*n**2
275 liwmin = 6 + 6*n + 5*n*lgn
276 ELSE IF( icompz.EQ.2 )
THEN
277 lwmin = 1 + 4*n + n**2
284 IF( lwork.LT.lwmin .AND. .NOT. lquery )
THEN
286 ELSE IF( liwork.LT.liwmin .AND. .NOT. lquery )
THEN
292 CALL
xerbla(
'DSTEDC', -info )
294 ELSE IF (lquery)
THEN
319 IF( icompz.EQ.0 )
THEN
320 CALL
dsterf( n, d, e, info )
327 IF( n.LE.smlsiz )
THEN
329 CALL
dsteqr( compz, n, d, e, z, ldz, work, info )
336 IF( icompz.EQ.1 )
THEN
342 IF( icompz.EQ.2 )
THEN
343 CALL
dlaset(
'Full', n, n, zero, one, z, ldz )
348 orgnrm =
dlanst(
'M', n, d, e )
359 IF( start.LE.n )
THEN
369 IF( finish.LT.n )
THEN
370 tiny = eps*sqrt( abs( d( finish ) ) )*
371 $ sqrt( abs( d( finish+1 ) ) )
372 IF( abs( e( finish ) ).GT.tiny )
THEN
380 m = finish - start + 1
385 IF( m.GT.smlsiz )
THEN
389 orgnrm =
dlanst(
'M', m, d( start ), e( start ) )
390 CALL
dlascl(
'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
392 CALL
dlascl(
'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
395 IF( icompz.EQ.1 )
THEN
400 CALL
dlaed0( icompz, n, m, d( start ), e( start ),
401 $ z( strtrw, start ), ldz, work( 1 ), n,
402 $ work( storez ), iwork, info )
404 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
405 $ mod( info, ( m+1 ) ) + start - 1
411 CALL
dlascl(
'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
415 IF( icompz.EQ.1 )
THEN
421 CALL
dsteqr(
'I', m, d( start ), e( start ), work, m,
422 $ work( m*m+1 ), info )
423 CALL
dlacpy(
'A', n, m, z( 1, start ), ldz,
424 $ work( storez ), n )
425 CALL
dgemm(
'N',
'N', n, m, m, one,
426 $ work( storez ), n, work, m, zero,
427 $ z( 1, start ), ldz )
428 ELSE IF( icompz.EQ.2 )
THEN
429 CALL
dsteqr(
'I', m, d( start ), e( start ),
430 $ z( start, start ), ldz, work, info )
432 CALL
dsterf( m, d( start ), e( start ), info )
435 info = start*( n+1 ) + finish
451 IF( icompz.EQ.0 )
THEN
455 CALL
dlasrt(
'I', n, d, info )
466 IF( d(
j ).LT.p )
THEN
474 CALL
dswap( n, z( 1, i ), 1, z( 1, k ), 1 )
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
LOGICAL function lsame(CA, CB)
LSAME
subroutine dlaed0(ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, WORK, IWORK, INFO)
DLAED0 used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dsterf(N, D, E, INFO)
DSTERF
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.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEBZ
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.
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...