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