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