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_STERF_HEADER
00038 #define TEMPLATE_LAPACK_STERF_HEADER
00039
00040 #include "template_lapack_common.h"
00041
00042 template<class Treal>
00043 int template_lapack_sterf(const integer *n, Treal *d__, Treal *e,
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 integer c__0 = 0;
00088 integer c__1 = 1;
00089 Treal c_b32 = 1.;
00090
00091
00092 integer i__1;
00093 Treal d__1, d__2, d__3;
00094
00095 Treal oldc;
00096 integer lend, jtot;
00097 Treal c__;
00098 integer i__, l, m;
00099 Treal p, gamma, r__, s, alpha, sigma, anorm;
00100 integer l1;
00101 Treal bb;
00102 integer iscale;
00103 Treal oldgam, safmin;
00104 Treal safmax;
00105 integer lendsv;
00106 Treal ssfmin;
00107 integer nmaxit;
00108 Treal ssfmax, rt1, rt2, eps, rte;
00109 integer lsv;
00110 Treal eps2;
00111
00112
00113 --e;
00114 --d__;
00115
00116
00117 *info = 0;
00118
00119
00120
00121 if (*n < 0) {
00122 *info = -1;
00123 i__1 = -(*info);
00124 template_blas_erbla("STERF ", &i__1);
00125 return 0;
00126 }
00127 if (*n <= 1) {
00128 return 0;
00129 }
00130
00131
00132
00133 eps = template_lapack_lamch("E", (Treal)0);
00134
00135 d__1 = eps;
00136 eps2 = d__1 * d__1;
00137 safmin = template_lapack_lamch("S", (Treal)0);
00138 safmax = 1. / safmin;
00139 ssfmax = template_blas_sqrt(safmax) / 3.;
00140 ssfmin = template_blas_sqrt(safmin) / eps2;
00141
00142
00143
00144 nmaxit = *n * 30;
00145 sigma = 0.;
00146 jtot = 0;
00147
00148
00149
00150
00151
00152 l1 = 1;
00153
00154 L10:
00155 if (l1 > *n) {
00156 goto L170;
00157 }
00158 if (l1 > 1) {
00159 e[l1 - 1] = 0.;
00160 }
00161 i__1 = *n - 1;
00162 for (m = l1; m <= i__1; ++m) {
00163 if ((d__3 = e[m], absMACRO(d__3)) <= template_blas_sqrt((d__1 = d__[m], absMACRO(d__1))) *
00164 template_blas_sqrt((d__2 = d__[m + 1], absMACRO(d__2))) * eps) {
00165 e[m] = 0.;
00166 goto L30;
00167 }
00168
00169 }
00170 m = *n;
00171
00172 L30:
00173 l = l1;
00174 lsv = l;
00175 lend = m;
00176 lendsv = lend;
00177 l1 = m + 1;
00178 if (lend == l) {
00179 goto L10;
00180 }
00181
00182
00183
00184 i__1 = lend - l + 1;
00185 anorm = template_lapack_lanst("I", &i__1, &d__[l], &e[l]);
00186 iscale = 0;
00187 if (anorm > ssfmax) {
00188 iscale = 1;
00189 i__1 = lend - l + 1;
00190 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
00191 info);
00192 i__1 = lend - l;
00193 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
00194 info);
00195 } else if (anorm < ssfmin) {
00196 iscale = 2;
00197 i__1 = lend - l + 1;
00198 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
00199 info);
00200 i__1 = lend - l;
00201 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
00202 info);
00203 }
00204
00205 i__1 = lend - 1;
00206 for (i__ = l; i__ <= i__1; ++i__) {
00207
00208 d__1 = e[i__];
00209 e[i__] = d__1 * d__1;
00210
00211 }
00212
00213
00214
00215 if ((d__1 = d__[lend], absMACRO(d__1)) < (d__2 = d__[l], absMACRO(d__2))) {
00216 lend = lsv;
00217 l = lendsv;
00218 }
00219
00220 if (lend >= l) {
00221
00222
00223
00224
00225
00226 L50:
00227 if (l != lend) {
00228 i__1 = lend - 1;
00229 for (m = l; m <= i__1; ++m) {
00230 if ((d__2 = e[m], absMACRO(d__2)) <= eps2 * (d__1 = d__[m] * d__[m
00231 + 1], absMACRO(d__1))) {
00232 goto L70;
00233 }
00234
00235 }
00236 }
00237 m = lend;
00238
00239 L70:
00240 if (m < lend) {
00241 e[m] = 0.;
00242 }
00243 p = d__[l];
00244 if (m == l) {
00245 goto L90;
00246 }
00247
00248
00249
00250
00251 if (m == l + 1) {
00252 rte = template_blas_sqrt(e[l]);
00253 template_lapack_lae2(&d__[l], &rte, &d__[l + 1], &rt1, &rt2);
00254 d__[l] = rt1;
00255 d__[l + 1] = rt2;
00256 e[l] = 0.;
00257 l += 2;
00258 if (l <= lend) {
00259 goto L50;
00260 }
00261 goto L150;
00262 }
00263
00264 if (jtot == nmaxit) {
00265 goto L150;
00266 }
00267 ++jtot;
00268
00269
00270
00271 rte = template_blas_sqrt(e[l]);
00272 sigma = (d__[l + 1] - p) / (rte * 2.);
00273 r__ = template_lapack_lapy2(&sigma, &c_b32);
00274 sigma = p - rte / (sigma + template_lapack_d_sign(&r__, &sigma));
00275
00276 c__ = 1.;
00277 s = 0.;
00278 gamma = d__[m] - sigma;
00279 p = gamma * gamma;
00280
00281
00282
00283 i__1 = l;
00284 for (i__ = m - 1; i__ >= i__1; --i__) {
00285 bb = e[i__];
00286 r__ = p + bb;
00287 if (i__ != m - 1) {
00288 e[i__ + 1] = s * r__;
00289 }
00290 oldc = c__;
00291 c__ = p / r__;
00292 s = bb / r__;
00293 oldgam = gamma;
00294 alpha = d__[i__];
00295 gamma = c__ * (alpha - sigma) - s * oldgam;
00296 d__[i__ + 1] = oldgam + (alpha - gamma);
00297 if (c__ != 0.) {
00298 p = gamma * gamma / c__;
00299 } else {
00300 p = oldc * bb;
00301 }
00302
00303 }
00304
00305 e[l] = s * p;
00306 d__[l] = sigma + gamma;
00307 goto L50;
00308
00309
00310
00311 L90:
00312 d__[l] = p;
00313
00314 ++l;
00315 if (l <= lend) {
00316 goto L50;
00317 }
00318 goto L150;
00319
00320 } else {
00321
00322
00323
00324
00325
00326 L100:
00327 i__1 = lend + 1;
00328 for (m = l; m >= i__1; --m) {
00329 if ((d__2 = e[m - 1], absMACRO(d__2)) <= eps2 * (d__1 = d__[m] * d__[m
00330 - 1], absMACRO(d__1))) {
00331 goto L120;
00332 }
00333
00334 }
00335 m = lend;
00336
00337 L120:
00338 if (m > lend) {
00339 e[m - 1] = 0.;
00340 }
00341 p = d__[l];
00342 if (m == l) {
00343 goto L140;
00344 }
00345
00346
00347
00348
00349 if (m == l - 1) {
00350 rte = template_blas_sqrt(e[l - 1]);
00351 template_lapack_lae2(&d__[l], &rte, &d__[l - 1], &rt1, &rt2);
00352 d__[l] = rt1;
00353 d__[l - 1] = rt2;
00354 e[l - 1] = 0.;
00355 l += -2;
00356 if (l >= lend) {
00357 goto L100;
00358 }
00359 goto L150;
00360 }
00361
00362 if (jtot == nmaxit) {
00363 goto L150;
00364 }
00365 ++jtot;
00366
00367
00368
00369 rte = template_blas_sqrt(e[l - 1]);
00370 sigma = (d__[l - 1] - p) / (rte * 2.);
00371 r__ = template_lapack_lapy2(&sigma, &c_b32);
00372 sigma = p - rte / (sigma + template_lapack_d_sign(&r__, &sigma));
00373
00374 c__ = 1.;
00375 s = 0.;
00376 gamma = d__[m] - sigma;
00377 p = gamma * gamma;
00378
00379
00380
00381 i__1 = l - 1;
00382 for (i__ = m; i__ <= i__1; ++i__) {
00383 bb = e[i__];
00384 r__ = p + bb;
00385 if (i__ != m) {
00386 e[i__ - 1] = s * r__;
00387 }
00388 oldc = c__;
00389 c__ = p / r__;
00390 s = bb / r__;
00391 oldgam = gamma;
00392 alpha = d__[i__ + 1];
00393 gamma = c__ * (alpha - sigma) - s * oldgam;
00394 d__[i__] = oldgam + (alpha - gamma);
00395 if (c__ != 0.) {
00396 p = gamma * gamma / c__;
00397 } else {
00398 p = oldc * bb;
00399 }
00400
00401 }
00402
00403 e[l - 1] = s * p;
00404 d__[l] = sigma + gamma;
00405 goto L100;
00406
00407
00408
00409 L140:
00410 d__[l] = p;
00411
00412 --l;
00413 if (l >= lend) {
00414 goto L100;
00415 }
00416 goto L150;
00417
00418 }
00419
00420
00421
00422 L150:
00423 if (iscale == 1) {
00424 i__1 = lendsv - lsv + 1;
00425 template_lapack_lascl("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
00426 n, info);
00427 }
00428 if (iscale == 2) {
00429 i__1 = lendsv - lsv + 1;
00430 template_lapack_lascl("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
00431 n, info);
00432 }
00433
00434
00435
00436
00437 if (jtot < nmaxit) {
00438 goto L10;
00439 }
00440 i__1 = *n - 1;
00441 for (i__ = 1; i__ <= i__1; ++i__) {
00442 if (e[i__] != 0.) {
00443 ++(*info);
00444 }
00445
00446 }
00447 goto L180;
00448
00449
00450
00451 L170:
00452 template_lapack_lasrt("I", n, &d__[1], info);
00453
00454 L180:
00455 return 0;
00456
00457
00458
00459 }
00460
00461 #endif