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