00001 /* sspevd.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static integer c__1 = 1; 00019 00020 /* Subroutine */ int sspevd_(char *jobz, char *uplo, integer *n, real *ap, 00021 real *w, real *z__, integer *ldz, real *work, integer *lwork, integer 00022 *iwork, integer *liwork, integer *info) 00023 { 00024 /* System generated locals */ 00025 integer z_dim1, z_offset, i__1; 00026 real r__1; 00027 00028 /* Builtin functions */ 00029 double sqrt(doublereal); 00030 00031 /* Local variables */ 00032 real eps; 00033 integer inde; 00034 real anrm, rmin, rmax, sigma; 00035 extern logical lsame_(char *, char *); 00036 integer iinfo; 00037 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *); 00038 integer lwmin; 00039 logical wantz; 00040 integer iscale; 00041 extern doublereal slamch_(char *); 00042 real safmin; 00043 extern /* Subroutine */ int xerbla_(char *, integer *); 00044 real bignum; 00045 integer indtau; 00046 extern /* Subroutine */ int sstedc_(char *, integer *, real *, real *, 00047 real *, integer *, real *, integer *, integer *, integer *, 00048 integer *); 00049 integer indwrk, liwmin; 00050 extern doublereal slansp_(char *, char *, integer *, real *, real *); 00051 extern /* Subroutine */ int ssterf_(integer *, real *, real *, integer *); 00052 integer llwork; 00053 real smlnum; 00054 extern /* Subroutine */ int ssptrd_(char *, integer *, real *, real *, 00055 real *, real *, integer *); 00056 logical lquery; 00057 extern /* Subroutine */ int sopmtr_(char *, char *, char *, integer *, 00058 integer *, real *, real *, real *, integer *, real *, integer *); 00059 00060 00061 /* -- LAPACK driver routine (version 3.2) -- */ 00062 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00063 /* November 2006 */ 00064 00065 /* .. Scalar Arguments .. */ 00066 /* .. */ 00067 /* .. Array Arguments .. */ 00068 /* .. */ 00069 00070 /* Purpose */ 00071 /* ======= */ 00072 00073 /* SSPEVD computes all the eigenvalues and, optionally, eigenvectors */ 00074 /* of a real symmetric matrix A in packed storage. If eigenvectors are */ 00075 /* desired, it uses a divide and conquer algorithm. */ 00076 00077 /* The divide and conquer algorithm makes very mild assumptions about */ 00078 /* floating point arithmetic. It will work on machines with a guard */ 00079 /* digit in add/subtract, or on those binary machines without guard */ 00080 /* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or */ 00081 /* Cray-2. It could conceivably fail on hexadecimal or decimal machines */ 00082 /* without guard digits, but we know of none. */ 00083 00084 /* Arguments */ 00085 /* ========= */ 00086 00087 /* JOBZ (input) CHARACTER*1 */ 00088 /* = 'N': Compute eigenvalues only; */ 00089 /* = 'V': Compute eigenvalues and eigenvectors. */ 00090 00091 /* UPLO (input) CHARACTER*1 */ 00092 /* = 'U': Upper triangle of A is stored; */ 00093 /* = 'L': Lower triangle of A is stored. */ 00094 00095 /* N (input) INTEGER */ 00096 /* The order of the matrix A. N >= 0. */ 00097 00098 /* AP (input/output) REAL array, dimension (N*(N+1)/2) */ 00099 /* On entry, the upper or lower triangle of the symmetric matrix */ 00100 /* A, packed columnwise in a linear array. The j-th column of A */ 00101 /* is stored in the array AP as follows: */ 00102 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */ 00103 /* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. */ 00104 00105 /* On exit, AP is overwritten by values generated during the */ 00106 /* reduction to tridiagonal form. If UPLO = 'U', the diagonal */ 00107 /* and first superdiagonal of the tridiagonal matrix T overwrite */ 00108 /* the corresponding elements of A, and if UPLO = 'L', the */ 00109 /* diagonal and first subdiagonal of T overwrite the */ 00110 /* corresponding elements of A. */ 00111 00112 /* W (output) REAL array, dimension (N) */ 00113 /* If INFO = 0, the eigenvalues in ascending order. */ 00114 00115 /* Z (output) REAL array, dimension (LDZ, N) */ 00116 /* If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal */ 00117 /* eigenvectors of the matrix A, with the i-th column of Z */ 00118 /* holding the eigenvector associated with W(i). */ 00119 /* If JOBZ = 'N', then Z is not referenced. */ 00120 00121 /* LDZ (input) INTEGER */ 00122 /* The leading dimension of the array Z. LDZ >= 1, and if */ 00123 /* JOBZ = 'V', LDZ >= max(1,N). */ 00124 00125 /* WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) */ 00126 /* On exit, if INFO = 0, WORK(1) returns the required LWORK. */ 00127 00128 /* LWORK (input) INTEGER */ 00129 /* The dimension of the array WORK. */ 00130 /* If N <= 1, LWORK must be at least 1. */ 00131 /* If JOBZ = 'N' and N > 1, LWORK must be at least 2*N. */ 00132 /* If JOBZ = 'V' and N > 1, LWORK must be at least */ 00133 /* 1 + 6*N + N**2. */ 00134 00135 /* If LWORK = -1, then a workspace query is assumed; the routine */ 00136 /* only calculates the required sizes of the WORK and IWORK */ 00137 /* arrays, returns these values as the first entries of the WORK */ 00138 /* and IWORK arrays, and no error message related to LWORK or */ 00139 /* LIWORK is issued by XERBLA. */ 00140 00141 /* IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) */ 00142 /* On exit, if INFO = 0, IWORK(1) returns the required LIWORK. */ 00143 00144 /* LIWORK (input) INTEGER */ 00145 /* The dimension of the array IWORK. */ 00146 /* If JOBZ = 'N' or N <= 1, LIWORK must be at least 1. */ 00147 /* If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N. */ 00148 00149 /* If LIWORK = -1, then a workspace query is assumed; the */ 00150 /* routine only calculates the required sizes of the WORK and */ 00151 /* IWORK arrays, returns these values as the first entries of */ 00152 /* the WORK and IWORK arrays, and no error message related to */ 00153 /* LWORK or LIWORK is issued by XERBLA. */ 00154 00155 /* INFO (output) INTEGER */ 00156 /* = 0: successful exit */ 00157 /* < 0: if INFO = -i, the i-th argument had an illegal value. */ 00158 /* > 0: if INFO = i, the algorithm failed to converge; i */ 00159 /* off-diagonal elements of an intermediate tridiagonal */ 00160 /* form did not converge to zero. */ 00161 00162 /* ===================================================================== */ 00163 00164 /* .. Parameters .. */ 00165 /* .. */ 00166 /* .. Local Scalars .. */ 00167 /* .. */ 00168 /* .. External Functions .. */ 00169 /* .. */ 00170 /* .. External Subroutines .. */ 00171 /* .. */ 00172 /* .. Intrinsic Functions .. */ 00173 /* .. */ 00174 /* .. Executable Statements .. */ 00175 00176 /* Test the input parameters. */ 00177 00178 /* Parameter adjustments */ 00179 --ap; 00180 --w; 00181 z_dim1 = *ldz; 00182 z_offset = 1 + z_dim1; 00183 z__ -= z_offset; 00184 --work; 00185 --iwork; 00186 00187 /* Function Body */ 00188 wantz = lsame_(jobz, "V"); 00189 lquery = *lwork == -1 || *liwork == -1; 00190 00191 *info = 0; 00192 if (! (wantz || lsame_(jobz, "N"))) { 00193 *info = -1; 00194 } else if (! (lsame_(uplo, "U") || lsame_(uplo, 00195 "L"))) { 00196 *info = -2; 00197 } else if (*n < 0) { 00198 *info = -3; 00199 } else if (*ldz < 1 || wantz && *ldz < *n) { 00200 *info = -7; 00201 } 00202 00203 if (*info == 0) { 00204 if (*n <= 1) { 00205 liwmin = 1; 00206 lwmin = 1; 00207 } else { 00208 if (wantz) { 00209 liwmin = *n * 5 + 3; 00210 /* Computing 2nd power */ 00211 i__1 = *n; 00212 lwmin = *n * 6 + 1 + i__1 * i__1; 00213 } else { 00214 liwmin = 1; 00215 lwmin = *n << 1; 00216 } 00217 } 00218 iwork[1] = liwmin; 00219 work[1] = (real) lwmin; 00220 00221 if (*lwork < lwmin && ! lquery) { 00222 *info = -9; 00223 } else if (*liwork < liwmin && ! lquery) { 00224 *info = -11; 00225 } 00226 } 00227 00228 if (*info != 0) { 00229 i__1 = -(*info); 00230 xerbla_("SSPEVD", &i__1); 00231 return 0; 00232 } else if (lquery) { 00233 return 0; 00234 } 00235 00236 /* Quick return if possible */ 00237 00238 if (*n == 0) { 00239 return 0; 00240 } 00241 00242 if (*n == 1) { 00243 w[1] = ap[1]; 00244 if (wantz) { 00245 z__[z_dim1 + 1] = 1.f; 00246 } 00247 return 0; 00248 } 00249 00250 /* Get machine constants. */ 00251 00252 safmin = slamch_("Safe minimum"); 00253 eps = slamch_("Precision"); 00254 smlnum = safmin / eps; 00255 bignum = 1.f / smlnum; 00256 rmin = sqrt(smlnum); 00257 rmax = sqrt(bignum); 00258 00259 /* Scale matrix to allowable range, if necessary. */ 00260 00261 anrm = slansp_("M", uplo, n, &ap[1], &work[1]); 00262 iscale = 0; 00263 if (anrm > 0.f && anrm < rmin) { 00264 iscale = 1; 00265 sigma = rmin / anrm; 00266 } else if (anrm > rmax) { 00267 iscale = 1; 00268 sigma = rmax / anrm; 00269 } 00270 if (iscale == 1) { 00271 i__1 = *n * (*n + 1) / 2; 00272 sscal_(&i__1, &sigma, &ap[1], &c__1); 00273 } 00274 00275 /* Call SSPTRD to reduce symmetric packed matrix to tridiagonal form. */ 00276 00277 inde = 1; 00278 indtau = inde + *n; 00279 ssptrd_(uplo, n, &ap[1], &w[1], &work[inde], &work[indtau], &iinfo); 00280 00281 /* For eigenvalues only, call SSTERF. For eigenvectors, first call */ 00282 /* SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the */ 00283 /* tridiagonal matrix, then call SOPMTR to multiply it by the */ 00284 /* Householder transformations represented in AP. */ 00285 00286 if (! wantz) { 00287 ssterf_(n, &w[1], &work[inde], info); 00288 } else { 00289 indwrk = indtau + *n; 00290 llwork = *lwork - indwrk + 1; 00291 sstedc_("I", n, &w[1], &work[inde], &z__[z_offset], ldz, &work[indwrk] 00292 , &llwork, &iwork[1], liwork, info); 00293 sopmtr_("L", uplo, "N", n, n, &ap[1], &work[indtau], &z__[z_offset], 00294 ldz, &work[indwrk], &iinfo); 00295 } 00296 00297 /* If matrix was scaled, then rescale eigenvalues appropriately. */ 00298 00299 if (iscale == 1) { 00300 r__1 = 1.f / sigma; 00301 sscal_(n, &r__1, &w[1], &c__1); 00302 } 00303 00304 work[1] = (real) lwmin; 00305 iwork[1] = liwmin; 00306 return 0; 00307 00308 /* End of SSPEVD */ 00309 00310 } /* sspevd_ */