00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #ifndef TEMPLATE_LAPACK_SYEV_HEADER
00038 #define TEMPLATE_LAPACK_SYEV_HEADER
00039
00040
00041 template<class Treal>
00042 int template_lapack_syev(const char *jobz, const char *uplo, const integer *n, Treal *a,
00043 const integer *lda, Treal *w, Treal *work, const integer *lwork,
00044 integer *info)
00045 {
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120 integer c__1 = 1;
00121 integer c_n1 = -1;
00122 integer c__0 = 0;
00123 Treal c_b17 = 1.;
00124
00125
00126 integer a_dim1, a_offset, i__1, i__2;
00127 Treal d__1;
00128
00129 integer inde;
00130 Treal anrm;
00131 integer imax;
00132 Treal rmin, rmax;
00133 Treal sigma;
00134 integer iinfo;
00135 logical lower, wantz;
00136 integer nb;
00137 integer iscale;
00138 Treal safmin;
00139 Treal bignum;
00140 integer indtau;
00141 integer indwrk;
00142 integer llwork;
00143 Treal smlnum;
00144 integer lwkopt;
00145 logical lquery;
00146 Treal eps;
00147 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00148
00149
00150 a_dim1 = *lda;
00151 a_offset = 1 + a_dim1 * 1;
00152 a -= a_offset;
00153 --w;
00154 --work;
00155
00156
00157 lwkopt = 0;
00158
00159 wantz = template_blas_lsame(jobz, "V");
00160 lower = template_blas_lsame(uplo, "L");
00161 lquery = *lwork == -1;
00162
00163 *info = 0;
00164 if (! (wantz || template_blas_lsame(jobz, "N"))) {
00165 *info = -1;
00166 } else if (! (lower || template_blas_lsame(uplo, "U"))) {
00167 *info = -2;
00168 } else if (*n < 0) {
00169 *info = -3;
00170 } else if (*lda < maxMACRO(1,*n)) {
00171 *info = -5;
00172 } else {
00173
00174 i__1 = 1, i__2 = *n * 3 - 1;
00175 if (*lwork < maxMACRO(i__1,i__2) && ! lquery) {
00176 *info = -8;
00177 }
00178 }
00179
00180 if (*info == 0) {
00181 nb = template_lapack_ilaenv(&c__1, "DSYTRD", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6,
00182 (ftnlen)1);
00183
00184 i__1 = 1, i__2 = (nb + 2) * *n;
00185 lwkopt = maxMACRO(i__1,i__2);
00186 work[1] = (Treal) lwkopt;
00187 }
00188
00189 if (*info != 0) {
00190 i__1 = -(*info);
00191 template_blas_erbla("SYEV ", &i__1);
00192 return 0;
00193 } else if (lquery) {
00194 return 0;
00195 }
00196
00197
00198
00199 if (*n == 0) {
00200 work[1] = 1.;
00201 return 0;
00202 }
00203
00204 if (*n == 1) {
00205 w[1] = a_ref(1, 1);
00206 work[1] = 3.;
00207 if (wantz) {
00208 a_ref(1, 1) = 1.;
00209 }
00210 return 0;
00211 }
00212
00213
00214
00215
00216
00217
00218 safmin = template_lapack_lamch("Safe minimum", (Treal)0);
00219
00220
00221 eps = template_lapack_lamch("Precision", (Treal)0);
00222
00223
00224
00225 smlnum = safmin / eps;
00226 bignum = 1. / smlnum;
00227 rmin = template_blas_sqrt(smlnum);
00228 rmax = template_blas_sqrt(bignum);
00229
00230
00231
00232 anrm = template_lapack_lansy("M", uplo, n, &a[a_offset], lda, &work[1]);
00233 iscale = 0;
00234 if (anrm > 0. && anrm < rmin) {
00235 iscale = 1;
00236 sigma = rmin / anrm;
00237 } else if (anrm > rmax) {
00238 iscale = 1;
00239 sigma = rmax / anrm;
00240 }
00241 if (iscale == 1) {
00242 template_lapack_lascl(uplo, &c__0, &c__0, &c_b17, &sigma, n, n, &a[a_offset], lda,
00243 info);
00244 }
00245
00246
00247
00248 inde = 1;
00249 indtau = inde + *n;
00250 indwrk = indtau + *n;
00251 llwork = *lwork - indwrk + 1;
00252 template_lapack_sytrd(uplo, n, &a[a_offset], lda, &w[1], &work[inde], &work[indtau], &
00253 work[indwrk], &llwork, &iinfo);
00254
00255
00256
00257
00258 if (! wantz) {
00259 template_lapack_sterf(n, &w[1], &work[inde], info);
00260 } else {
00261 template_lapack_orgtr(uplo, n, &a[a_offset], lda, &work[indtau], &work[indwrk], &
00262 llwork, &iinfo);
00263 template_lapack_steqr(jobz, n, &w[1], &work[inde], &a[a_offset], lda, &work[indtau],
00264 info);
00265 }
00266
00267
00268
00269 if (iscale == 1) {
00270 if (*info == 0) {
00271 imax = *n;
00272 } else {
00273 imax = *info - 1;
00274 }
00275 d__1 = 1. / sigma;
00276 template_blas_scal(&imax, &d__1, &w[1], &c__1);
00277 }
00278
00279
00280
00281 work[1] = (Treal) lwkopt;
00282
00283 return 0;
00284
00285
00286
00287 }
00288
00289 #undef a_ref
00290
00291
00292 #endif