00001 SUBROUTINE SSTEVD( 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 REAL D( * ), E( * ), WORK( * ), Z( LDZ, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * SSTEVD 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) REAL 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) REAL 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) REAL 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) REAL 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 REAL ZERO, ONE 00102 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00103 * .. 00104 * .. Local Scalars .. 00105 LOGICAL LQUERY, WANTZ 00106 INTEGER ISCALE, LIWMIN, LWMIN 00107 REAL BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM, 00108 $ TNRM 00109 * .. 00110 * .. External Functions .. 00111 LOGICAL LSAME 00112 REAL SLAMCH, SLANST 00113 EXTERNAL LSAME, SLAMCH, SLANST 00114 * .. 00115 * .. External Subroutines .. 00116 EXTERNAL SSCAL, SSTEDC, SSTERF, 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( 'SSTEVD', -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 = SLAMCH( 'Safe minimum' ) 00176 EPS = SLAMCH( '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 = SLANST( '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 SSCAL( N, SIGMA, D, 1 ) 00195 CALL SSCAL( N-1, SIGMA, E( 1 ), 1 ) 00196 END IF 00197 * 00198 * For eigenvalues only, call SSTERF. For eigenvalues and 00199 * eigenvectors, call SSTEDC. 00200 * 00201 IF( .NOT.WANTZ ) THEN 00202 CALL SSTERF( N, D, E, INFO ) 00203 ELSE 00204 CALL SSTEDC( '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 SSCAL( N, ONE / SIGMA, D, 1 ) 00212 * 00213 WORK( 1 ) = LWMIN 00214 IWORK( 1 ) = LIWMIN 00215 * 00216 RETURN 00217 * 00218 * End of SSTEVD 00219 * 00220 END
1.7.2