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