LAPACK: SRC/dsteqr.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          COMPZ
00010       INTEGER            INFO, LDZ, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
00020 *  symmetric tridiagonal matrix using the implicit QL or QR method.
00021 *  The eigenvectors of a full or band symmetric matrix can also be found
00022 *  if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
00023 *  tridiagonal form.
00024 *
00025 *  Arguments
00026 *  =========
00027 *
00028 *  COMPZ   (input) CHARACTER*1
00029 *          = 'N':  Compute eigenvalues only.
00030 *          = 'V':  Compute eigenvalues and eigenvectors of the original
00031 *                  symmetric matrix.  On entry, Z must contain the
00032 *                  orthogonal matrix used to reduce the original matrix
00033 *                  to tridiagonal form.
00034 *          = 'I':  Compute eigenvalues and eigenvectors of the
00035 *                  tridiagonal matrix.  Z is initialized to the identity
00036 *                  matrix.
00037 *
00038 *  N       (input) INTEGER
00039 *          The order of the matrix.  N >= 0.
00040 *
00041 *  D       (input/output) DOUBLE PRECISION array, dimension (N)
00042 *          On entry, the diagonal elements of the tridiagonal matrix.
00043 *          On exit, if INFO = 0, the eigenvalues in ascending order.
00044 *
00045 *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
00046 *          On entry, the (n-1) subdiagonal elements of the tridiagonal
00047 *          matrix.
00048 *          On exit, E has been destroyed.
00049 *
00050 *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
00051 *          On entry, if  COMPZ = 'V', then Z contains the orthogonal
00052 *          matrix used in the reduction to tridiagonal form.
00053 *          On exit, if INFO = 0, then if  COMPZ = 'V', Z contains the
00054 *          orthonormal eigenvectors of the original symmetric matrix,
00055 *          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
00056 *          of the symmetric tridiagonal matrix.
00057 *          If COMPZ = 'N', then Z is not referenced.
00058 *
00059 *  LDZ     (input) INTEGER
00060 *          The leading dimension of the array Z.  LDZ >= 1, and if
00061 *          eigenvectors are desired, then  LDZ >= max(1,N).
00062 *
00063 *  WORK    (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
00064 *          If COMPZ = 'N', then WORK is not referenced.
00065 *
00066 *  INFO    (output) INTEGER
00067 *          = 0:  successful exit
00068 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00069 *          > 0:  the algorithm has failed to find all the eigenvalues in
00070 *                a total of 30*N iterations; if INFO = i, then i
00071 *                elements of E have not converged to zero; on exit, D
00072 *                and E contain the elements of a symmetric tridiagonal
00073 *                matrix which is orthogonally similar to the original
00074 *                matrix.
00075 *
00076 *  =====================================================================
00077 *
00078 *     .. Parameters ..
00079       DOUBLE PRECISION   ZERO, ONE, TWO, THREE
00080       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
00081      $                   THREE = 3.0D0 )
00082       INTEGER            MAXIT
00083       PARAMETER          ( MAXIT = 30 )
00084 *     ..
00085 *     .. Local Scalars ..
00086       INTEGER            I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
00087      $                   LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
00088      $                   NM1, NMAXIT
00089       DOUBLE PRECISION   ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
00090      $                   S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
00091 *     ..
00092 *     .. External Functions ..
00093       LOGICAL            LSAME
00094       DOUBLE PRECISION   DLAMCH, DLANST, DLAPY2
00095       EXTERNAL           LSAME, DLAMCH, DLANST, DLAPY2
00096 *     ..
00097 *     .. External Subroutines ..
00098       EXTERNAL           DLAE2, DLAEV2, DLARTG, DLASCL, DLASET, DLASR,
00099      $                   DLASRT, DSWAP, XERBLA
00100 *     ..
00101 *     .. Intrinsic Functions ..
00102       INTRINSIC          ABS, MAX, SIGN, SQRT
00103 *     ..
00104 *     .. Executable Statements ..
00105 *
00106 *     Test the input parameters.
00107 *
00108       INFO = 0
00109 *
00110       IF( LSAME( COMPZ, 'N' ) ) THEN
00111          ICOMPZ = 0
00112       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
00113          ICOMPZ = 1
00114       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
00115          ICOMPZ = 2
00116       ELSE
00117          ICOMPZ = -1
00118       END IF
00119       IF( ICOMPZ.LT.0 ) THEN
00120          INFO = -1
00121       ELSE IF( N.LT.0 ) THEN
00122          INFO = -2
00123       ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
00124      $         N ) ) ) THEN
00125          INFO = -6
00126       END IF
00127       IF( INFO.NE.0 ) THEN
00128          CALL XERBLA( 'DSTEQR', -INFO )
00129          RETURN
00130       END IF
00131 *
00132 *     Quick return if possible
00133 *
00134       IF( N.EQ.0 )
00135      $   RETURN
00136 *
00137       IF( N.EQ.1 ) THEN
00138          IF( ICOMPZ.EQ.2 )
00139      $      Z( 1, 1 ) = ONE
00140          RETURN
00141       END IF
00142 *
00143 *     Determine the unit roundoff and over/underflow thresholds.
00144 *
00145       EPS = DLAMCH( 'E' )
00146       EPS2 = EPS**2
00147       SAFMIN = DLAMCH( 'S' )
00148       SAFMAX = ONE / SAFMIN
00149       SSFMAX = SQRT( SAFMAX ) / THREE
00150       SSFMIN = SQRT( SAFMIN ) / EPS2
00151 *
00152 *     Compute the eigenvalues and eigenvectors of the tridiagonal
00153 *     matrix.
00154 *
00155       IF( ICOMPZ.EQ.2 )
00156      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
00157 *
00158       NMAXIT = N*MAXIT
00159       JTOT = 0
00160 *
00161 *     Determine where the matrix splits and choose QL or QR iteration
00162 *     for each block, according to whether top or bottom diagonal
00163 *     element is smaller.
00164 *
00165       L1 = 1
00166       NM1 = N - 1
00167 *
00168    10 CONTINUE
00169       IF( L1.GT.N )
00170      $   GO TO 160
00171       IF( L1.GT.1 )
00172      $   E( L1-1 ) = ZERO
00173       IF( L1.LE.NM1 ) THEN
00174          DO 20 M = L1, NM1
00175             TST = ABS( E( M ) )
00176             IF( TST.EQ.ZERO )
00177      $         GO TO 30
00178             IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
00179      $          1 ) ) ) )*EPS ) THEN
00180                E( M ) = ZERO
00181                GO TO 30
00182             END IF
00183    20    CONTINUE
00184       END IF
00185       M = N
00186 *
00187    30 CONTINUE
00188       L = L1
00189       LSV = L
00190       LEND = M
00191       LENDSV = LEND
00192       L1 = M + 1
00193       IF( LEND.EQ.L )
00194      $   GO TO 10
00195 *
00196 *     Scale submatrix in rows and columns L to LEND
00197 *
00198       ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) )
00199       ISCALE = 0
00200       IF( ANORM.EQ.ZERO )
00201      $   GO TO 10
00202       IF( ANORM.GT.SSFMAX ) THEN
00203          ISCALE = 1
00204          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
00205      $                INFO )
00206          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
00207      $                INFO )
00208       ELSE IF( ANORM.LT.SSFMIN ) THEN
00209          ISCALE = 2
00210          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
00211      $                INFO )
00212          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
00213      $                INFO )
00214       END IF
00215 *
00216 *     Choose between QL and QR iteration
00217 *
00218       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
00219          LEND = LSV
00220          L = LENDSV
00221       END IF
00222 *
00223       IF( LEND.GT.L ) THEN
00224 *
00225 *        QL Iteration
00226 *
00227 *        Look for small subdiagonal element.
00228 *
00229    40    CONTINUE
00230          IF( L.NE.LEND ) THEN
00231             LENDM1 = LEND - 1
00232             DO 50 M = L, LENDM1
00233                TST = ABS( E( M ) )**2
00234                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
00235      $             SAFMIN )GO TO 60
00236    50       CONTINUE
00237          END IF
00238 *
00239          M = LEND
00240 *
00241    60    CONTINUE
00242          IF( M.LT.LEND )
00243      $      E( M ) = ZERO
00244          P = D( L )
00245          IF( M.EQ.L )
00246      $      GO TO 80
00247 *
00248 *        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
00249 *        to compute its eigensystem.
00250 *
00251          IF( M.EQ.L+1 ) THEN
00252             IF( ICOMPZ.GT.0 ) THEN
00253                CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
00254                WORK( L ) = C
00255                WORK( N-1+L ) = S
00256                CALL DLASR( 'R', 'V', 'B', N, 2, WORK( L ),
00257      $                     WORK( N-1+L ), Z( 1, L ), LDZ )
00258             ELSE
00259                CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 )
00260             END IF
00261             D( L ) = RT1
00262             D( L+1 ) = RT2
00263             E( L ) = ZERO
00264             L = L + 2
00265             IF( L.LE.LEND )
00266      $         GO TO 40
00267             GO TO 140
00268          END IF
00269 *
00270          IF( JTOT.EQ.NMAXIT )
00271      $      GO TO 140
00272          JTOT = JTOT + 1
00273 *
00274 *        Form shift.
00275 *
00276          G = ( D( L+1 )-P ) / ( TWO*E( L ) )
00277          R = DLAPY2( G, ONE )
00278          G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
00279 *
00280          S = ONE
00281          C = ONE
00282          P = ZERO
00283 *
00284 *        Inner loop
00285 *
00286          MM1 = M - 1
00287          DO 70 I = MM1, L, -1
00288             F = S*E( I )
00289             B = C*E( I )
00290             CALL DLARTG( G, F, C, S, R )
00291             IF( I.NE.M-1 )
00292      $         E( I+1 ) = R
00293             G = D( I+1 ) - P
00294             R = ( D( I )-G )*S + TWO*C*B
00295             P = S*R
00296             D( I+1 ) = G + P
00297             G = C*R - B
00298 *
00299 *           If eigenvectors are desired, then save rotations.
00300 *
00301             IF( ICOMPZ.GT.0 ) THEN
00302                WORK( I ) = C
00303                WORK( N-1+I ) = -S
00304             END IF
00305 *
00306    70    CONTINUE
00307 *
00308 *        If eigenvectors are desired, then apply saved rotations.
00309 *
00310          IF( ICOMPZ.GT.0 ) THEN
00311             MM = M - L + 1
00312             CALL DLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ),
00313      $                  Z( 1, L ), LDZ )
00314          END IF
00315 *
00316          D( L ) = D( L ) - P
00317          E( L ) = G
00318          GO TO 40
00319 *
00320 *        Eigenvalue found.
00321 *
00322    80    CONTINUE
00323          D( L ) = P
00324 *
00325          L = L + 1
00326          IF( L.LE.LEND )
00327      $      GO TO 40
00328          GO TO 140
00329 *
00330       ELSE
00331 *
00332 *        QR Iteration
00333 *
00334 *        Look for small superdiagonal element.
00335 *
00336    90    CONTINUE
00337          IF( L.NE.LEND ) THEN
00338             LENDP1 = LEND + 1
00339             DO 100 M = L, LENDP1, -1
00340                TST = ABS( E( M-1 ) )**2
00341                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
00342      $             SAFMIN )GO TO 110
00343   100       CONTINUE
00344          END IF
00345 *
00346          M = LEND
00347 *
00348   110    CONTINUE
00349          IF( M.GT.LEND )
00350      $      E( M-1 ) = ZERO
00351          P = D( L )
00352          IF( M.EQ.L )
00353      $      GO TO 130
00354 *
00355 *        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
00356 *        to compute its eigensystem.
00357 *
00358          IF( M.EQ.L-1 ) THEN
00359             IF( ICOMPZ.GT.0 ) THEN
00360                CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
00361                WORK( M ) = C
00362                WORK( N-1+M ) = S
00363                CALL DLASR( 'R', 'V', 'F', N, 2, WORK( M ),
00364      $                     WORK( N-1+M ), Z( 1, L-1 ), LDZ )
00365             ELSE
00366                CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
00367             END IF
00368             D( L-1 ) = RT1
00369             D( L ) = RT2
00370             E( L-1 ) = ZERO
00371             L = L - 2
00372             IF( L.GE.LEND )
00373      $         GO TO 90
00374             GO TO 140
00375          END IF
00376 *
00377          IF( JTOT.EQ.NMAXIT )
00378      $      GO TO 140
00379          JTOT = JTOT + 1
00380 *
00381 *        Form shift.
00382 *
00383          G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
00384          R = DLAPY2( G, ONE )
00385          G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
00386 *
00387          S = ONE
00388          C = ONE
00389          P = ZERO
00390 *
00391 *        Inner loop
00392 *
00393          LM1 = L - 1
00394          DO 120 I = M, LM1
00395             F = S*E( I )
00396             B = C*E( I )
00397             CALL DLARTG( G, F, C, S, R )
00398             IF( I.NE.M )
00399      $         E( I-1 ) = R
00400             G = D( I ) - P
00401             R = ( D( I+1 )-G )*S + TWO*C*B
00402             P = S*R
00403             D( I ) = G + P
00404             G = C*R - B
00405 *
00406 *           If eigenvectors are desired, then save rotations.
00407 *
00408             IF( ICOMPZ.GT.0 ) THEN
00409                WORK( I ) = C
00410                WORK( N-1+I ) = S
00411             END IF
00412 *
00413   120    CONTINUE
00414 *
00415 *        If eigenvectors are desired, then apply saved rotations.
00416 *
00417          IF( ICOMPZ.GT.0 ) THEN
00418             MM = L - M + 1
00419             CALL DLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ),
00420      $                  Z( 1, M ), LDZ )
00421          END IF
00422 *
00423          D( L ) = D( L ) - P
00424          E( LM1 ) = G
00425          GO TO 90
00426 *
00427 *        Eigenvalue found.
00428 *
00429   130    CONTINUE
00430          D( L ) = P
00431 *
00432          L = L - 1
00433          IF( L.GE.LEND )
00434      $      GO TO 90
00435          GO TO 140
00436 *
00437       END IF
00438 *
00439 *     Undo scaling if necessary
00440 *
00441   140 CONTINUE
00442       IF( ISCALE.EQ.1 ) THEN
00443          CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
00444      $                D( LSV ), N, INFO )
00445          CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
00446      $                N, INFO )
00447       ELSE IF( ISCALE.EQ.2 ) THEN
00448          CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
00449      $                D( LSV ), N, INFO )
00450          CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
00451      $                N, INFO )
00452       END IF
00453 *
00454 *     Check for no convergence to an eigenvalue after a total
00455 *     of N*MAXIT iterations.
00456 *
00457       IF( JTOT.LT.NMAXIT )
00458      $   GO TO 10
00459       DO 150 I = 1, N - 1
00460          IF( E( I ).NE.ZERO )
00461      $      INFO = INFO + 1
00462   150 CONTINUE
00463       GO TO 190
00464 *
00465 *     Order eigenvalues and eigenvectors.
00466 *
00467   160 CONTINUE
00468       IF( ICOMPZ.EQ.0 ) THEN
00469 *
00470 *        Use Quick Sort
00471 *
00472          CALL DLASRT( 'I', N, D, INFO )
00473 *
00474       ELSE
00475 *
00476 *        Use Selection Sort to minimize swaps of eigenvectors
00477 *
00478          DO 180 II = 2, N
00479             I = II - 1
00480             K = I
00481             P = D( I )
00482             DO 170 J = II, N
00483                IF( D( J ).LT.P ) THEN
00484                   K = J
00485                   P = D( J )
00486                END IF
00487   170       CONTINUE
00488             IF( K.NE.I ) THEN
00489                D( K ) = D( I )
00490                D( I ) = P
00491                CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
00492             END IF
00493   180    CONTINUE
00494       END IF
00495 *
00496   190 CONTINUE
00497       RETURN
00498 *
00499 *     End of DSTEQR
00500 *
00501       END