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_STEQR_HEADER
00038 #define TEMPLATE_LAPACK_STEQR_HEADER
00039
00040
00041 template<class Treal>
00042 int template_lapack_steqr(const char *compz, const integer *n, Treal *d__,
00043 Treal *e, Treal *z__, const integer *ldz, Treal *work,
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 Treal c_b9 = 0.;
00120 Treal c_b10 = 1.;
00121 integer c__0 = 0;
00122 integer c__1 = 1;
00123 integer c__2 = 2;
00124
00125
00126 integer z_dim1, z_offset, i__1, i__2;
00127 Treal d__1, d__2;
00128
00129 integer lend, jtot;
00130 Treal b, c__, f, g;
00131 integer i__, j, k, l, m;
00132 Treal p, r__, s;
00133 Treal anorm;
00134 integer l1;
00135 integer lendm1, lendp1;
00136 integer ii;
00137 integer mm, iscale;
00138 Treal safmin;
00139 Treal safmax;
00140 integer lendsv;
00141 Treal ssfmin;
00142 integer nmaxit, icompz;
00143 Treal ssfmax;
00144 integer lm1, mm1, nm1;
00145 Treal rt1, rt2, eps;
00146 integer lsv;
00147 Treal tst, eps2;
00148 #define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
00149
00150
00151 --d__;
00152 --e;
00153 z_dim1 = *ldz;
00154 z_offset = 1 + z_dim1 * 1;
00155 z__ -= z_offset;
00156 --work;
00157
00158
00159 *info = 0;
00160
00161 if (template_blas_lsame(compz, "N")) {
00162 icompz = 0;
00163 } else if (template_blas_lsame(compz, "V")) {
00164 icompz = 1;
00165 } else if (template_blas_lsame(compz, "I")) {
00166 icompz = 2;
00167 } else {
00168 icompz = -1;
00169 }
00170 if (icompz < 0) {
00171 *info = -1;
00172 } else if (*n < 0) {
00173 *info = -2;
00174 } else if (*ldz < 1 || (icompz > 0 && *ldz < maxMACRO(1,*n) ) ) {
00175 *info = -6;
00176 }
00177 if (*info != 0) {
00178 i__1 = -(*info);
00179 template_blas_erbla("STEQR ", &i__1);
00180 return 0;
00181 }
00182
00183
00184
00185 if (*n == 0) {
00186 return 0;
00187 }
00188
00189 if (*n == 1) {
00190 if (icompz == 2) {
00191 z___ref(1, 1) = 1.;
00192 }
00193 return 0;
00194 }
00195
00196
00197
00198 eps = template_lapack_lamch("E", (Treal)0);
00199
00200 d__1 = eps;
00201 eps2 = d__1 * d__1;
00202 safmin = template_lapack_lamch("S", (Treal)0);
00203 safmax = 1. / safmin;
00204 ssfmax = template_blas_sqrt(safmax) / 3.;
00205 ssfmin = template_blas_sqrt(safmin) / eps2;
00206
00207
00208
00209
00210 if (icompz == 2) {
00211 template_lapack_laset("Full", n, n, &c_b9, &c_b10, &z__[z_offset], ldz);
00212 }
00213
00214 nmaxit = *n * 30;
00215 jtot = 0;
00216
00217
00218
00219
00220
00221 l1 = 1;
00222 nm1 = *n - 1;
00223
00224 L10:
00225 if (l1 > *n) {
00226 goto L160;
00227 }
00228 if (l1 > 1) {
00229 e[l1 - 1] = 0.;
00230 }
00231 if (l1 <= nm1) {
00232 i__1 = nm1;
00233 for (m = l1; m <= i__1; ++m) {
00234 tst = (d__1 = e[m], absMACRO(d__1));
00235 if (tst == 0.) {
00236 goto L30;
00237 }
00238 if (tst <= template_blas_sqrt((d__1 = d__[m], absMACRO(d__1))) * template_blas_sqrt((d__2 = d__[m
00239 + 1], absMACRO(d__2))) * eps) {
00240 e[m] = 0.;
00241 goto L30;
00242 }
00243
00244 }
00245 }
00246 m = *n;
00247
00248 L30:
00249 l = l1;
00250 lsv = l;
00251 lend = m;
00252 lendsv = lend;
00253 l1 = m + 1;
00254 if (lend == l) {
00255 goto L10;
00256 }
00257
00258
00259
00260 i__1 = lend - l + 1;
00261 anorm = template_lapack_lanst("I", &i__1, &d__[l], &e[l]);
00262 iscale = 0;
00263 if (anorm == 0.) {
00264 goto L10;
00265 }
00266 if (anorm > ssfmax) {
00267 iscale = 1;
00268 i__1 = lend - l + 1;
00269 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
00270 info);
00271 i__1 = lend - l;
00272 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
00273 info);
00274 } else if (anorm < ssfmin) {
00275 iscale = 2;
00276 i__1 = lend - l + 1;
00277 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
00278 info);
00279 i__1 = lend - l;
00280 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
00281 info);
00282 }
00283
00284
00285
00286 if ((d__1 = d__[lend], absMACRO(d__1)) < (d__2 = d__[l], absMACRO(d__2))) {
00287 lend = lsv;
00288 l = lendsv;
00289 }
00290
00291 if (lend > l) {
00292
00293
00294
00295
00296
00297 L40:
00298 if (l != lend) {
00299 lendm1 = lend - 1;
00300 i__1 = lendm1;
00301 for (m = l; m <= i__1; ++m) {
00302
00303 d__2 = (d__1 = e[m], absMACRO(d__1));
00304 tst = d__2 * d__2;
00305 if (tst <= eps2 * (d__1 = d__[m], absMACRO(d__1)) * (d__2 = d__[m
00306 + 1], absMACRO(d__2)) + safmin) {
00307 goto L60;
00308 }
00309
00310 }
00311 }
00312
00313 m = lend;
00314
00315 L60:
00316 if (m < lend) {
00317 e[m] = 0.;
00318 }
00319 p = d__[l];
00320 if (m == l) {
00321 goto L80;
00322 }
00323
00324
00325
00326
00327 if (m == l + 1) {
00328 if (icompz > 0) {
00329 template_lapack_laev2(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
00330 work[l] = c__;
00331 work[*n - 1 + l] = s;
00332 template_lapack_lasr("R", "V", "B", n, &c__2, &work[l], &work[*n - 1 + l], &
00333 z___ref(1, l), ldz);
00334 } else {
00335 template_lapack_lae2(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
00336 }
00337 d__[l] = rt1;
00338 d__[l + 1] = rt2;
00339 e[l] = 0.;
00340 l += 2;
00341 if (l <= lend) {
00342 goto L40;
00343 }
00344 goto L140;
00345 }
00346
00347 if (jtot == nmaxit) {
00348 goto L140;
00349 }
00350 ++jtot;
00351
00352
00353
00354 g = (d__[l + 1] - p) / (e[l] * 2.);
00355 r__ = template_lapack_lapy2(&g, &c_b10);
00356 g = d__[m] - p + e[l] / (g + template_lapack_d_sign(&r__, &g));
00357
00358 s = 1.;
00359 c__ = 1.;
00360 p = 0.;
00361
00362
00363
00364 mm1 = m - 1;
00365 i__1 = l;
00366 for (i__ = mm1; i__ >= i__1; --i__) {
00367 f = s * e[i__];
00368 b = c__ * e[i__];
00369 template_lapack_lartg(&g, &f, &c__, &s, &r__);
00370 if (i__ != m - 1) {
00371 e[i__ + 1] = r__;
00372 }
00373 g = d__[i__ + 1] - p;
00374 r__ = (d__[i__] - g) * s + c__ * 2. * b;
00375 p = s * r__;
00376 d__[i__ + 1] = g + p;
00377 g = c__ * r__ - b;
00378
00379
00380
00381 if (icompz > 0) {
00382 work[i__] = c__;
00383 work[*n - 1 + i__] = -s;
00384 }
00385
00386
00387 }
00388
00389
00390
00391 if (icompz > 0) {
00392 mm = m - l + 1;
00393 template_lapack_lasr("R", "V", "B", n, &mm, &work[l], &work[*n - 1 + l], &
00394 z___ref(1, l), ldz);
00395 }
00396
00397 d__[l] -= p;
00398 e[l] = g;
00399 goto L40;
00400
00401
00402
00403 L80:
00404 d__[l] = p;
00405
00406 ++l;
00407 if (l <= lend) {
00408 goto L40;
00409 }
00410 goto L140;
00411
00412 } else {
00413
00414
00415
00416
00417
00418 L90:
00419 if (l != lend) {
00420 lendp1 = lend + 1;
00421 i__1 = lendp1;
00422 for (m = l; m >= i__1; --m) {
00423
00424 d__2 = (d__1 = e[m - 1], absMACRO(d__1));
00425 tst = d__2 * d__2;
00426 if (tst <= eps2 * (d__1 = d__[m], absMACRO(d__1)) * (d__2 = d__[m
00427 - 1], absMACRO(d__2)) + safmin) {
00428 goto L110;
00429 }
00430
00431 }
00432 }
00433
00434 m = lend;
00435
00436 L110:
00437 if (m > lend) {
00438 e[m - 1] = 0.;
00439 }
00440 p = d__[l];
00441 if (m == l) {
00442 goto L130;
00443 }
00444
00445
00446
00447
00448 if (m == l - 1) {
00449 if (icompz > 0) {
00450 template_lapack_laev2(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
00451 ;
00452 work[m] = c__;
00453 work[*n - 1 + m] = s;
00454 template_lapack_lasr("R", "V", "F", n, &c__2, &work[m], &work[*n - 1 + m], &
00455 z___ref(1, l - 1), ldz);
00456 } else {
00457 template_lapack_lae2(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
00458 }
00459 d__[l - 1] = rt1;
00460 d__[l] = rt2;
00461 e[l - 1] = 0.;
00462 l += -2;
00463 if (l >= lend) {
00464 goto L90;
00465 }
00466 goto L140;
00467 }
00468
00469 if (jtot == nmaxit) {
00470 goto L140;
00471 }
00472 ++jtot;
00473
00474
00475
00476 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
00477 r__ = template_lapack_lapy2(&g, &c_b10);
00478 g = d__[m] - p + e[l - 1] / (g + template_lapack_d_sign(&r__, &g));
00479
00480 s = 1.;
00481 c__ = 1.;
00482 p = 0.;
00483
00484
00485
00486 lm1 = l - 1;
00487 i__1 = lm1;
00488 for (i__ = m; i__ <= i__1; ++i__) {
00489 f = s * e[i__];
00490 b = c__ * e[i__];
00491 template_lapack_lartg(&g, &f, &c__, &s, &r__);
00492 if (i__ != m) {
00493 e[i__ - 1] = r__;
00494 }
00495 g = d__[i__] - p;
00496 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
00497 p = s * r__;
00498 d__[i__] = g + p;
00499 g = c__ * r__ - b;
00500
00501
00502
00503 if (icompz > 0) {
00504 work[i__] = c__;
00505 work[*n - 1 + i__] = s;
00506 }
00507
00508
00509 }
00510
00511
00512
00513 if (icompz > 0) {
00514 mm = l - m + 1;
00515 template_lapack_lasr("R", "V", "F", n, &mm, &work[m], &work[*n - 1 + m], &
00516 z___ref(1, m), ldz);
00517 }
00518
00519 d__[l] -= p;
00520 e[lm1] = g;
00521 goto L90;
00522
00523
00524
00525 L130:
00526 d__[l] = p;
00527
00528 --l;
00529 if (l >= lend) {
00530 goto L90;
00531 }
00532 goto L140;
00533
00534 }
00535
00536
00537
00538 L140:
00539 if (iscale == 1) {
00540 i__1 = lendsv - lsv + 1;
00541 template_lapack_lascl("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
00542 n, info);
00543 i__1 = lendsv - lsv;
00544 template_lapack_lascl("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
00545 info);
00546 } else if (iscale == 2) {
00547 i__1 = lendsv - lsv + 1;
00548 template_lapack_lascl("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
00549 n, info);
00550 i__1 = lendsv - lsv;
00551 template_lapack_lascl("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
00552 info);
00553 }
00554
00555
00556
00557
00558 if (jtot < nmaxit) {
00559 goto L10;
00560 }
00561 i__1 = *n - 1;
00562 for (i__ = 1; i__ <= i__1; ++i__) {
00563 if (e[i__] != 0.) {
00564 ++(*info);
00565 }
00566
00567 }
00568 goto L190;
00569
00570
00571
00572 L160:
00573 if (icompz == 0) {
00574
00575
00576
00577 template_lapack_lasrt("I", n, &d__[1], info);
00578
00579 } else {
00580
00581
00582
00583 i__1 = *n;
00584 for (ii = 2; ii <= i__1; ++ii) {
00585 i__ = ii - 1;
00586 k = i__;
00587 p = d__[i__];
00588 i__2 = *n;
00589 for (j = ii; j <= i__2; ++j) {
00590 if (d__[j] < p) {
00591 k = j;
00592 p = d__[j];
00593 }
00594
00595 }
00596 if (k != i__) {
00597 d__[k] = d__[i__];
00598 d__[i__] = p;
00599 template_blas_swap(n, &z___ref(1, i__), &c__1, &z___ref(1, k), &c__1);
00600 }
00601
00602 }
00603 }
00604
00605 L190:
00606 return 0;
00607
00608
00609
00610 }
00611
00612 #undef z___ref
00613
00614
00615 #endif