LAPACK: SRC/csteqr.f Source File

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