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