LAPACK: SRC/dstedc.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
00002      $                   LIWORK, INFO )
00003 *
00004 *  -- LAPACK driver routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          COMPZ
00011       INTEGER            INFO, LDZ, LIWORK, LWORK, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            IWORK( * )
00015       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
00022 *  symmetric tridiagonal matrix using the divide and conquer method.
00023 *  The eigenvectors of a full or band real symmetric matrix can also be
00024 *  found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
00025 *  matrix to tridiagonal form.
00026 *
00027 *  This code makes very mild assumptions about floating point
00028 *  arithmetic. It will work on machines with a guard digit in
00029 *  add/subtract, or on those binary machines without guard digits
00030 *  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
00031 *  It could conceivably fail on hexadecimal or decimal machines
00032 *  without guard digits, but we know of none.  See DLAED3 for details.
00033 *
00034 *  Arguments
00035 *  =========
00036 *
00037 *  COMPZ   (input) CHARACTER*1
00038 *          = 'N':  Compute eigenvalues only.
00039 *          = 'I':  Compute eigenvectors of tridiagonal matrix also.
00040 *          = 'V':  Compute eigenvectors of original dense symmetric
00041 *                  matrix also.  On entry, Z contains the orthogonal
00042 *                  matrix used to reduce the original matrix to
00043 *                  tridiagonal form.
00044 *
00045 *  N       (input) INTEGER
00046 *          The dimension of the symmetric tridiagonal matrix.  N >= 0.
00047 *
00048 *  D       (input/output) DOUBLE PRECISION array, dimension (N)
00049 *          On entry, the diagonal elements of the tridiagonal matrix.
00050 *          On exit, if INFO = 0, the eigenvalues in ascending order.
00051 *
00052 *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
00053 *          On entry, the subdiagonal elements of the tridiagonal matrix.
00054 *          On exit, E has been destroyed.
00055 *
00056 *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
00057 *          On entry, if COMPZ = 'V', then Z contains the orthogonal
00058 *          matrix used in the reduction to tridiagonal form.
00059 *          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
00060 *          orthonormal eigenvectors of the original symmetric matrix,
00061 *          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
00062 *          of the symmetric tridiagonal matrix.
00063 *          If  COMPZ = 'N', then Z is not referenced.
00064 *
00065 *  LDZ     (input) INTEGER
00066 *          The leading dimension of the array Z.  LDZ >= 1.
00067 *          If eigenvectors are desired, then LDZ >= max(1,N).
00068 *
00069 *  WORK    (workspace/output) DOUBLE PRECISION array,
00070 *                                         dimension (LWORK)
00071 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00072 *
00073 *  LWORK   (input) INTEGER
00074 *          The dimension of the array WORK.
00075 *          If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
00076 *          If COMPZ = 'V' and N > 1 then LWORK must be at least
00077 *                         ( 1 + 3*N + 2*N*lg N + 3*N**2 ),
00078 *                         where lg( N ) = smallest integer k such
00079 *                         that 2**k >= N.
00080 *          If COMPZ = 'I' and N > 1 then LWORK must be at least
00081 *                         ( 1 + 4*N + N**2 ).
00082 *          Note that for COMPZ = 'I' or 'V', then if N is less than or
00083 *          equal to the minimum divide size, usually 25, then LWORK need
00084 *          only be max(1,2*(N-1)).
00085 *
00086 *          If LWORK = -1, then a workspace query is assumed; the routine
00087 *          only calculates the optimal size of the WORK array, returns
00088 *          this value as the first entry of the WORK array, and no error
00089 *          message related to LWORK is issued by XERBLA.
00090 *
00091 *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
00092 *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
00093 *
00094 *  LIWORK  (input) INTEGER
00095 *          The dimension of the array IWORK.
00096 *          If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
00097 *          If COMPZ = 'V' and N > 1 then LIWORK must be at least
00098 *                         ( 6 + 6*N + 5*N*lg N ).
00099 *          If COMPZ = 'I' and N > 1 then LIWORK must be at least
00100 *                         ( 3 + 5*N ).
00101 *          Note that for COMPZ = 'I' or 'V', then if N is less than or
00102 *          equal to the minimum divide size, usually 25, then LIWORK
00103 *          need only be 1.
00104 *
00105 *          If LIWORK = -1, then a workspace query is assumed; the
00106 *          routine only calculates the optimal size of the IWORK array,
00107 *          returns this value as the first entry of the IWORK array, and
00108 *          no error message related to LIWORK is issued by XERBLA.
00109 *
00110 *  INFO    (output) INTEGER
00111 *          = 0:  successful exit.
00112 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00113 *          > 0:  The algorithm failed to compute an eigenvalue while
00114 *                working on the submatrix lying in rows and columns
00115 *                INFO/(N+1) through mod(INFO,N+1).
00116 *
00117 *  Further Details
00118 *  ===============
00119 *
00120 *  Based on contributions by
00121 *     Jeff Rutter, Computer Science Division, University of California
00122 *     at Berkeley, USA
00123 *  Modified by Francoise Tisseur, University of Tennessee.
00124 *
00125 *  =====================================================================
00126 *
00127 *     .. Parameters ..
00128       DOUBLE PRECISION   ZERO, ONE, TWO
00129       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
00130 *     ..
00131 *     .. Local Scalars ..
00132       LOGICAL            LQUERY
00133       INTEGER            FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
00134      $                   LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
00135       DOUBLE PRECISION   EPS, ORGNRM, P, TINY
00136 *     ..
00137 *     .. External Functions ..
00138       LOGICAL            LSAME
00139       INTEGER            ILAENV
00140       DOUBLE PRECISION   DLAMCH, DLANST
00141       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANST
00142 *     ..
00143 *     .. External Subroutines ..
00144       EXTERNAL           DGEMM, DLACPY, DLAED0, DLASCL, DLASET, DLASRT,
00145      $                   DSTEQR, DSTERF, DSWAP, XERBLA
00146 *     ..
00147 *     .. Intrinsic Functions ..
00148       INTRINSIC          ABS, DBLE, INT, LOG, MAX, MOD, SQRT
00149 *     ..
00150 *     .. Executable Statements ..
00151 *
00152 *     Test the input parameters.
00153 *
00154       INFO = 0
00155       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
00156 *
00157       IF( LSAME( COMPZ, 'N' ) ) THEN
00158          ICOMPZ = 0
00159       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
00160          ICOMPZ = 1
00161       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
00162          ICOMPZ = 2
00163       ELSE
00164          ICOMPZ = -1
00165       END IF
00166       IF( ICOMPZ.LT.0 ) THEN
00167          INFO = -1
00168       ELSE IF( N.LT.0 ) THEN
00169          INFO = -2
00170       ELSE IF( ( LDZ.LT.1 ) .OR.
00171      $         ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
00172          INFO = -6
00173       END IF
00174 *
00175       IF( INFO.EQ.0 ) THEN
00176 *
00177 *        Compute the workspace requirements
00178 *
00179          SMLSIZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
00180          IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
00181             LIWMIN = 1
00182             LWMIN = 1
00183          ELSE IF( N.LE.SMLSIZ ) THEN
00184             LIWMIN = 1
00185             LWMIN = 2*( N - 1 )
00186          ELSE
00187             LGN = INT( LOG( DBLE( N ) )/LOG( TWO ) )
00188             IF( 2**LGN.LT.N )
00189      $         LGN = LGN + 1
00190             IF( 2**LGN.LT.N )
00191      $         LGN = LGN + 1
00192             IF( ICOMPZ.EQ.1 ) THEN
00193                LWMIN = 1 + 3*N + 2*N*LGN + 3*N**2
00194                LIWMIN = 6 + 6*N + 5*N*LGN
00195             ELSE IF( ICOMPZ.EQ.2 ) THEN
00196                LWMIN = 1 + 4*N + N**2
00197                LIWMIN = 3 + 5*N
00198             END IF
00199          END IF
00200          WORK( 1 ) = LWMIN
00201          IWORK( 1 ) = LIWMIN
00202 *
00203          IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
00204             INFO = -8
00205          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
00206             INFO = -10
00207          END IF
00208       END IF
00209 *
00210       IF( INFO.NE.0 ) THEN
00211          CALL XERBLA( 'DSTEDC', -INFO )
00212          RETURN
00213       ELSE IF (LQUERY) THEN
00214          RETURN
00215       END IF
00216 *
00217 *     Quick return if possible
00218 *
00219       IF( N.EQ.0 )
00220      $   RETURN
00221       IF( N.EQ.1 ) THEN
00222          IF( ICOMPZ.NE.0 )
00223      $      Z( 1, 1 ) = ONE
00224          RETURN
00225       END IF
00226 *
00227 *     If the following conditional clause is removed, then the routine
00228 *     will use the Divide and Conquer routine to compute only the
00229 *     eigenvalues, which requires (3N + 3N**2) real workspace and
00230 *     (2 + 5N + 2N lg(N)) integer workspace.
00231 *     Since on many architectures DSTERF is much faster than any other
00232 *     algorithm for finding eigenvalues only, it is used here
00233 *     as the default. If the conditional clause is removed, then
00234 *     information on the size of workspace needs to be changed.
00235 *
00236 *     If COMPZ = 'N', use DSTERF to compute the eigenvalues.
00237 *
00238       IF( ICOMPZ.EQ.0 ) THEN
00239          CALL DSTERF( N, D, E, INFO )
00240          GO TO 50
00241       END IF
00242 *
00243 *     If N is smaller than the minimum divide size (SMLSIZ+1), then
00244 *     solve the problem with another solver.
00245 *
00246       IF( N.LE.SMLSIZ ) THEN
00247 *
00248          CALL DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
00249 *
00250       ELSE
00251 *
00252 *        If COMPZ = 'V', the Z matrix must be stored elsewhere for later
00253 *        use.
00254 *
00255          IF( ICOMPZ.EQ.1 ) THEN
00256             STOREZ = 1 + N*N
00257          ELSE
00258             STOREZ = 1
00259          END IF
00260 *
00261          IF( ICOMPZ.EQ.2 ) THEN
00262             CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
00263          END IF
00264 *
00265 *        Scale.
00266 *
00267          ORGNRM = DLANST( 'M', N, D, E )
00268          IF( ORGNRM.EQ.ZERO )
00269      $      GO TO 50
00270 *
00271          EPS = DLAMCH( 'Epsilon' )
00272 *
00273          START = 1
00274 *
00275 *        while ( START <= N )
00276 *
00277    10    CONTINUE
00278          IF( START.LE.N ) THEN
00279 *
00280 *           Let FINISH be the position of the next subdiagonal entry
00281 *           such that E( FINISH ) <= TINY or FINISH = N if no such
00282 *           subdiagonal exists.  The matrix identified by the elements
00283 *           between START and FINISH constitutes an independent
00284 *           sub-problem.
00285 *
00286             FINISH = START
00287    20       CONTINUE
00288             IF( FINISH.LT.N ) THEN
00289                TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
00290      $                    SQRT( ABS( D( FINISH+1 ) ) )
00291                IF( ABS( E( FINISH ) ).GT.TINY ) THEN
00292                   FINISH = FINISH + 1
00293                   GO TO 20
00294                END IF
00295             END IF
00296 *
00297 *           (Sub) Problem determined.  Compute its size and solve it.
00298 *
00299             M = FINISH - START + 1
00300             IF( M.EQ.1 ) THEN
00301                START = FINISH + 1
00302                GO TO 10
00303             END IF
00304             IF( M.GT.SMLSIZ ) THEN
00305 *
00306 *              Scale.
00307 *
00308                ORGNRM = DLANST( 'M', M, D( START ), E( START ) )
00309                CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
00310      $                      INFO )
00311                CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
00312      $                      M-1, INFO )
00313 *
00314                IF( ICOMPZ.EQ.1 ) THEN
00315                   STRTRW = 1
00316                ELSE
00317                   STRTRW = START
00318                END IF
00319                CALL DLAED0( ICOMPZ, N, M, D( START ), E( START ),
00320      $                      Z( STRTRW, START ), LDZ, WORK( 1 ), N,
00321      $                      WORK( STOREZ ), IWORK, INFO )
00322                IF( INFO.NE.0 ) THEN
00323                   INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
00324      $                   MOD( INFO, ( M+1 ) ) + START - 1
00325                   GO TO 50
00326                END IF
00327 *
00328 *              Scale back.
00329 *
00330                CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
00331      $                      INFO )
00332 *
00333             ELSE
00334                IF( ICOMPZ.EQ.1 ) THEN
00335 *
00336 *                 Since QR won't update a Z matrix which is larger than
00337 *                 the length of D, we must solve the sub-problem in a
00338 *                 workspace and then multiply back into Z.
00339 *
00340                   CALL DSTEQR( 'I', M, D( START ), E( START ), WORK, M,
00341      $                         WORK( M*M+1 ), INFO )
00342                   CALL DLACPY( 'A', N, M, Z( 1, START ), LDZ,
00343      $                         WORK( STOREZ ), N )
00344                   CALL DGEMM( 'N', 'N', N, M, M, ONE,
00345      $                        WORK( STOREZ ), N, WORK, M, ZERO,
00346      $                        Z( 1, START ), LDZ )
00347                ELSE IF( ICOMPZ.EQ.2 ) THEN
00348                   CALL DSTEQR( 'I', M, D( START ), E( START ),
00349      $                         Z( START, START ), LDZ, WORK, INFO )
00350                ELSE
00351                   CALL DSTERF( M, D( START ), E( START ), INFO )
00352                END IF
00353                IF( INFO.NE.0 ) THEN
00354                   INFO = START*( N+1 ) + FINISH
00355                   GO TO 50
00356                END IF
00357             END IF
00358 *
00359             START = FINISH + 1
00360             GO TO 10
00361          END IF
00362 *
00363 *        endwhile
00364 *
00365 *        If the problem split any number of times, then the eigenvalues
00366 *        will not be properly ordered.  Here we permute the eigenvalues
00367 *        (and the associated eigenvectors) into ascending order.
00368 *
00369          IF( M.NE.N ) THEN
00370             IF( ICOMPZ.EQ.0 ) THEN
00371 *
00372 *              Use Quick Sort
00373 *
00374                CALL DLASRT( 'I', N, D, INFO )
00375 *
00376             ELSE
00377 *
00378 *              Use Selection Sort to minimize swaps of eigenvectors
00379 *
00380                DO 40 II = 2, N
00381                   I = II - 1
00382                   K = I
00383                   P = D( I )
00384                   DO 30 J = II, N
00385                      IF( D( J ).LT.P ) THEN
00386                         K = J
00387                         P = D( J )
00388                      END IF
00389    30             CONTINUE
00390                   IF( K.NE.I ) THEN
00391                      D( K ) = D( I )
00392                      D( I ) = P
00393                      CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
00394                   END IF
00395    40          CONTINUE
00396             END IF
00397          END IF
00398       END IF
00399 *
00400    50 CONTINUE
00401       WORK( 1 ) = LWMIN
00402       IWORK( 1 ) = LIWMIN
00403 *
00404       RETURN
00405 *
00406 *     End of DSTEDC
00407 *
00408       END