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