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_LASQ2_HEADER
00038 #define TEMPLATE_LAPACK_LASQ2_HEADER
00039
00040
00041 #include "template_lapack_lasq3.h"
00042
00043
00044 template<class Treal>
00045 int template_lapack_lasq2(integer *n, Treal *z__, integer *info)
00046 {
00047
00048 integer i__1, i__2, i__3;
00049 Treal d__1, d__2;
00050
00051
00052
00053 Treal d__, e, g;
00054 integer k;
00055 Treal s, t;
00056 integer i0, i4, n0;
00057 Treal dn;
00058 integer pp;
00059 Treal dn1, dn2, dee, eps, tau, tol;
00060 integer ipn4;
00061 Treal tol2;
00062 logical ieee;
00063 integer nbig;
00064 Treal dmin__, emin, emax;
00065 integer kmin, ndiv, iter;
00066 Treal qmin, temp, qmax, zmax;
00067 integer splt;
00068 Treal dmin1, dmin2;
00069 integer nfail;
00070 Treal desig, trace, sigma;
00071 integer iinfo, ttype;
00072 Treal deemin;
00073 integer iwhila, iwhilb;
00074 Treal oldemn, safmin;
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 integer c__1 = 1;
00164 integer c__2 = 2;
00165 integer c__10 = 10;
00166 integer c__3 = 3;
00167 integer c__4 = 4;
00168 integer c__11 = 11;
00169
00170
00171 --z__;
00172
00173
00174 *info = 0;
00175 eps = template_lapack_lamch("Precision", (Treal)0);
00176 safmin = template_lapack_lamch("Safe minimum", (Treal)0);
00177 tol = eps * 100.;
00178
00179 d__1 = tol;
00180 tol2 = d__1 * d__1;
00181
00182 if (*n < 0) {
00183 *info = -1;
00184 template_blas_erbla("DLASQ2", &c__1);
00185 return 0;
00186 } else if (*n == 0) {
00187 return 0;
00188 } else if (*n == 1) {
00189
00190
00191
00192 if (z__[1] < 0.) {
00193 *info = -201;
00194 template_blas_erbla("DLASQ2", &c__2);
00195 }
00196 return 0;
00197 } else if (*n == 2) {
00198
00199
00200
00201 if (z__[2] < 0. || z__[3] < 0.) {
00202 *info = -2;
00203 template_blas_erbla("DLASQ2", &c__2);
00204 return 0;
00205 } else if (z__[3] > z__[1]) {
00206 d__ = z__[3];
00207 z__[3] = z__[1];
00208 z__[1] = d__;
00209 }
00210 z__[5] = z__[1] + z__[2] + z__[3];
00211 if (z__[2] > z__[3] * tol2) {
00212 t = (z__[1] - z__[3] + z__[2]) * .5;
00213 s = z__[3] * (z__[2] / t);
00214 if (s <= t) {
00215 s = z__[3] * (z__[2] / (t * (template_blas_sqrt(s / t + 1.) + 1.)));
00216 } else {
00217 s = z__[3] * (z__[2] / (t + template_blas_sqrt(t) * template_blas_sqrt(t + s)));
00218 }
00219 t = z__[1] + (s + z__[2]);
00220 z__[3] *= z__[1] / t;
00221 z__[1] = t;
00222 }
00223 z__[2] = z__[3];
00224 z__[6] = z__[2] + z__[1];
00225 return 0;
00226 }
00227
00228
00229
00230 z__[*n * 2] = 0.;
00231 emin = z__[2];
00232 qmax = 0.;
00233 zmax = 0.;
00234 d__ = 0.;
00235 e = 0.;
00236
00237 i__1 = ( *n - 1 ) << 1;
00238 for (k = 1; k <= i__1; k += 2) {
00239 if (z__[k] < 0.) {
00240 *info = -(k + 200);
00241 template_blas_erbla("DLASQ2", &c__2);
00242 return 0;
00243 } else if (z__[k + 1] < 0.) {
00244 *info = -(k + 201);
00245 template_blas_erbla("DLASQ2", &c__2);
00246 return 0;
00247 }
00248 d__ += z__[k];
00249 e += z__[k + 1];
00250
00251 d__1 = qmax, d__2 = z__[k];
00252 qmax = maxMACRO(d__1,d__2);
00253
00254 d__1 = emin, d__2 = z__[k + 1];
00255 emin = minMACRO(d__1,d__2);
00256
00257 d__1 = maxMACRO(qmax,zmax), d__2 = z__[k + 1];
00258 zmax = maxMACRO(d__1,d__2);
00259
00260 }
00261 if (z__[(*n << 1) - 1] < 0.) {
00262 *info = -((*n << 1) + 199);
00263 template_blas_erbla("DLASQ2", &c__2);
00264 return 0;
00265 }
00266 d__ += z__[(*n << 1) - 1];
00267
00268 d__1 = qmax, d__2 = z__[(*n << 1) - 1];
00269 qmax = maxMACRO(d__1,d__2);
00270 zmax = maxMACRO(qmax,zmax);
00271
00272
00273
00274 if (e == 0.) {
00275 i__1 = *n;
00276 for (k = 2; k <= i__1; ++k) {
00277 z__[k] = z__[(k << 1) - 1];
00278
00279 }
00280 template_lapack_lasrt("D", n, &z__[1], &iinfo);
00281 z__[(*n << 1) - 1] = d__;
00282 return 0;
00283 }
00284
00285 trace = d__ + e;
00286
00287
00288
00289 if (trace == 0.) {
00290 z__[(*n << 1) - 1] = 0.;
00291 return 0;
00292 }
00293
00294
00295
00296 ieee = template_lapack_ilaenv(&c__10, "DLASQ2", "N", &c__1, &c__2, &c__3, &c__4, (ftnlen)6, (ftnlen)1) == 1 && template_lapack_ilaenv(&c__11, "DLASQ2", "N", &c__1, &c__2,
00297 &c__3, &c__4, (ftnlen)6, (ftnlen)1) == 1;
00298
00299
00300
00301 for (k = *n << 1; k >= 2; k += -2) {
00302 z__[k * 2] = 0.;
00303 z__[(k << 1) - 1] = z__[k];
00304 z__[(k << 1) - 2] = 0.;
00305 z__[(k << 1) - 3] = z__[k - 1];
00306
00307 }
00308
00309 i0 = 1;
00310 n0 = *n;
00311
00312
00313
00314 if (z__[(i0 << 2) - 3] * 1.5 < z__[(n0 << 2) - 3]) {
00315 ipn4 = ( i0 + n0 ) << 2;
00316 i__1 = ( i0 + n0 - 1 ) << 1;
00317 for (i4 = i0 << 2; i4 <= i__1; i4 += 4) {
00318 temp = z__[i4 - 3];
00319 z__[i4 - 3] = z__[ipn4 - i4 - 3];
00320 z__[ipn4 - i4 - 3] = temp;
00321 temp = z__[i4 - 1];
00322 z__[i4 - 1] = z__[ipn4 - i4 - 5];
00323 z__[ipn4 - i4 - 5] = temp;
00324
00325 }
00326 }
00327
00328
00329
00330 pp = 0;
00331
00332 for (k = 1; k <= 2; ++k) {
00333
00334 d__ = z__[(n0 << 2) + pp - 3];
00335 i__1 = (i0 << 2) + pp;
00336 for (i4 = ( ( n0 - 1 ) << 2) + pp; i4 >= i__1; i4 += -4) {
00337 if (z__[i4 - 1] <= tol2 * d__) {
00338 z__[i4 - 1] = -0.;
00339 d__ = z__[i4 - 3];
00340 } else {
00341 d__ = z__[i4 - 3] * (d__ / (d__ + z__[i4 - 1]));
00342 }
00343
00344 }
00345
00346
00347
00348 emin = z__[(i0 << 2) + pp + 1];
00349 d__ = z__[(i0 << 2) + pp - 3];
00350 i__1 = ( ( n0 - 1 ) << 2) + pp;
00351 for (i4 = (i0 << 2) + pp; i4 <= i__1; i4 += 4) {
00352 z__[i4 - (pp << 1) - 2] = d__ + z__[i4 - 1];
00353 if (z__[i4 - 1] <= tol2 * d__) {
00354 z__[i4 - 1] = -0.;
00355 z__[i4 - (pp << 1) - 2] = d__;
00356 z__[i4 - (pp << 1)] = 0.;
00357 d__ = z__[i4 + 1];
00358 } else if (safmin * z__[i4 + 1] < z__[i4 - (pp << 1) - 2] &&
00359 safmin * z__[i4 - (pp << 1) - 2] < z__[i4 + 1]) {
00360 temp = z__[i4 + 1] / z__[i4 - (pp << 1) - 2];
00361 z__[i4 - (pp << 1)] = z__[i4 - 1] * temp;
00362 d__ *= temp;
00363 } else {
00364 z__[i4 - (pp << 1)] = z__[i4 + 1] * (z__[i4 - 1] / z__[i4 - (
00365 pp << 1) - 2]);
00366 d__ = z__[i4 + 1] * (d__ / z__[i4 - (pp << 1) - 2]);
00367 }
00368
00369 d__1 = emin, d__2 = z__[i4 - (pp << 1)];
00370 emin = minMACRO(d__1,d__2);
00371
00372 }
00373 z__[(n0 << 2) - pp - 2] = d__;
00374
00375
00376
00377 qmax = z__[(i0 << 2) - pp - 2];
00378 i__1 = (n0 << 2) - pp - 2;
00379 for (i4 = (i0 << 2) - pp + 2; i4 <= i__1; i4 += 4) {
00380
00381 d__1 = qmax, d__2 = z__[i4];
00382 qmax = maxMACRO(d__1,d__2);
00383
00384 }
00385
00386
00387
00388 pp = 1 - pp;
00389
00390 }
00391
00392
00393
00394 ttype = 0;
00395 dmin1 = 0.;
00396 dmin2 = 0.;
00397 dn = 0.;
00398 dn1 = 0.;
00399 dn2 = 0.;
00400 g = 0.;
00401 tau = 0.;
00402
00403 iter = 2;
00404 nfail = 0;
00405 ndiv = ( n0 - i0 ) << 1;
00406
00407 i__1 = *n + 1;
00408 for (iwhila = 1; iwhila <= i__1; ++iwhila) {
00409 if (n0 < 1) {
00410 goto L170;
00411 }
00412
00413
00414
00415
00416
00417
00418 desig = 0.;
00419 if (n0 == *n) {
00420 sigma = 0.;
00421 } else {
00422 sigma = -z__[(n0 << 2) - 1];
00423 }
00424 if (sigma < 0.) {
00425 *info = 1;
00426 return 0;
00427 }
00428
00429
00430
00431
00432 emax = 0.;
00433 if (n0 > i0) {
00434 emin = (d__1 = z__[(n0 << 2) - 5], absMACRO(d__1));
00435 } else {
00436 emin = 0.;
00437 }
00438 qmin = z__[(n0 << 2) - 3];
00439 qmax = qmin;
00440 for (i4 = n0 << 2; i4 >= 8; i4 += -4) {
00441 if (z__[i4 - 5] <= 0.) {
00442 goto L100;
00443 }
00444 if (qmin >= emax * 4.) {
00445
00446 d__1 = qmin, d__2 = z__[i4 - 3];
00447 qmin = minMACRO(d__1,d__2);
00448
00449 d__1 = emax, d__2 = z__[i4 - 5];
00450 emax = maxMACRO(d__1,d__2);
00451 }
00452
00453 d__1 = qmax, d__2 = z__[i4 - 7] + z__[i4 - 5];
00454 qmax = maxMACRO(d__1,d__2);
00455
00456 d__1 = emin, d__2 = z__[i4 - 5];
00457 emin = minMACRO(d__1,d__2);
00458
00459 }
00460 i4 = 4;
00461
00462 L100:
00463 i0 = i4 / 4;
00464 pp = 0;
00465
00466 if (n0 - i0 > 1) {
00467 dee = z__[(i0 << 2) - 3];
00468 deemin = dee;
00469 kmin = i0;
00470 i__2 = (n0 << 2) - 3;
00471 for (i4 = (i0 << 2) + 1; i4 <= i__2; i4 += 4) {
00472 dee = z__[i4] * (dee / (dee + z__[i4 - 2]));
00473 if (dee <= deemin) {
00474 deemin = dee;
00475 kmin = (i4 + 3) / 4;
00476 }
00477
00478 }
00479 if ( ( kmin - i0 ) << 1 < n0 - kmin && deemin <= z__[(n0 << 2) - 3] *
00480 .5) {
00481 ipn4 = ( i0 + n0 ) << 2;
00482 pp = 2;
00483 i__2 = ( i0 + n0 - 1 ) << 1;
00484 for (i4 = i0 << 2; i4 <= i__2; i4 += 4) {
00485 temp = z__[i4 - 3];
00486 z__[i4 - 3] = z__[ipn4 - i4 - 3];
00487 z__[ipn4 - i4 - 3] = temp;
00488 temp = z__[i4 - 2];
00489 z__[i4 - 2] = z__[ipn4 - i4 - 2];
00490 z__[ipn4 - i4 - 2] = temp;
00491 temp = z__[i4 - 1];
00492 z__[i4 - 1] = z__[ipn4 - i4 - 5];
00493 z__[ipn4 - i4 - 5] = temp;
00494 temp = z__[i4];
00495 z__[i4] = z__[ipn4 - i4 - 4];
00496 z__[ipn4 - i4 - 4] = temp;
00497
00498 }
00499 }
00500 }
00501
00502
00503
00504
00505 d__1 = 0., d__2 = qmin - template_blas_sqrt(qmin) * 2. * template_blas_sqrt(emax);
00506 dmin__ = -maxMACRO(d__1,d__2);
00507
00508
00509
00510
00511
00512
00513
00514 nbig = (n0 - i0 + 1) * 30;
00515 i__2 = nbig;
00516 for (iwhilb = 1; iwhilb <= i__2; ++iwhilb) {
00517 if (i0 > n0) {
00518 goto L150;
00519 }
00520
00521
00522
00523 template_lapack_lasq3(&i0, &n0, &z__[1], &pp, &dmin__, &sigma, &desig, &qmax, &
00524 nfail, &iter, &ndiv, &ieee, &ttype, &dmin1, &dmin2, &dn, &
00525 dn1, &dn2, &g, &tau);
00526
00527 pp = 1 - pp;
00528
00529
00530
00531 if (pp == 0 && n0 - i0 >= 3) {
00532 if (z__[n0 * 4] <= tol2 * qmax || z__[(n0 << 2) - 1] <= tol2 *
00533 sigma) {
00534 splt = i0 - 1;
00535 qmax = z__[(i0 << 2) - 3];
00536 emin = z__[(i0 << 2) - 1];
00537 oldemn = z__[i0 * 4];
00538 i__3 = ( n0 - 3 ) << 2;
00539 for (i4 = i0 << 2; i4 <= i__3; i4 += 4) {
00540 if (z__[i4] <= tol2 * z__[i4 - 3] || z__[i4 - 1] <=
00541 tol2 * sigma) {
00542 z__[i4 - 1] = -sigma;
00543 splt = i4 / 4;
00544 qmax = 0.;
00545 emin = z__[i4 + 3];
00546 oldemn = z__[i4 + 4];
00547 } else {
00548
00549 d__1 = qmax, d__2 = z__[i4 + 1];
00550 qmax = maxMACRO(d__1,d__2);
00551
00552 d__1 = emin, d__2 = z__[i4 - 1];
00553 emin = minMACRO(d__1,d__2);
00554
00555 d__1 = oldemn, d__2 = z__[i4];
00556 oldemn = minMACRO(d__1,d__2);
00557 }
00558
00559 }
00560 z__[(n0 << 2) - 1] = emin;
00561 z__[n0 * 4] = oldemn;
00562 i0 = splt + 1;
00563 }
00564 }
00565
00566
00567 }
00568
00569 *info = 2;
00570 return 0;
00571
00572
00573
00574 L150:
00575
00576
00577 ;
00578 }
00579
00580 *info = 3;
00581 return 0;
00582
00583
00584
00585 L170:
00586
00587
00588
00589 i__1 = *n;
00590 for (k = 2; k <= i__1; ++k) {
00591 z__[k] = z__[(k << 2) - 3];
00592
00593 }
00594
00595
00596
00597 template_lapack_lasrt("D", n, &z__[1], &iinfo);
00598
00599 e = 0.;
00600 for (k = *n; k >= 1; --k) {
00601 e += z__[k];
00602
00603 }
00604
00605
00606
00607 z__[(*n << 1) + 1] = trace;
00608 z__[(*n << 1) + 2] = e;
00609 z__[(*n << 1) + 3] = (Treal) iter;
00610
00611 i__1 = *n;
00612 z__[(*n << 1) + 4] = (Treal) ndiv / (Treal) (i__1 * i__1);
00613 z__[(*n << 1) + 5] = nfail * 100. / (Treal) iter;
00614 return 0;
00615
00616
00617
00618 }
00619
00620 #endif