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