00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015
00016
00017
00018 static integer c__1 = 1;
00019 static integer c__0 = 0;
00020 static integer c_n1 = -1;
00021
00022 int cgeevx_(char *balanc, char *jobvl, char *jobvr, char *
00023 sense, integer *n, complex *a, integer *lda, complex *w, complex *vl,
00024 integer *ldvl, complex *vr, integer *ldvr, integer *ilo, integer *ihi,
00025 real *scale, real *abnrm, real *rconde, real *rcondv, complex *work,
00026 integer *lwork, real *rwork, integer *info)
00027 {
00028
00029 integer a_dim1, a_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1,
00030 i__2, i__3;
00031 real r__1, r__2;
00032 complex q__1, q__2;
00033
00034
00035 double sqrt(doublereal), r_imag(complex *);
00036 void r_cnjg(complex *, complex *);
00037
00038
00039 integer i__, k;
00040 char job[1];
00041 real scl, dum[1], eps;
00042 complex tmp;
00043 char side[1];
00044 real anrm;
00045 integer ierr, itau, iwrk, nout;
00046 extern int cscal_(integer *, complex *, complex *,
00047 integer *);
00048 integer icond;
00049 extern logical lsame_(char *, char *);
00050 extern doublereal scnrm2_(integer *, complex *, integer *);
00051 extern int cgebak_(char *, char *, integer *, integer *,
00052 integer *, real *, integer *, complex *, integer *, integer *), cgebal_(char *, integer *, complex *, integer *,
00053 integer *, integer *, real *, integer *), slabad_(real *,
00054 real *);
00055 logical scalea;
00056 extern doublereal clange_(char *, integer *, integer *, complex *,
00057 integer *, real *);
00058 real cscale;
00059 extern int cgehrd_(integer *, integer *, integer *,
00060 complex *, integer *, complex *, complex *, integer *, integer *),
00061 clascl_(char *, integer *, integer *, real *, real *, integer *,
00062 integer *, complex *, integer *, integer *);
00063 extern doublereal slamch_(char *);
00064 extern int csscal_(integer *, real *, complex *, integer
00065 *), clacpy_(char *, integer *, integer *, complex *, integer *,
00066 complex *, integer *), xerbla_(char *, integer *);
00067 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
00068 integer *, integer *);
00069 logical select[1];
00070 real bignum;
00071 extern int slascl_(char *, integer *, integer *, real *,
00072 real *, integer *, integer *, real *, integer *, integer *);
00073 extern integer isamax_(integer *, real *, integer *);
00074 extern int chseqr_(char *, char *, integer *, integer *,
00075 integer *, complex *, integer *, complex *, complex *, integer *,
00076 complex *, integer *, integer *), ctrevc_(char *,
00077 char *, logical *, integer *, complex *, integer *, complex *,
00078 integer *, complex *, integer *, integer *, integer *, complex *,
00079 real *, integer *), cunghr_(integer *, integer *,
00080 integer *, complex *, integer *, complex *, complex *, integer *,
00081 integer *), ctrsna_(char *, char *, logical *, integer *, complex
00082 *, integer *, complex *, integer *, complex *, integer *, real *,
00083 real *, integer *, integer *, complex *, integer *, real *,
00084 integer *);
00085 integer minwrk, maxwrk;
00086 logical wantvl, wntsnb;
00087 integer hswork;
00088 logical wntsne;
00089 real smlnum;
00090 logical lquery, wantvr, wntsnn, wntsnv;
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
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284 a_dim1 = *lda;
00285 a_offset = 1 + a_dim1;
00286 a -= a_offset;
00287 --w;
00288 vl_dim1 = *ldvl;
00289 vl_offset = 1 + vl_dim1;
00290 vl -= vl_offset;
00291 vr_dim1 = *ldvr;
00292 vr_offset = 1 + vr_dim1;
00293 vr -= vr_offset;
00294 --scale;
00295 --rconde;
00296 --rcondv;
00297 --work;
00298 --rwork;
00299
00300
00301 *info = 0;
00302 lquery = *lwork == -1;
00303 wantvl = lsame_(jobvl, "V");
00304 wantvr = lsame_(jobvr, "V");
00305 wntsnn = lsame_(sense, "N");
00306 wntsne = lsame_(sense, "E");
00307 wntsnv = lsame_(sense, "V");
00308 wntsnb = lsame_(sense, "B");
00309 if (! (lsame_(balanc, "N") || lsame_(balanc, "S") || lsame_(balanc, "P")
00310 || lsame_(balanc, "B"))) {
00311 *info = -1;
00312 } else if (! wantvl && ! lsame_(jobvl, "N")) {
00313 *info = -2;
00314 } else if (! wantvr && ! lsame_(jobvr, "N")) {
00315 *info = -3;
00316 } else if (! (wntsnn || wntsne || wntsnb || wntsnv) || (wntsne || wntsnb)
00317 && ! (wantvl && wantvr)) {
00318 *info = -4;
00319 } else if (*n < 0) {
00320 *info = -5;
00321 } else if (*lda < max(1,*n)) {
00322 *info = -7;
00323 } else if (*ldvl < 1 || wantvl && *ldvl < *n) {
00324 *info = -10;
00325 } else if (*ldvr < 1 || wantvr && *ldvr < *n) {
00326 *info = -12;
00327 }
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340 if (*info == 0) {
00341 if (*n == 0) {
00342 minwrk = 1;
00343 maxwrk = 1;
00344 } else {
00345 maxwrk = *n + *n * ilaenv_(&c__1, "CGEHRD", " ", n, &c__1, n, &
00346 c__0);
00347
00348 if (wantvl) {
00349 chseqr_("S", "V", n, &c__1, n, &a[a_offset], lda, &w[1], &vl[
00350 vl_offset], ldvl, &work[1], &c_n1, info);
00351 } else if (wantvr) {
00352 chseqr_("S", "V", n, &c__1, n, &a[a_offset], lda, &w[1], &vr[
00353 vr_offset], ldvr, &work[1], &c_n1, info);
00354 } else {
00355 if (wntsnn) {
00356 chseqr_("E", "N", n, &c__1, n, &a[a_offset], lda, &w[1], &
00357 vr[vr_offset], ldvr, &work[1], &c_n1, info);
00358 } else {
00359 chseqr_("S", "N", n, &c__1, n, &a[a_offset], lda, &w[1], &
00360 vr[vr_offset], ldvr, &work[1], &c_n1, info);
00361 }
00362 }
00363 hswork = work[1].r;
00364
00365 if (! wantvl && ! wantvr) {
00366 minwrk = *n << 1;
00367 if (! (wntsnn || wntsne)) {
00368
00369 i__1 = minwrk, i__2 = *n * *n + (*n << 1);
00370 minwrk = max(i__1,i__2);
00371 }
00372 maxwrk = max(maxwrk,hswork);
00373 if (! (wntsnn || wntsne)) {
00374
00375 i__1 = maxwrk, i__2 = *n * *n + (*n << 1);
00376 maxwrk = max(i__1,i__2);
00377 }
00378 } else {
00379 minwrk = *n << 1;
00380 if (! (wntsnn || wntsne)) {
00381
00382 i__1 = minwrk, i__2 = *n * *n + (*n << 1);
00383 minwrk = max(i__1,i__2);
00384 }
00385 maxwrk = max(maxwrk,hswork);
00386
00387 i__1 = maxwrk, i__2 = *n + (*n - 1) * ilaenv_(&c__1, "CUNGHR",
00388 " ", n, &c__1, n, &c_n1);
00389 maxwrk = max(i__1,i__2);
00390 if (! (wntsnn || wntsne)) {
00391
00392 i__1 = maxwrk, i__2 = *n * *n + (*n << 1);
00393 maxwrk = max(i__1,i__2);
00394 }
00395
00396 i__1 = maxwrk, i__2 = *n << 1;
00397 maxwrk = max(i__1,i__2);
00398 }
00399 maxwrk = max(maxwrk,minwrk);
00400 }
00401 work[1].r = (real) maxwrk, work[1].i = 0.f;
00402
00403 if (*lwork < minwrk && ! lquery) {
00404 *info = -20;
00405 }
00406 }
00407
00408 if (*info != 0) {
00409 i__1 = -(*info);
00410 xerbla_("CGEEVX", &i__1);
00411 return 0;
00412 } else if (lquery) {
00413 return 0;
00414 }
00415
00416
00417
00418 if (*n == 0) {
00419 return 0;
00420 }
00421
00422
00423
00424 eps = slamch_("P");
00425 smlnum = slamch_("S");
00426 bignum = 1.f / smlnum;
00427 slabad_(&smlnum, &bignum);
00428 smlnum = sqrt(smlnum) / eps;
00429 bignum = 1.f / smlnum;
00430
00431
00432
00433 icond = 0;
00434 anrm = clange_("M", n, n, &a[a_offset], lda, dum);
00435 scalea = FALSE_;
00436 if (anrm > 0.f && anrm < smlnum) {
00437 scalea = TRUE_;
00438 cscale = smlnum;
00439 } else if (anrm > bignum) {
00440 scalea = TRUE_;
00441 cscale = bignum;
00442 }
00443 if (scalea) {
00444 clascl_("G", &c__0, &c__0, &anrm, &cscale, n, n, &a[a_offset], lda, &
00445 ierr);
00446 }
00447
00448
00449
00450 cgebal_(balanc, n, &a[a_offset], lda, ilo, ihi, &scale[1], &ierr);
00451 *abnrm = clange_("1", n, n, &a[a_offset], lda, dum);
00452 if (scalea) {
00453 dum[0] = *abnrm;
00454 slascl_("G", &c__0, &c__0, &cscale, &anrm, &c__1, &c__1, dum, &c__1, &
00455 ierr);
00456 *abnrm = dum[0];
00457 }
00458
00459
00460
00461
00462
00463 itau = 1;
00464 iwrk = itau + *n;
00465 i__1 = *lwork - iwrk + 1;
00466 cgehrd_(n, ilo, ihi, &a[a_offset], lda, &work[itau], &work[iwrk], &i__1, &
00467 ierr);
00468
00469 if (wantvl) {
00470
00471
00472
00473
00474 *(unsigned char *)side = 'L';
00475 clacpy_("L", n, n, &a[a_offset], lda, &vl[vl_offset], ldvl)
00476 ;
00477
00478
00479
00480
00481
00482 i__1 = *lwork - iwrk + 1;
00483 cunghr_(n, ilo, ihi, &vl[vl_offset], ldvl, &work[itau], &work[iwrk], &
00484 i__1, &ierr);
00485
00486
00487
00488
00489
00490 iwrk = itau;
00491 i__1 = *lwork - iwrk + 1;
00492 chseqr_("S", "V", n, ilo, ihi, &a[a_offset], lda, &w[1], &vl[
00493 vl_offset], ldvl, &work[iwrk], &i__1, info);
00494
00495 if (wantvr) {
00496
00497
00498
00499
00500 *(unsigned char *)side = 'B';
00501 clacpy_("F", n, n, &vl[vl_offset], ldvl, &vr[vr_offset], ldvr);
00502 }
00503
00504 } else if (wantvr) {
00505
00506
00507
00508
00509 *(unsigned char *)side = 'R';
00510 clacpy_("L", n, n, &a[a_offset], lda, &vr[vr_offset], ldvr)
00511 ;
00512
00513
00514
00515
00516
00517 i__1 = *lwork - iwrk + 1;
00518 cunghr_(n, ilo, ihi, &vr[vr_offset], ldvr, &work[itau], &work[iwrk], &
00519 i__1, &ierr);
00520
00521
00522
00523
00524
00525 iwrk = itau;
00526 i__1 = *lwork - iwrk + 1;
00527 chseqr_("S", "V", n, ilo, ihi, &a[a_offset], lda, &w[1], &vr[
00528 vr_offset], ldvr, &work[iwrk], &i__1, info);
00529
00530 } else {
00531
00532
00533
00534
00535 if (wntsnn) {
00536 *(unsigned char *)job = 'E';
00537 } else {
00538 *(unsigned char *)job = 'S';
00539 }
00540
00541
00542
00543
00544 iwrk = itau;
00545 i__1 = *lwork - iwrk + 1;
00546 chseqr_(job, "N", n, ilo, ihi, &a[a_offset], lda, &w[1], &vr[
00547 vr_offset], ldvr, &work[iwrk], &i__1, info);
00548 }
00549
00550
00551
00552 if (*info > 0) {
00553 goto L50;
00554 }
00555
00556 if (wantvl || wantvr) {
00557
00558
00559
00560
00561
00562 ctrevc_(side, "B", select, n, &a[a_offset], lda, &vl[vl_offset], ldvl,
00563 &vr[vr_offset], ldvr, n, &nout, &work[iwrk], &rwork[1], &
00564 ierr);
00565 }
00566
00567
00568
00569
00570
00571 if (! wntsnn) {
00572 ctrsna_(sense, "A", select, n, &a[a_offset], lda, &vl[vl_offset],
00573 ldvl, &vr[vr_offset], ldvr, &rconde[1], &rcondv[1], n, &nout,
00574 &work[iwrk], n, &rwork[1], &icond);
00575 }
00576
00577 if (wantvl) {
00578
00579
00580
00581 cgebak_(balanc, "L", n, ilo, ihi, &scale[1], n, &vl[vl_offset], ldvl,
00582 &ierr);
00583
00584
00585
00586 i__1 = *n;
00587 for (i__ = 1; i__ <= i__1; ++i__) {
00588 scl = 1.f / scnrm2_(n, &vl[i__ * vl_dim1 + 1], &c__1);
00589 csscal_(n, &scl, &vl[i__ * vl_dim1 + 1], &c__1);
00590 i__2 = *n;
00591 for (k = 1; k <= i__2; ++k) {
00592 i__3 = k + i__ * vl_dim1;
00593
00594 r__1 = vl[i__3].r;
00595
00596 r__2 = r_imag(&vl[k + i__ * vl_dim1]);
00597 rwork[k] = r__1 * r__1 + r__2 * r__2;
00598
00599 }
00600 k = isamax_(n, &rwork[1], &c__1);
00601 r_cnjg(&q__2, &vl[k + i__ * vl_dim1]);
00602 r__1 = sqrt(rwork[k]);
00603 q__1.r = q__2.r / r__1, q__1.i = q__2.i / r__1;
00604 tmp.r = q__1.r, tmp.i = q__1.i;
00605 cscal_(n, &tmp, &vl[i__ * vl_dim1 + 1], &c__1);
00606 i__2 = k + i__ * vl_dim1;
00607 i__3 = k + i__ * vl_dim1;
00608 r__1 = vl[i__3].r;
00609 q__1.r = r__1, q__1.i = 0.f;
00610 vl[i__2].r = q__1.r, vl[i__2].i = q__1.i;
00611
00612 }
00613 }
00614
00615 if (wantvr) {
00616
00617
00618
00619 cgebak_(balanc, "R", n, ilo, ihi, &scale[1], n, &vr[vr_offset], ldvr,
00620 &ierr);
00621
00622
00623
00624 i__1 = *n;
00625 for (i__ = 1; i__ <= i__1; ++i__) {
00626 scl = 1.f / scnrm2_(n, &vr[i__ * vr_dim1 + 1], &c__1);
00627 csscal_(n, &scl, &vr[i__ * vr_dim1 + 1], &c__1);
00628 i__2 = *n;
00629 for (k = 1; k <= i__2; ++k) {
00630 i__3 = k + i__ * vr_dim1;
00631
00632 r__1 = vr[i__3].r;
00633
00634 r__2 = r_imag(&vr[k + i__ * vr_dim1]);
00635 rwork[k] = r__1 * r__1 + r__2 * r__2;
00636
00637 }
00638 k = isamax_(n, &rwork[1], &c__1);
00639 r_cnjg(&q__2, &vr[k + i__ * vr_dim1]);
00640 r__1 = sqrt(rwork[k]);
00641 q__1.r = q__2.r / r__1, q__1.i = q__2.i / r__1;
00642 tmp.r = q__1.r, tmp.i = q__1.i;
00643 cscal_(n, &tmp, &vr[i__ * vr_dim1 + 1], &c__1);
00644 i__2 = k + i__ * vr_dim1;
00645 i__3 = k + i__ * vr_dim1;
00646 r__1 = vr[i__3].r;
00647 q__1.r = r__1, q__1.i = 0.f;
00648 vr[i__2].r = q__1.r, vr[i__2].i = q__1.i;
00649
00650 }
00651 }
00652
00653
00654
00655 L50:
00656 if (scalea) {
00657 i__1 = *n - *info;
00658
00659 i__3 = *n - *info;
00660 i__2 = max(i__3,1);
00661 clascl_("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &w[*info + 1]
00662 , &i__2, &ierr);
00663 if (*info == 0) {
00664 if ((wntsnv || wntsnb) && icond == 0) {
00665 slascl_("G", &c__0, &c__0, &cscale, &anrm, n, &c__1, &rcondv[
00666 1], n, &ierr);
00667 }
00668 } else {
00669 i__1 = *ilo - 1;
00670 clascl_("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &w[1], n,
00671 &ierr);
00672 }
00673 }
00674
00675 work[1].r = (real) maxwrk, work[1].i = 0.f;
00676 return 0;
00677
00678
00679
00680 }