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