LAPACK 3.3.1
Linear Algebra PACKage

dstevd.f

Go to the documentation of this file.
00001       SUBROUTINE DSTEVD( JOBZ, 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          JOBZ
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 *  DSTEVD computes all eigenvalues and, optionally, eigenvectors of a
00022 *  real symmetric tridiagonal matrix. If eigenvectors are desired, it
00023 *  uses a divide and conquer algorithm.
00024 *
00025 *  The divide and conquer algorithm makes very mild assumptions about
00026 *  floating point arithmetic. It will work on machines with a guard
00027 *  digit in add/subtract, or on those binary machines without guard
00028 *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
00029 *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
00030 *  without guard digits, but we know of none.
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  JOBZ    (input) CHARACTER*1
00036 *          = 'N':  Compute eigenvalues only;
00037 *          = 'V':  Compute eigenvalues and eigenvectors.
00038 *
00039 *  N       (input) INTEGER
00040 *          The order of the matrix.  N >= 0.
00041 *
00042 *  D       (input/output) DOUBLE PRECISION array, dimension (N)
00043 *          On entry, the n diagonal elements of the tridiagonal matrix
00044 *          A.
00045 *          On exit, if INFO = 0, the eigenvalues in ascending order.
00046 *
00047 *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
00048 *          On entry, the (n-1) subdiagonal elements of the tridiagonal
00049 *          matrix A, stored in elements 1 to N-1 of E.
00050 *          On exit, the contents of E are destroyed.
00051 *
00052 *  Z       (output) DOUBLE PRECISION array, dimension (LDZ, N)
00053 *          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
00054 *          eigenvectors of the matrix A, with the i-th column of Z
00055 *          holding the eigenvector associated with D(i).
00056 *          If JOBZ = 'N', then Z is not referenced.
00057 *
00058 *  LDZ     (input) INTEGER
00059 *          The leading dimension of the array Z.  LDZ >= 1, and if
00060 *          JOBZ = 'V', LDZ >= max(1,N).
00061 *
00062 *  WORK    (workspace/output) DOUBLE PRECISION array,
00063 *                                         dimension (LWORK)
00064 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00065 *
00066 *  LWORK   (input) INTEGER
00067 *          The dimension of the array WORK.
00068 *          If JOBZ  = 'N' or N <= 1 then LWORK must be at least 1.
00069 *          If JOBZ  = 'V' and N > 1 then LWORK must be at least
00070 *                         ( 1 + 4*N + N**2 ).
00071 *
00072 *          If LWORK = -1, then a workspace query is assumed; the routine
00073 *          only calculates the optimal sizes of the WORK and IWORK
00074 *          arrays, returns these values as the first entries of the WORK
00075 *          and IWORK arrays, and no error message related to LWORK or
00076 *          LIWORK is issued by XERBLA.
00077 *
00078 *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
00079 *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
00080 *
00081 *  LIWORK  (input) INTEGER
00082 *          The dimension of the array IWORK.
00083 *          If JOBZ  = 'N' or N <= 1 then LIWORK must be at least 1.
00084 *          If JOBZ  = 'V' and N > 1 then LIWORK must be at least 3+5*N.
00085 *
00086 *          If LIWORK = -1, then a workspace query is assumed; the
00087 *          routine only calculates the optimal sizes of the WORK and
00088 *          IWORK arrays, returns these values as the first entries of
00089 *          the WORK and IWORK arrays, and no error message related to
00090 *          LWORK or LIWORK is issued by XERBLA.
00091 *
00092 *  INFO    (output) INTEGER
00093 *          = 0:  successful exit
00094 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00095 *          > 0:  if INFO = i, the algorithm failed to converge; i
00096 *                off-diagonal elements of E did not converge to zero.
00097 *
00098 *  =====================================================================
00099 *
00100 *     .. Parameters ..
00101       DOUBLE PRECISION   ZERO, ONE
00102       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00103 *     ..
00104 *     .. Local Scalars ..
00105       LOGICAL            LQUERY, WANTZ
00106       INTEGER            ISCALE, LIWMIN, LWMIN
00107       DOUBLE PRECISION   BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
00108      $                   TNRM
00109 *     ..
00110 *     .. External Functions ..
00111       LOGICAL            LSAME
00112       DOUBLE PRECISION   DLAMCH, DLANST
00113       EXTERNAL           LSAME, DLAMCH, DLANST
00114 *     ..
00115 *     .. External Subroutines ..
00116       EXTERNAL           DSCAL, DSTEDC, DSTERF, XERBLA
00117 *     ..
00118 *     .. Intrinsic Functions ..
00119       INTRINSIC          SQRT
00120 *     ..
00121 *     .. Executable Statements ..
00122 *
00123 *     Test the input parameters.
00124 *
00125       WANTZ = LSAME( JOBZ, 'V' )
00126       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
00127 *
00128       INFO = 0
00129       LIWMIN = 1
00130       LWMIN = 1
00131       IF( N.GT.1 .AND. WANTZ ) THEN
00132          LWMIN = 1 + 4*N + N**2
00133          LIWMIN = 3 + 5*N
00134       END IF
00135 *
00136       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
00137          INFO = -1
00138       ELSE IF( N.LT.0 ) THEN
00139          INFO = -2
00140       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
00141          INFO = -6
00142       END IF
00143 *
00144       IF( INFO.EQ.0 ) THEN
00145          WORK( 1 ) = LWMIN
00146          IWORK( 1 ) = LIWMIN
00147 *
00148          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00149             INFO = -8
00150          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
00151             INFO = -10
00152          END IF
00153       END IF
00154 *
00155       IF( INFO.NE.0 ) THEN
00156          CALL XERBLA( 'DSTEVD', -INFO )
00157          RETURN
00158       ELSE IF( LQUERY ) THEN
00159          RETURN
00160       END IF
00161 *
00162 *     Quick return if possible
00163 *
00164       IF( N.EQ.0 )
00165      $   RETURN
00166 *
00167       IF( N.EQ.1 ) THEN
00168          IF( WANTZ )
00169      $      Z( 1, 1 ) = ONE
00170          RETURN
00171       END IF
00172 *
00173 *     Get machine constants.
00174 *
00175       SAFMIN = DLAMCH( 'Safe minimum' )
00176       EPS = DLAMCH( 'Precision' )
00177       SMLNUM = SAFMIN / EPS
00178       BIGNUM = ONE / SMLNUM
00179       RMIN = SQRT( SMLNUM )
00180       RMAX = SQRT( BIGNUM )
00181 *
00182 *     Scale matrix to allowable range, if necessary.
00183 *
00184       ISCALE = 0
00185       TNRM = DLANST( 'M', N, D, E )
00186       IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
00187          ISCALE = 1
00188          SIGMA = RMIN / TNRM
00189       ELSE IF( TNRM.GT.RMAX ) THEN
00190          ISCALE = 1
00191          SIGMA = RMAX / TNRM
00192       END IF
00193       IF( ISCALE.EQ.1 ) THEN
00194          CALL DSCAL( N, SIGMA, D, 1 )
00195          CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
00196       END IF
00197 *
00198 *     For eigenvalues only, call DSTERF.  For eigenvalues and
00199 *     eigenvectors, call DSTEDC.
00200 *
00201       IF( .NOT.WANTZ ) THEN
00202          CALL DSTERF( N, D, E, INFO )
00203       ELSE
00204          CALL DSTEDC( 'I', N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK,
00205      $                INFO )
00206       END IF
00207 *
00208 *     If matrix was scaled, then rescale eigenvalues appropriately.
00209 *
00210       IF( ISCALE.EQ.1 )
00211      $   CALL DSCAL( N, ONE / SIGMA, D, 1 )
00212 *
00213       WORK( 1 ) = LWMIN
00214       IWORK( 1 ) = LIWMIN
00215 *
00216       RETURN
00217 *
00218 *     End of DSTEVD
00219 *
00220       END
 All Files Functions