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_HGEQZ_HEADER
00036 #define TEMPLATE_LAPACK_HGEQZ_HEADER
00037
00038
00039 template<class Treal>
00040 int template_lapack_hgeqz(const char *job, const char *compq, const char *compz, const integer *n,
00041 const integer *ilo, const integer *ihi, Treal *a, const integer *lda, Treal *
00042 b, const integer *ldb, Treal *alphar, Treal *alphai, Treal *
00043 beta, Treal *q, const integer *ldq, Treal *z__, const integer *ldz,
00044 Treal *work, const integer *lwork, 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
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
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248 Treal c_b12 = 0.;
00249 Treal c_b13 = 1.;
00250 integer c__1 = 1;
00251 integer c__3 = 3;
00252
00253
00254 integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1,
00255 z_offset, i__1, i__2, i__3, i__4;
00256 Treal d__1, d__2, d__3, d__4;
00257
00258 Treal ad11l, ad12l, ad21l, ad22l, ad32l, wabs, atol, btol,
00259 temp;
00260 Treal temp2, s1inv, c__;
00261 integer j;
00262 Treal s, t, v[3], scale;
00263 integer iiter, ilast, jiter;
00264 Treal anorm, bnorm;
00265 integer maxit;
00266 Treal tempi, tempr, s1, s2, u1, u2;
00267 logical ilazr2;
00268 Treal a11, a12, a21, a22, b11, b22, c12, c21;
00269 integer jc;
00270 Treal an, bn, cl, cq, cr;
00271 integer in;
00272 Treal ascale, bscale, u12, w11;
00273 integer jr;
00274 Treal cz, sl, w12, w21, w22, wi;
00275 Treal sr;
00276 Treal vs, wr;
00277 Treal safmin;
00278 Treal safmax;
00279 Treal eshift;
00280 logical ilschr;
00281 Treal b1a, b2a;
00282 integer icompq, ilastm;
00283 Treal a1i;
00284 integer ischur;
00285 Treal a2i, b1i;
00286 logical ilazro;
00287 integer icompz, ifirst;
00288 Treal b2i;
00289 integer ifrstm;
00290 Treal a1r;
00291 integer istart;
00292 logical ilpivt;
00293 Treal a2r, b1r, b2r;
00294 logical lquery;
00295 Treal wr2, ad11, ad12, ad21, ad22, c11i, c22i;
00296 integer jch;
00297 Treal c11r, c22r, u12l;
00298 logical ilq;
00299 Treal tau, sqi;
00300 logical ilz;
00301 Treal ulp, sqr, szi, szr;
00302 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00303 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
00304 #define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1]
00305 #define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
00306
00307
00308 a_dim1 = *lda;
00309 a_offset = 1 + a_dim1 * 1;
00310 a -= a_offset;
00311 b_dim1 = *ldb;
00312 b_offset = 1 + b_dim1 * 1;
00313 b -= b_offset;
00314 --alphar;
00315 --alphai;
00316 --beta;
00317 q_dim1 = *ldq;
00318 q_offset = 1 + q_dim1 * 1;
00319 q -= q_offset;
00320 z_dim1 = *ldz;
00321 z_offset = 1 + z_dim1 * 1;
00322 z__ -= z_offset;
00323 --work;
00324
00325
00326 ilschr = ilq = ilz = 0;
00327
00328 if (template_blas_lsame(job, "E")) {
00329 ilschr = FALSE_;
00330 ischur = 1;
00331 } else if (template_blas_lsame(job, "S")) {
00332 ilschr = TRUE_;
00333 ischur = 2;
00334 } else {
00335 ischur = 0;
00336 }
00337
00338 if (template_blas_lsame(compq, "N")) {
00339 ilq = FALSE_;
00340 icompq = 1;
00341 } else if (template_blas_lsame(compq, "V")) {
00342 ilq = TRUE_;
00343 icompq = 2;
00344 } else if (template_blas_lsame(compq, "I")) {
00345 ilq = TRUE_;
00346 icompq = 3;
00347 } else {
00348 icompq = 0;
00349 }
00350
00351 if (template_blas_lsame(compz, "N")) {
00352 ilz = FALSE_;
00353 icompz = 1;
00354 } else if (template_blas_lsame(compz, "V")) {
00355 ilz = TRUE_;
00356 icompz = 2;
00357 } else if (template_blas_lsame(compz, "I")) {
00358 ilz = TRUE_;
00359 icompz = 3;
00360 } else {
00361 icompz = 0;
00362 }
00363
00364
00365
00366 *info = 0;
00367 work[1] = (Treal) maxMACRO(1,*n);
00368 lquery = *lwork == -1;
00369 if (ischur == 0) {
00370 *info = -1;
00371 } else if (icompq == 0) {
00372 *info = -2;
00373 } else if (icompz == 0) {
00374 *info = -3;
00375 } else if (*n < 0) {
00376 *info = -4;
00377 } else if (*ilo < 1) {
00378 *info = -5;
00379 } else if (*ihi > *n || *ihi < *ilo - 1) {
00380 *info = -6;
00381 } else if (*lda < *n) {
00382 *info = -8;
00383 } else if (*ldb < *n) {
00384 *info = -10;
00385 } else if (*ldq < 1 || ( ilq && *ldq < *n ) ) {
00386 *info = -15;
00387 } else if (*ldz < 1 || ( ilz && *ldz < *n ) ) {
00388 *info = -17;
00389 } else if (*lwork < maxMACRO(1,*n) && ! lquery) {
00390 *info = -19;
00391 }
00392 if (*info != 0) {
00393 i__1 = -(*info);
00394 template_blas_erbla("HGEQZ ", &i__1);
00395 return 0;
00396 } else if (lquery) {
00397 return 0;
00398 }
00399
00400
00401
00402 if (*n <= 0) {
00403 work[1] = 1.;
00404 return 0;
00405 }
00406
00407
00408
00409 if (icompq == 3) {
00410 template_lapack_laset("Full", n, n, &c_b12, &c_b13, &q[q_offset], ldq);
00411 }
00412 if (icompz == 3) {
00413 template_lapack_laset("Full", n, n, &c_b12, &c_b13, &z__[z_offset], ldz);
00414 }
00415
00416
00417
00418 in = *ihi + 1 - *ilo;
00419 safmin = template_lapack_lamch("S", (Treal)0);
00420 safmax = 1. / safmin;
00421 ulp = template_lapack_lamch("E", (Treal)0) * template_lapack_lamch("B", (Treal)0);
00422 anorm = dlanhs_("F", &in, &a_ref(*ilo, *ilo), lda, &work[1]);
00423 bnorm = dlanhs_("F", &in, &b_ref(*ilo, *ilo), ldb, &work[1]);
00424
00425 d__1 = safmin, d__2 = ulp * anorm;
00426 atol = maxMACRO(d__1,d__2);
00427
00428 d__1 = safmin, d__2 = ulp * bnorm;
00429 btol = maxMACRO(d__1,d__2);
00430 ascale = 1. / maxMACRO(safmin,anorm);
00431 bscale = 1. / maxMACRO(safmin,bnorm);
00432
00433
00434
00435 i__1 = *n;
00436 for (j = *ihi + 1; j <= i__1; ++j) {
00437 if (b_ref(j, j) < 0.) {
00438 if (ilschr) {
00439 i__2 = j;
00440 for (jr = 1; jr <= i__2; ++jr) {
00441 a_ref(jr, j) = -a_ref(jr, j);
00442 b_ref(jr, j) = -b_ref(jr, j);
00443
00444 }
00445 } else {
00446 a_ref(j, j) = -a_ref(j, j);
00447 b_ref(j, j) = -b_ref(j, j);
00448 }
00449 if (ilz) {
00450 i__2 = *n;
00451 for (jr = 1; jr <= i__2; ++jr) {
00452 z___ref(jr, j) = -z___ref(jr, j);
00453
00454 }
00455 }
00456 }
00457 alphar[j] = a_ref(j, j);
00458 alphai[j] = 0.;
00459 beta[j] = b_ref(j, j);
00460
00461 }
00462
00463
00464
00465 if (*ihi < *ilo) {
00466 goto L380;
00467 }
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484 ilast = *ihi;
00485 if (ilschr) {
00486 ifrstm = 1;
00487 ilastm = *n;
00488 } else {
00489 ifrstm = *ilo;
00490 ilastm = *ihi;
00491 }
00492 iiter = 0;
00493 eshift = 0.;
00494 maxit = (*ihi - *ilo + 1) * 30;
00495
00496 i__1 = maxit;
00497 for (jiter = 1; jiter <= i__1; ++jiter) {
00498
00499
00500
00501
00502
00503
00504
00505 if (ilast == *ilo) {
00506
00507
00508
00509 goto L80;
00510 } else {
00511 if ((d__1 = a_ref(ilast, ilast - 1), absMACRO(d__1)) <= atol) {
00512 a_ref(ilast, ilast - 1) = 0.;
00513 goto L80;
00514 }
00515 }
00516
00517 if ((d__1 = b_ref(ilast, ilast), absMACRO(d__1)) <= btol) {
00518 b_ref(ilast, ilast) = 0.;
00519 goto L70;
00520 }
00521
00522
00523
00524 i__2 = *ilo;
00525 for (j = ilast - 1; j >= i__2; --j) {
00526
00527
00528
00529 if (j == *ilo) {
00530 ilazro = TRUE_;
00531 } else {
00532 if ((d__1 = a_ref(j, j - 1), absMACRO(d__1)) <= atol) {
00533 a_ref(j, j - 1) = 0.;
00534 ilazro = TRUE_;
00535 } else {
00536 ilazro = FALSE_;
00537 }
00538 }
00539
00540
00541
00542 if ((d__1 = b_ref(j, j), absMACRO(d__1)) < btol) {
00543 b_ref(j, j) = 0.;
00544
00545
00546
00547 ilazr2 = FALSE_;
00548 if (! ilazro) {
00549 temp = (d__1 = a_ref(j, j - 1), absMACRO(d__1));
00550 temp2 = (d__1 = a_ref(j, j), absMACRO(d__1));
00551 tempr = maxMACRO(temp,temp2);
00552 if (tempr < 1. && tempr != 0.) {
00553 temp /= tempr;
00554 temp2 /= tempr;
00555 }
00556 if (temp * (ascale * (d__1 = a_ref(j + 1, j), absMACRO(d__1)))
00557 <= temp2 * (ascale * atol)) {
00558 ilazr2 = TRUE_;
00559 }
00560 }
00561
00562
00563
00564
00565
00566
00567
00568 if (ilazro || ilazr2) {
00569 i__3 = ilast - 1;
00570 for (jch = j; jch <= i__3; ++jch) {
00571 temp = a_ref(jch, jch);
00572 template_lapack_lartg(&temp, &a_ref(jch + 1, jch), &c__, &s, &a_ref(
00573 jch, jch));
00574 a_ref(jch + 1, jch) = 0.;
00575 i__4 = ilastm - jch;
00576 template_blas_rot(&i__4, &a_ref(jch, jch + 1), lda, &a_ref(jch +
00577 1, jch + 1), lda, &c__, &s);
00578 i__4 = ilastm - jch;
00579 template_blas_rot(&i__4, &b_ref(jch, jch + 1), ldb, &b_ref(jch +
00580 1, jch + 1), ldb, &c__, &s);
00581 if (ilq) {
00582 template_blas_rot(n, &q_ref(1, jch), &c__1, &q_ref(1, jch + 1)
00583 , &c__1, &c__, &s);
00584 }
00585 if (ilazr2) {
00586 a_ref(jch, jch - 1) = a_ref(jch, jch - 1) * c__;
00587 }
00588 ilazr2 = FALSE_;
00589 if ((d__1 = b_ref(jch + 1, jch + 1), absMACRO(d__1)) >=
00590 btol) {
00591 if (jch + 1 >= ilast) {
00592 goto L80;
00593 } else {
00594 ifirst = jch + 1;
00595 goto L110;
00596 }
00597 }
00598 b_ref(jch + 1, jch + 1) = 0.;
00599
00600 }
00601 goto L70;
00602 } else {
00603
00604
00605
00606
00607 i__3 = ilast - 1;
00608 for (jch = j; jch <= i__3; ++jch) {
00609 temp = b_ref(jch, jch + 1);
00610 template_lapack_lartg(&temp, &b_ref(jch + 1, jch + 1), &c__, &s, &
00611 b_ref(jch, jch + 1));
00612 b_ref(jch + 1, jch + 1) = 0.;
00613 if (jch < ilastm - 1) {
00614 i__4 = ilastm - jch - 1;
00615 template_blas_rot(&i__4, &b_ref(jch, jch + 2), ldb, &b_ref(
00616 jch + 1, jch + 2), ldb, &c__, &s);
00617 }
00618 i__4 = ilastm - jch + 2;
00619 template_blas_rot(&i__4, &a_ref(jch, jch - 1), lda, &a_ref(jch +
00620 1, jch - 1), lda, &c__, &s);
00621 if (ilq) {
00622 template_blas_rot(n, &q_ref(1, jch), &c__1, &q_ref(1, jch + 1)
00623 , &c__1, &c__, &s);
00624 }
00625 temp = a_ref(jch + 1, jch);
00626 template_lapack_lartg(&temp, &a_ref(jch + 1, jch - 1), &c__, &s, &
00627 a_ref(jch + 1, jch));
00628 a_ref(jch + 1, jch - 1) = 0.;
00629 i__4 = jch + 1 - ifrstm;
00630 template_blas_rot(&i__4, &a_ref(ifrstm, jch), &c__1, &a_ref(
00631 ifrstm, jch - 1), &c__1, &c__, &s);
00632 i__4 = jch - ifrstm;
00633 template_blas_rot(&i__4, &b_ref(ifrstm, jch), &c__1, &b_ref(
00634 ifrstm, jch - 1), &c__1, &c__, &s);
00635 if (ilz) {
00636 template_blas_rot(n, &z___ref(1, jch), &c__1, &z___ref(1, jch
00637 - 1), &c__1, &c__, &s);
00638 }
00639
00640 }
00641 goto L70;
00642 }
00643 } else if (ilazro) {
00644
00645
00646
00647 ifirst = j;
00648 goto L110;
00649 }
00650
00651
00652
00653
00654 }
00655
00656
00657
00658 *info = *n + 1;
00659 goto L420;
00660
00661
00662
00663
00664 L70:
00665 temp = a_ref(ilast, ilast);
00666 template_lapack_lartg(&temp, &a_ref(ilast, ilast - 1), &c__, &s, &a_ref(ilast,
00667 ilast));
00668 a_ref(ilast, ilast - 1) = 0.;
00669 i__2 = ilast - ifrstm;
00670 template_blas_rot(&i__2, &a_ref(ifrstm, ilast), &c__1, &a_ref(ifrstm, ilast - 1),
00671 &c__1, &c__, &s);
00672 i__2 = ilast - ifrstm;
00673 template_blas_rot(&i__2, &b_ref(ifrstm, ilast), &c__1, &b_ref(ifrstm, ilast - 1),
00674 &c__1, &c__, &s);
00675 if (ilz) {
00676 template_blas_rot(n, &z___ref(1, ilast), &c__1, &z___ref(1, ilast - 1), &c__1,
00677 &c__, &s);
00678 }
00679
00680
00681
00682
00683 L80:
00684 if (b_ref(ilast, ilast) < 0.) {
00685 if (ilschr) {
00686 i__2 = ilast;
00687 for (j = ifrstm; j <= i__2; ++j) {
00688 a_ref(j, ilast) = -a_ref(j, ilast);
00689 b_ref(j, ilast) = -b_ref(j, ilast);
00690
00691 }
00692 } else {
00693 a_ref(ilast, ilast) = -a_ref(ilast, ilast);
00694 b_ref(ilast, ilast) = -b_ref(ilast, ilast);
00695 }
00696 if (ilz) {
00697 i__2 = *n;
00698 for (j = 1; j <= i__2; ++j) {
00699 z___ref(j, ilast) = -z___ref(j, ilast);
00700
00701 }
00702 }
00703 }
00704 alphar[ilast] = a_ref(ilast, ilast);
00705 alphai[ilast] = 0.;
00706 beta[ilast] = b_ref(ilast, ilast);
00707
00708
00709
00710 --ilast;
00711 if (ilast < *ilo) {
00712 goto L380;
00713 }
00714
00715
00716
00717 iiter = 0;
00718 eshift = 0.;
00719 if (! ilschr) {
00720 ilastm = ilast;
00721 if (ifrstm > ilast) {
00722 ifrstm = *ilo;
00723 }
00724 }
00725 goto L350;
00726
00727
00728
00729
00730
00731
00732 L110:
00733 ++iiter;
00734 if (! ilschr) {
00735 ifrstm = ifirst;
00736 }
00737
00738
00739
00740
00741
00742
00743
00744 if (iiter / 10 * 10 == iiter) {
00745
00746
00747
00748
00749 if ((Treal) maxit * safmin * (d__1 = a_ref(ilast - 1, ilast),
00750 absMACRO(d__1)) < (d__2 = b_ref(ilast - 1, ilast - 1), absMACRO(
00751 d__2))) {
00752 eshift += a_ref(ilast - 1, ilast) / b_ref(ilast - 1, ilast -
00753 1);
00754 } else {
00755 eshift += 1. / (safmin * (Treal) maxit);
00756 }
00757 s1 = 1.;
00758 wr = eshift;
00759
00760 } else {
00761
00762
00763
00764
00765
00766 d__1 = safmin * 100.;
00767 template_lapack_lag2(&a_ref(ilast - 1, ilast - 1), lda, &b_ref(ilast - 1, ilast
00768 - 1), ldb, &d__1, &s1, &s2, &wr, &wr2, &wi);
00769
00770
00771
00772 d__3 = 1., d__4 = absMACRO(wr), d__3 = maxMACRO(d__3,d__4), d__4 = absMACRO(wi);
00773 d__1 = s1, d__2 = safmin * maxMACRO(d__3,d__4);
00774 temp = maxMACRO(d__1,d__2);
00775 if (wi != 0.) {
00776 goto L200;
00777 }
00778 }
00779
00780
00781
00782 temp = minMACRO(ascale,1.) * (safmax * .5);
00783 if (s1 > temp) {
00784 scale = temp / s1;
00785 } else {
00786 scale = 1.;
00787 }
00788
00789 temp = minMACRO(bscale,1.) * (safmax * .5);
00790 if (absMACRO(wr) > temp) {
00791
00792 d__1 = scale, d__2 = temp / absMACRO(wr);
00793 scale = minMACRO(d__1,d__2);
00794 }
00795 s1 = scale * s1;
00796 wr = scale * wr;
00797
00798
00799
00800 i__2 = ifirst + 1;
00801 for (j = ilast - 1; j >= i__2; --j) {
00802 istart = j;
00803 temp = (d__1 = s1 * a_ref(j, j - 1), absMACRO(d__1));
00804 temp2 = (d__1 = s1 * a_ref(j, j) - wr * b_ref(j, j), absMACRO(d__1));
00805 tempr = maxMACRO(temp,temp2);
00806 if (tempr < 1. && tempr != 0.) {
00807 temp /= tempr;
00808 temp2 /= tempr;
00809 }
00810 if ((d__1 = ascale * a_ref(j + 1, j) * temp, absMACRO(d__1)) <= ascale
00811 * atol * temp2) {
00812 goto L130;
00813 }
00814
00815 }
00816
00817 istart = ifirst;
00818 L130:
00819
00820
00821
00822
00823
00824 temp = s1 * a_ref(istart, istart) - wr * b_ref(istart, istart);
00825 temp2 = s1 * a_ref(istart + 1, istart);
00826 template_lapack_lartg(&temp, &temp2, &c__, &s, &tempr);
00827
00828
00829
00830 i__2 = ilast - 1;
00831 for (j = istart; j <= i__2; ++j) {
00832 if (j > istart) {
00833 temp = a_ref(j, j - 1);
00834 template_lapack_lartg(&temp, &a_ref(j + 1, j - 1), &c__, &s, &a_ref(j, j -
00835 1));
00836 a_ref(j + 1, j - 1) = 0.;
00837 }
00838
00839 i__3 = ilastm;
00840 for (jc = j; jc <= i__3; ++jc) {
00841 temp = c__ * a_ref(j, jc) + s * a_ref(j + 1, jc);
00842 a_ref(j + 1, jc) = -s * a_ref(j, jc) + c__ * a_ref(j + 1, jc);
00843 a_ref(j, jc) = temp;
00844 temp2 = c__ * b_ref(j, jc) + s * b_ref(j + 1, jc);
00845 b_ref(j + 1, jc) = -s * b_ref(j, jc) + c__ * b_ref(j + 1, jc);
00846 b_ref(j, jc) = temp2;
00847
00848 }
00849 if (ilq) {
00850 i__3 = *n;
00851 for (jr = 1; jr <= i__3; ++jr) {
00852 temp = c__ * q_ref(jr, j) + s * q_ref(jr, j + 1);
00853 q_ref(jr, j + 1) = -s * q_ref(jr, j) + c__ * q_ref(jr, j
00854 + 1);
00855 q_ref(jr, j) = temp;
00856
00857 }
00858 }
00859
00860 temp = b_ref(j + 1, j + 1);
00861 template_lapack_lartg(&temp, &b_ref(j + 1, j), &c__, &s, &b_ref(j + 1, j + 1));
00862 b_ref(j + 1, j) = 0.;
00863
00864
00865 i__4 = j + 2;
00866 i__3 = minMACRO(i__4,ilast);
00867 for (jr = ifrstm; jr <= i__3; ++jr) {
00868 temp = c__ * a_ref(jr, j + 1) + s * a_ref(jr, j);
00869 a_ref(jr, j) = -s * a_ref(jr, j + 1) + c__ * a_ref(jr, j);
00870 a_ref(jr, j + 1) = temp;
00871
00872 }
00873 i__3 = j;
00874 for (jr = ifrstm; jr <= i__3; ++jr) {
00875 temp = c__ * b_ref(jr, j + 1) + s * b_ref(jr, j);
00876 b_ref(jr, j) = -s * b_ref(jr, j + 1) + c__ * b_ref(jr, j);
00877 b_ref(jr, j + 1) = temp;
00878
00879 }
00880 if (ilz) {
00881 i__3 = *n;
00882 for (jr = 1; jr <= i__3; ++jr) {
00883 temp = c__ * z___ref(jr, j + 1) + s * z___ref(jr, j);
00884 z___ref(jr, j) = -s * z___ref(jr, j + 1) + c__ * z___ref(
00885 jr, j);
00886 z___ref(jr, j + 1) = temp;
00887
00888 }
00889 }
00890
00891 }
00892
00893 goto L350;
00894
00895
00896
00897
00898
00899
00900
00901
00902 L200:
00903 if (ifirst + 1 == ilast) {
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913 template_lapack_lasv2(&b_ref(ilast - 1, ilast - 1), &b_ref(ilast - 1, ilast), &
00914 b_ref(ilast, ilast), &b22, &b11, &sr, &cr, &sl, &cl);
00915
00916 if (b11 < 0.) {
00917 cr = -cr;
00918 sr = -sr;
00919 b11 = -b11;
00920 b22 = -b22;
00921 }
00922
00923 i__2 = ilastm + 1 - ifirst;
00924 template_blas_rot(&i__2, &a_ref(ilast - 1, ilast - 1), lda, &a_ref(ilast,
00925 ilast - 1), lda, &cl, &sl);
00926 i__2 = ilast + 1 - ifrstm;
00927 template_blas_rot(&i__2, &a_ref(ifrstm, ilast - 1), &c__1, &a_ref(ifrstm,
00928 ilast), &c__1, &cr, &sr);
00929
00930 if (ilast < ilastm) {
00931 i__2 = ilastm - ilast;
00932 template_blas_rot(&i__2, &b_ref(ilast - 1, ilast + 1), ldb, &b_ref(ilast,
00933 ilast + 1), lda, &cl, &sl);
00934 }
00935 if (ifrstm < ilast - 1) {
00936 i__2 = ifirst - ifrstm;
00937 template_blas_rot(&i__2, &b_ref(ifrstm, ilast - 1), &c__1, &b_ref(ifrstm,
00938 ilast), &c__1, &cr, &sr);
00939 }
00940
00941 if (ilq) {
00942 template_blas_rot(n, &q_ref(1, ilast - 1), &c__1, &q_ref(1, ilast), &c__1,
00943 &cl, &sl);
00944 }
00945 if (ilz) {
00946 template_blas_rot(n, &z___ref(1, ilast - 1), &c__1, &z___ref(1, ilast), &
00947 c__1, &cr, &sr);
00948 }
00949
00950 b_ref(ilast - 1, ilast - 1) = b11;
00951 b_ref(ilast - 1, ilast) = 0.;
00952 b_ref(ilast, ilast - 1) = 0.;
00953 b_ref(ilast, ilast) = b22;
00954
00955
00956
00957 if (b22 < 0.) {
00958 i__2 = ilast;
00959 for (j = ifrstm; j <= i__2; ++j) {
00960 a_ref(j, ilast) = -a_ref(j, ilast);
00961 b_ref(j, ilast) = -b_ref(j, ilast);
00962
00963 }
00964
00965 if (ilz) {
00966 i__2 = *n;
00967 for (j = 1; j <= i__2; ++j) {
00968 z___ref(j, ilast) = -z___ref(j, ilast);
00969
00970 }
00971 }
00972 }
00973
00974
00975
00976
00977
00978 d__1 = safmin * 100.;
00979 template_lapack_lag2(&a_ref(ilast - 1, ilast - 1), lda, &b_ref(ilast - 1, ilast
00980 - 1), ldb, &d__1, &s1, &temp, &wr, &temp2, &wi);
00981
00982
00983
00984
00985 if (wi == 0.) {
00986 goto L350;
00987 }
00988 s1inv = 1. / s1;
00989
00990
00991
00992 a11 = a_ref(ilast - 1, ilast - 1);
00993 a21 = a_ref(ilast, ilast - 1);
00994 a12 = a_ref(ilast - 1, ilast);
00995 a22 = a_ref(ilast, ilast);
00996
00997
00998
00999
01000
01001
01002
01003 c11r = s1 * a11 - wr * b11;
01004 c11i = -wi * b11;
01005 c12 = s1 * a12;
01006 c21 = s1 * a21;
01007 c22r = s1 * a22 - wr * b22;
01008 c22i = -wi * b22;
01009
01010 if (absMACRO(c11r) + absMACRO(c11i) + absMACRO(c12) > absMACRO(c21) + absMACRO(c22r) + absMACRO(
01011 c22i)) {
01012 t = template_lapack_lapy3(&c12, &c11r, &c11i);
01013 cz = c12 / t;
01014 szr = -c11r / t;
01015 szi = -c11i / t;
01016 } else {
01017 cz = template_lapack_lapy2(&c22r, &c22i);
01018 if (cz <= safmin) {
01019 cz = 0.;
01020 szr = 1.;
01021 szi = 0.;
01022 } else {
01023 tempr = c22r / cz;
01024 tempi = c22i / cz;
01025 t = template_lapack_lapy2(&cz, &c21);
01026 cz /= t;
01027 szr = -c21 * tempr / t;
01028 szi = c21 * tempi / t;
01029 }
01030 }
01031
01032
01033
01034
01035
01036
01037
01038 an = absMACRO(a11) + absMACRO(a12) + absMACRO(a21) + absMACRO(a22);
01039 bn = absMACRO(b11) + absMACRO(b22);
01040 wabs = absMACRO(wr) + absMACRO(wi);
01041 if (s1 * an > wabs * bn) {
01042 cq = cz * b11;
01043 sqr = szr * b22;
01044 sqi = -szi * b22;
01045 } else {
01046 a1r = cz * a11 + szr * a12;
01047 a1i = szi * a12;
01048 a2r = cz * a21 + szr * a22;
01049 a2i = szi * a22;
01050 cq = template_lapack_lapy2(&a1r, &a1i);
01051 if (cq <= safmin) {
01052 cq = 0.;
01053 sqr = 1.;
01054 sqi = 0.;
01055 } else {
01056 tempr = a1r / cq;
01057 tempi = a1i / cq;
01058 sqr = tempr * a2r + tempi * a2i;
01059 sqi = tempi * a2r - tempr * a2i;
01060 }
01061 }
01062 t = template_lapack_lapy3(&cq, &sqr, &sqi);
01063 cq /= t;
01064 sqr /= t;
01065 sqi /= t;
01066
01067
01068
01069 tempr = sqr * szr - sqi * szi;
01070 tempi = sqr * szi + sqi * szr;
01071 b1r = cq * cz * b11 + tempr * b22;
01072 b1i = tempi * b22;
01073 b1a = template_lapack_lapy2(&b1r, &b1i);
01074 b2r = cq * cz * b22 + tempr * b11;
01075 b2i = -tempi * b11;
01076 b2a = template_lapack_lapy2(&b2r, &b2i);
01077
01078
01079
01080 beta[ilast - 1] = b1a;
01081 beta[ilast] = b2a;
01082 alphar[ilast - 1] = wr * b1a * s1inv;
01083 alphai[ilast - 1] = wi * b1a * s1inv;
01084 alphar[ilast] = wr * b2a * s1inv;
01085 alphai[ilast] = -(wi * b2a) * s1inv;
01086
01087
01088
01089 ilast = ifirst - 1;
01090 if (ilast < *ilo) {
01091 goto L380;
01092 }
01093
01094
01095
01096 iiter = 0;
01097 eshift = 0.;
01098 if (! ilschr) {
01099 ilastm = ilast;
01100 if (ifrstm > ilast) {
01101 ifrstm = *ilo;
01102 }
01103 }
01104 goto L350;
01105 } else {
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119 ad11 = ascale * a_ref(ilast - 1, ilast - 1) / (bscale * b_ref(
01120 ilast - 1, ilast - 1));
01121 ad21 = ascale * a_ref(ilast, ilast - 1) / (bscale * b_ref(ilast -
01122 1, ilast - 1));
01123 ad12 = ascale * a_ref(ilast - 1, ilast) / (bscale * b_ref(ilast,
01124 ilast));
01125 ad22 = ascale * a_ref(ilast, ilast) / (bscale * b_ref(ilast,
01126 ilast));
01127 u12 = b_ref(ilast - 1, ilast) / b_ref(ilast, ilast);
01128 ad11l = ascale * a_ref(ifirst, ifirst) / (bscale * b_ref(ifirst,
01129 ifirst));
01130 ad21l = ascale * a_ref(ifirst + 1, ifirst) / (bscale * b_ref(
01131 ifirst, ifirst));
01132 ad12l = ascale * a_ref(ifirst, ifirst + 1) / (bscale * b_ref(
01133 ifirst + 1, ifirst + 1));
01134 ad22l = ascale * a_ref(ifirst + 1, ifirst + 1) / (bscale * b_ref(
01135 ifirst + 1, ifirst + 1));
01136 ad32l = ascale * a_ref(ifirst + 2, ifirst + 1) / (bscale * b_ref(
01137 ifirst + 1, ifirst + 1));
01138 u12l = b_ref(ifirst, ifirst + 1) / b_ref(ifirst + 1, ifirst + 1);
01139
01140 v[0] = (ad11 - ad11l) * (ad22 - ad11l) - ad12 * ad21 + ad21 * u12
01141 * ad11l + (ad12l - ad11l * u12l) * ad21l;
01142 v[1] = (ad22l - ad11l - ad21l * u12l - (ad11 - ad11l) - (ad22 -
01143 ad11l) + ad21 * u12) * ad21l;
01144 v[2] = ad32l * ad21l;
01145
01146 istart = ifirst;
01147
01148 template_lapack_larfg(&c__3, v, &v[1], &c__1, &tau);
01149 v[0] = 1.;
01150
01151
01152
01153 i__2 = ilast - 2;
01154 for (j = istart; j <= i__2; ++j) {
01155
01156
01157
01158
01159
01160 if (j > istart) {
01161 v[0] = a_ref(j, j - 1);
01162 v[1] = a_ref(j + 1, j - 1);
01163 v[2] = a_ref(j + 2, j - 1);
01164
01165 template_lapack_larfg(&c__3, &a_ref(j, j - 1), &v[1], &c__1, &tau);
01166 v[0] = 1.;
01167 a_ref(j + 1, j - 1) = 0.;
01168 a_ref(j + 2, j - 1) = 0.;
01169 }
01170
01171 i__3 = ilastm;
01172 for (jc = j; jc <= i__3; ++jc) {
01173 temp = tau * (a_ref(j, jc) + v[1] * a_ref(j + 1, jc) + v[
01174 2] * a_ref(j + 2, jc));
01175 a_ref(j, jc) = a_ref(j, jc) - temp;
01176 a_ref(j + 1, jc) = a_ref(j + 1, jc) - temp * v[1];
01177 a_ref(j + 2, jc) = a_ref(j + 2, jc) - temp * v[2];
01178 temp2 = tau * (b_ref(j, jc) + v[1] * b_ref(j + 1, jc) + v[
01179 2] * b_ref(j + 2, jc));
01180 b_ref(j, jc) = b_ref(j, jc) - temp2;
01181 b_ref(j + 1, jc) = b_ref(j + 1, jc) - temp2 * v[1];
01182 b_ref(j + 2, jc) = b_ref(j + 2, jc) - temp2 * v[2];
01183
01184 }
01185 if (ilq) {
01186 i__3 = *n;
01187 for (jr = 1; jr <= i__3; ++jr) {
01188 temp = tau * (q_ref(jr, j) + v[1] * q_ref(jr, j + 1)
01189 + v[2] * q_ref(jr, j + 2));
01190 q_ref(jr, j) = q_ref(jr, j) - temp;
01191 q_ref(jr, j + 1) = q_ref(jr, j + 1) - temp * v[1];
01192 q_ref(jr, j + 2) = q_ref(jr, j + 2) - temp * v[2];
01193
01194 }
01195 }
01196
01197
01198
01199
01200
01201 ilpivt = FALSE_;
01202
01203 d__3 = (d__1 = b_ref(j + 1, j + 1), absMACRO(d__1)), d__4 = (d__2 =
01204 b_ref(j + 1, j + 2), absMACRO(d__2));
01205 temp = maxMACRO(d__3,d__4);
01206
01207 d__3 = (d__1 = b_ref(j + 2, j + 1), absMACRO(d__1)), d__4 = (d__2 =
01208 b_ref(j + 2, j + 2), absMACRO(d__2));
01209 temp2 = maxMACRO(d__3,d__4);
01210 if (maxMACRO(temp,temp2) < safmin) {
01211 scale = 0.;
01212 u1 = 1.;
01213 u2 = 0.;
01214 goto L250;
01215 } else if (temp >= temp2) {
01216 w11 = b_ref(j + 1, j + 1);
01217 w21 = b_ref(j + 2, j + 1);
01218 w12 = b_ref(j + 1, j + 2);
01219 w22 = b_ref(j + 2, j + 2);
01220 u1 = b_ref(j + 1, j);
01221 u2 = b_ref(j + 2, j);
01222 } else {
01223 w21 = b_ref(j + 1, j + 1);
01224 w11 = b_ref(j + 2, j + 1);
01225 w22 = b_ref(j + 1, j + 2);
01226 w12 = b_ref(j + 2, j + 2);
01227 u2 = b_ref(j + 1, j);
01228 u1 = b_ref(j + 2, j);
01229 }
01230
01231
01232
01233 if (absMACRO(w12) > absMACRO(w11)) {
01234 ilpivt = TRUE_;
01235 temp = w12;
01236 temp2 = w22;
01237 w12 = w11;
01238 w22 = w21;
01239 w11 = temp;
01240 w21 = temp2;
01241 }
01242
01243
01244
01245 temp = w21 / w11;
01246 u2 -= temp * u1;
01247 w22 -= temp * w12;
01248 w21 = 0.;
01249
01250
01251
01252 scale = 1.;
01253 if (absMACRO(w22) < safmin) {
01254 scale = 0.;
01255 u2 = 1.;
01256 u1 = -w12 / w11;
01257 goto L250;
01258 }
01259 if (absMACRO(w22) < absMACRO(u2)) {
01260 scale = (d__1 = w22 / u2, absMACRO(d__1));
01261 }
01262 if (absMACRO(w11) < absMACRO(u1)) {
01263
01264 d__2 = scale, d__3 = (d__1 = w11 / u1, absMACRO(d__1));
01265 scale = minMACRO(d__2,d__3);
01266 }
01267
01268
01269
01270 u2 = scale * u2 / w22;
01271 u1 = (scale * u1 - w12 * u2) / w11;
01272
01273 L250:
01274 if (ilpivt) {
01275 temp = u2;
01276 u2 = u1;
01277 u1 = temp;
01278 }
01279
01280
01281
01282
01283 d__1 = scale;
01284
01285 d__2 = u1;
01286
01287 d__3 = u2;
01288 t = template_blas_sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
01289 tau = scale / t + 1.;
01290 vs = -1. / (scale + t);
01291 v[0] = 1.;
01292 v[1] = vs * u1;
01293 v[2] = vs * u2;
01294
01295
01296
01297
01298 i__4 = j + 3;
01299 i__3 = minMACRO(i__4,ilast);
01300 for (jr = ifrstm; jr <= i__3; ++jr) {
01301 temp = tau * (a_ref(jr, j) + v[1] * a_ref(jr, j + 1) + v[
01302 2] * a_ref(jr, j + 2));
01303 a_ref(jr, j) = a_ref(jr, j) - temp;
01304 a_ref(jr, j + 1) = a_ref(jr, j + 1) - temp * v[1];
01305 a_ref(jr, j + 2) = a_ref(jr, j + 2) - temp * v[2];
01306
01307 }
01308 i__3 = j + 2;
01309 for (jr = ifrstm; jr <= i__3; ++jr) {
01310 temp = tau * (b_ref(jr, j) + v[1] * b_ref(jr, j + 1) + v[
01311 2] * b_ref(jr, j + 2));
01312 b_ref(jr, j) = b_ref(jr, j) - temp;
01313 b_ref(jr, j + 1) = b_ref(jr, j + 1) - temp * v[1];
01314 b_ref(jr, j + 2) = b_ref(jr, j + 2) - temp * v[2];
01315
01316 }
01317 if (ilz) {
01318 i__3 = *n;
01319 for (jr = 1; jr <= i__3; ++jr) {
01320 temp = tau * (z___ref(jr, j) + v[1] * z___ref(jr, j +
01321 1) + v[2] * z___ref(jr, j + 2));
01322 z___ref(jr, j) = z___ref(jr, j) - temp;
01323 z___ref(jr, j + 1) = z___ref(jr, j + 1) - temp * v[1];
01324 z___ref(jr, j + 2) = z___ref(jr, j + 2) - temp * v[2];
01325
01326 }
01327 }
01328 b_ref(j + 1, j) = 0.;
01329 b_ref(j + 2, j) = 0.;
01330
01331 }
01332
01333
01334
01335
01336
01337 j = ilast - 1;
01338 temp = a_ref(j, j - 1);
01339 template_lapack_lartg(&temp, &a_ref(j + 1, j - 1), &c__, &s, &a_ref(j, j - 1));
01340 a_ref(j + 1, j - 1) = 0.;
01341
01342 i__2 = ilastm;
01343 for (jc = j; jc <= i__2; ++jc) {
01344 temp = c__ * a_ref(j, jc) + s * a_ref(j + 1, jc);
01345 a_ref(j + 1, jc) = -s * a_ref(j, jc) + c__ * a_ref(j + 1, jc);
01346 a_ref(j, jc) = temp;
01347 temp2 = c__ * b_ref(j, jc) + s * b_ref(j + 1, jc);
01348 b_ref(j + 1, jc) = -s * b_ref(j, jc) + c__ * b_ref(j + 1, jc);
01349 b_ref(j, jc) = temp2;
01350
01351 }
01352 if (ilq) {
01353 i__2 = *n;
01354 for (jr = 1; jr <= i__2; ++jr) {
01355 temp = c__ * q_ref(jr, j) + s * q_ref(jr, j + 1);
01356 q_ref(jr, j + 1) = -s * q_ref(jr, j) + c__ * q_ref(jr, j
01357 + 1);
01358 q_ref(jr, j) = temp;
01359
01360 }
01361 }
01362
01363
01364
01365 temp = b_ref(j + 1, j + 1);
01366 template_lapack_lartg(&temp, &b_ref(j + 1, j), &c__, &s, &b_ref(j + 1, j + 1));
01367 b_ref(j + 1, j) = 0.;
01368
01369 i__2 = ilast;
01370 for (jr = ifrstm; jr <= i__2; ++jr) {
01371 temp = c__ * a_ref(jr, j + 1) + s * a_ref(jr, j);
01372 a_ref(jr, j) = -s * a_ref(jr, j + 1) + c__ * a_ref(jr, j);
01373 a_ref(jr, j + 1) = temp;
01374
01375 }
01376 i__2 = ilast - 1;
01377 for (jr = ifrstm; jr <= i__2; ++jr) {
01378 temp = c__ * b_ref(jr, j + 1) + s * b_ref(jr, j);
01379 b_ref(jr, j) = -s * b_ref(jr, j + 1) + c__ * b_ref(jr, j);
01380 b_ref(jr, j + 1) = temp;
01381
01382 }
01383 if (ilz) {
01384 i__2 = *n;
01385 for (jr = 1; jr <= i__2; ++jr) {
01386 temp = c__ * z___ref(jr, j + 1) + s * z___ref(jr, j);
01387 z___ref(jr, j) = -s * z___ref(jr, j + 1) + c__ * z___ref(
01388 jr, j);
01389 z___ref(jr, j + 1) = temp;
01390
01391 }
01392 }
01393
01394
01395
01396 }
01397
01398 goto L350;
01399
01400
01401
01402 L350:
01403
01404 ;
01405 }
01406
01407
01408
01409
01410 *info = ilast;
01411 goto L420;
01412
01413
01414
01415 L380:
01416
01417
01418
01419 i__1 = *ilo - 1;
01420 for (j = 1; j <= i__1; ++j) {
01421 if (b_ref(j, j) < 0.) {
01422 if (ilschr) {
01423 i__2 = j;
01424 for (jr = 1; jr <= i__2; ++jr) {
01425 a_ref(jr, j) = -a_ref(jr, j);
01426 b_ref(jr, j) = -b_ref(jr, j);
01427
01428 }
01429 } else {
01430 a_ref(j, j) = -a_ref(j, j);
01431 b_ref(j, j) = -b_ref(j, j);
01432 }
01433 if (ilz) {
01434 i__2 = *n;
01435 for (jr = 1; jr <= i__2; ++jr) {
01436 z___ref(jr, j) = -z___ref(jr, j);
01437
01438 }
01439 }
01440 }
01441 alphar[j] = a_ref(j, j);
01442 alphai[j] = 0.;
01443 beta[j] = b_ref(j, j);
01444
01445 }
01446
01447
01448
01449 *info = 0;
01450
01451
01452
01453 L420:
01454 work[1] = (Treal) (*n);
01455 return 0;
01456
01457
01458
01459 }
01460
01461 #undef z___ref
01462 #undef q_ref
01463 #undef b_ref
01464 #undef a_ref
01465
01466
01467 #endif