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