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 #ifndef TEMPLATE_LAPACK_SYEV_HEADER
00036 #define TEMPLATE_LAPACK_SYEV_HEADER
00037
00038
00039 template<class Treal>
00040 int template_lapack_syev(const char *jobz, const char *uplo, const integer *n, Treal *a,
00041 const integer *lda, Treal *w, Treal *work, const integer *lwork,
00042 integer *info)
00043 {
00044
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 integer c__1 = 1;
00119 integer c_n1 = -1;
00120 integer c__0 = 0;
00121 Treal c_b17 = 1.;
00122
00123
00124 integer a_dim1, a_offset, i__1, i__2;
00125 Treal d__1;
00126
00127 integer inde;
00128 Treal anrm;
00129 integer imax;
00130 Treal rmin, rmax;
00131 Treal sigma;
00132 integer iinfo;
00133 logical lower, wantz;
00134 integer nb;
00135 integer iscale;
00136 Treal safmin;
00137 Treal bignum;
00138 integer indtau;
00139 integer indwrk;
00140 integer llwork;
00141 Treal smlnum;
00142 integer lwkopt;
00143 logical lquery;
00144 Treal eps;
00145 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00146
00147
00148 a_dim1 = *lda;
00149 a_offset = 1 + a_dim1 * 1;
00150 a -= a_offset;
00151 --w;
00152 --work;
00153
00154
00155 lwkopt = 0;
00156
00157 wantz = template_blas_lsame(jobz, "V");
00158 lower = template_blas_lsame(uplo, "L");
00159 lquery = *lwork == -1;
00160
00161 *info = 0;
00162 if (! (wantz || template_blas_lsame(jobz, "N"))) {
00163 *info = -1;
00164 } else if (! (lower || template_blas_lsame(uplo, "U"))) {
00165 *info = -2;
00166 } else if (*n < 0) {
00167 *info = -3;
00168 } else if (*lda < maxMACRO(1,*n)) {
00169 *info = -5;
00170 } else {
00171
00172 i__1 = 1, i__2 = *n * 3 - 1;
00173 if (*lwork < maxMACRO(i__1,i__2) && ! lquery) {
00174 *info = -8;
00175 }
00176 }
00177
00178 if (*info == 0) {
00179 nb = template_lapack_ilaenv(&c__1, "DSYTRD", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6,
00180 (ftnlen)1);
00181
00182 i__1 = 1, i__2 = (nb + 2) * *n;
00183 lwkopt = maxMACRO(i__1,i__2);
00184 work[1] = (Treal) lwkopt;
00185 }
00186
00187 if (*info != 0) {
00188 i__1 = -(*info);
00189 template_blas_erbla("SYEV ", &i__1);
00190 return 0;
00191 } else if (lquery) {
00192 return 0;
00193 }
00194
00195
00196
00197 if (*n == 0) {
00198 work[1] = 1.;
00199 return 0;
00200 }
00201
00202 if (*n == 1) {
00203 w[1] = a_ref(1, 1);
00204 work[1] = 3.;
00205 if (wantz) {
00206 a_ref(1, 1) = 1.;
00207 }
00208 return 0;
00209 }
00210
00211
00212
00213
00214
00215
00216 safmin = template_lapack_lamch("Safe minimum", (Treal)0);
00217
00218
00219 eps = template_lapack_lamch("Precision", (Treal)0);
00220
00221
00222
00223 smlnum = safmin / eps;
00224 bignum = 1. / smlnum;
00225 rmin = template_blas_sqrt(smlnum);
00226 rmax = template_blas_sqrt(bignum);
00227
00228
00229
00230 anrm = template_lapack_lansy("M", uplo, n, &a[a_offset], lda, &work[1]);
00231 iscale = 0;
00232 if (anrm > 0. && anrm < rmin) {
00233 iscale = 1;
00234 sigma = rmin / anrm;
00235 } else if (anrm > rmax) {
00236 iscale = 1;
00237 sigma = rmax / anrm;
00238 }
00239 if (iscale == 1) {
00240 template_lapack_lascl(uplo, &c__0, &c__0, &c_b17, &sigma, n, n, &a[a_offset], lda,
00241 info);
00242 }
00243
00244
00245
00246 inde = 1;
00247 indtau = inde + *n;
00248 indwrk = indtau + *n;
00249 llwork = *lwork - indwrk + 1;
00250 template_lapack_sytrd(uplo, n, &a[a_offset], lda, &w[1], &work[inde], &work[indtau], &
00251 work[indwrk], &llwork, &iinfo);
00252
00253
00254
00255
00256 if (! wantz) {
00257 template_lapack_sterf(n, &w[1], &work[inde], info);
00258 } else {
00259 template_lapack_orgtr(uplo, n, &a[a_offset], lda, &work[indtau], &work[indwrk], &
00260 llwork, &iinfo);
00261 template_lapack_steqr(jobz, n, &w[1], &work[inde], &a[a_offset], lda, &work[indtau],
00262 info);
00263 }
00264
00265
00266
00267 if (iscale == 1) {
00268 if (*info == 0) {
00269 imax = *n;
00270 } else {
00271 imax = *info - 1;
00272 }
00273 d__1 = 1. / sigma;
00274 template_blas_scal(&imax, &d__1, &w[1], &c__1);
00275 }
00276
00277
00278
00279 work[1] = (Treal) lwkopt;
00280
00281 return 0;
00282
00283
00284
00285 }
00286
00287 #undef a_ref
00288
00289
00290 #endif