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