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