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