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_STEVX_HEADER
00036 #define TEMPLATE_LAPACK_STEVX_HEADER
00037
00038
00039 template<class Treal>
00040 int template_lapack_stevx(const char *jobz, const char *range, const integer *n, Treal *
00041 d__, Treal *e, const Treal *vl, const Treal *vu, const integer *il,
00042 const integer *iu, const Treal *abstol, integer *m, Treal *w,
00043 Treal *z__, const integer *ldz, Treal *work, integer *iwork,
00044 integer *ifail, 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
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172 integer c__1 = 1;
00173
00174
00175 integer z_dim1, z_offset, i__1, i__2;
00176 Treal d__1, d__2;
00177
00178 integer imax;
00179 Treal rmin, rmax, tnrm;
00180 integer itmp1, i__, j;
00181 Treal sigma;
00182 char order[1];
00183 logical wantz;
00184 integer jj;
00185 logical alleig, indeig;
00186 integer iscale, indibl;
00187 logical valeig;
00188 Treal safmin;
00189 Treal bignum;
00190 integer indisp;
00191 integer indiwo;
00192 integer indwrk;
00193 integer nsplit;
00194 Treal smlnum, eps, vll, vuu, tmp1;
00195 #define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
00196
00197
00198 --d__;
00199 --e;
00200 --w;
00201 z_dim1 = *ldz;
00202 z_offset = 1 + z_dim1 * 1;
00203 z__ -= z_offset;
00204 --work;
00205 --iwork;
00206 --ifail;
00207
00208
00209 wantz = template_blas_lsame(jobz, "V");
00210 alleig = template_blas_lsame(range, "A");
00211 valeig = template_blas_lsame(range, "V");
00212 indeig = template_blas_lsame(range, "I");
00213
00214 *info = 0;
00215 if (! (wantz || template_blas_lsame(jobz, "N"))) {
00216 *info = -1;
00217 } else if (! (alleig || valeig || indeig)) {
00218 *info = -2;
00219 } else if (*n < 0) {
00220 *info = -3;
00221 } else {
00222 if (valeig) {
00223 if (*n > 0 && *vu <= *vl) {
00224 *info = -7;
00225 }
00226 } else if (indeig) {
00227 if (*il < 1 || *il > maxMACRO(1,*n)) {
00228 *info = -8;
00229 } else if (*iu < minMACRO(*n,*il) || *iu > *n) {
00230 *info = -9;
00231 }
00232 }
00233 }
00234 if (*info == 0) {
00235 if (*ldz < 1 || wantz && *ldz < *n) {
00236 *info = -14;
00237 }
00238 }
00239
00240 if (*info != 0) {
00241 i__1 = -(*info);
00242 template_blas_erbla("STEVX ", &i__1);
00243 return 0;
00244 }
00245
00246
00247
00248 *m = 0;
00249 if (*n == 0) {
00250 return 0;
00251 }
00252
00253 if (*n == 1) {
00254 if (alleig || indeig) {
00255 *m = 1;
00256 w[1] = d__[1];
00257 } else {
00258 if (*vl < d__[1] && *vu >= d__[1]) {
00259 *m = 1;
00260 w[1] = d__[1];
00261 }
00262 }
00263 if (wantz) {
00264 z___ref(1, 1) = 1.;
00265 }
00266 return 0;
00267 }
00268
00269
00270
00271 safmin = template_lapack_lamch("Safe minimum", (Treal)0);
00272 eps = template_lapack_lamch("Precision", (Treal)0);
00273 smlnum = safmin / eps;
00274 bignum = 1. / smlnum;
00275 rmin = template_blas_sqrt(smlnum);
00276
00277 d__1 = template_blas_sqrt(bignum), d__2 = 1. / template_blas_sqrt(template_blas_sqrt(safmin));
00278 rmax = minMACRO(d__1,d__2);
00279
00280
00281
00282 iscale = 0;
00283 if (valeig) {
00284 vll = *vl;
00285 vuu = *vu;
00286 } else {
00287 vll = 0.;
00288 vuu = 0.;
00289 }
00290 tnrm = template_lapack_lanst("M", n, &d__[1], &e[1]);
00291 if (tnrm > 0. && tnrm < rmin) {
00292 iscale = 1;
00293 sigma = rmin / tnrm;
00294 } else if (tnrm > rmax) {
00295 iscale = 1;
00296 sigma = rmax / tnrm;
00297 }
00298 if (iscale == 1) {
00299 template_blas_scal(n, &sigma, &d__[1], &c__1);
00300 i__1 = *n - 1;
00301 template_blas_scal(&i__1, &sigma, &e[1], &c__1);
00302 if (valeig) {
00303 vll = *vl * sigma;
00304 vuu = *vu * sigma;
00305 }
00306 }
00307
00308
00309
00310
00311
00312 if ((alleig || indeig && *il == 1 && *iu == *n) && *abstol <= 0.) {
00313 template_blas_copy(n, &d__[1], &c__1, &w[1], &c__1);
00314 i__1 = *n - 1;
00315 template_blas_copy(&i__1, &e[1], &c__1, &work[1], &c__1);
00316 indwrk = *n + 1;
00317 if (! wantz) {
00318 template_lapack_sterf(n, &w[1], &work[1], info);
00319 } else {
00320 template_lapack_steqr("I", n, &w[1], &work[1], &z__[z_offset], ldz, &work[
00321 indwrk], info);
00322 if (*info == 0) {
00323 i__1 = *n;
00324 for (i__ = 1; i__ <= i__1; ++i__) {
00325 ifail[i__] = 0;
00326
00327 }
00328 }
00329 }
00330 if (*info == 0) {
00331 *m = *n;
00332 goto L20;
00333 }
00334 *info = 0;
00335 }
00336
00337
00338
00339 if (wantz) {
00340 *(unsigned char *)order = 'B';
00341 } else {
00342 *(unsigned char *)order = 'E';
00343 }
00344 indwrk = 1;
00345 indibl = 1;
00346 indisp = indibl + *n;
00347 indiwo = indisp + *n;
00348 template_lapack_stebz(range, order, n, &vll, &vuu, il, iu, abstol, &d__[1], &e[1], m, &
00349 nsplit, &w[1], &iwork[indibl], &iwork[indisp], &work[indwrk], &
00350 iwork[indiwo], info);
00351
00352 if (wantz) {
00353 template_lapack_stein(n, &d__[1], &e[1], m, &w[1], &iwork[indibl], &iwork[indisp], &
00354 z__[z_offset], ldz, &work[indwrk], &iwork[indiwo], &ifail[1],
00355 info);
00356 }
00357
00358
00359
00360 L20:
00361 if (iscale == 1) {
00362 if (*info == 0) {
00363 imax = *m;
00364 } else {
00365 imax = *info - 1;
00366 }
00367 d__1 = 1. / sigma;
00368 template_blas_scal(&imax, &d__1, &w[1], &c__1);
00369 }
00370
00371
00372
00373
00374 if (wantz) {
00375 i__1 = *m - 1;
00376 for (j = 1; j <= i__1; ++j) {
00377 i__ = 0;
00378 tmp1 = w[j];
00379 i__2 = *m;
00380 for (jj = j + 1; jj <= i__2; ++jj) {
00381 if (w[jj] < tmp1) {
00382 i__ = jj;
00383 tmp1 = w[jj];
00384 }
00385
00386 }
00387
00388 if (i__ != 0) {
00389 itmp1 = iwork[indibl + i__ - 1];
00390 w[i__] = w[j];
00391 iwork[indibl + i__ - 1] = iwork[indibl + j - 1];
00392 w[j] = tmp1;
00393 iwork[indibl + j - 1] = itmp1;
00394 template_blas_swap(n, &z___ref(1, i__), &c__1, &z___ref(1, j), &c__1);
00395 if (*info != 0) {
00396 itmp1 = ifail[i__];
00397 ifail[i__] = ifail[j];
00398 ifail[j] = itmp1;
00399 }
00400 }
00401
00402 }
00403 }
00404
00405 return 0;
00406
00407
00408
00409 }
00410
00411 #undef z___ref
00412
00413
00414 #endif