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 dgeesx_(char *jobvs, char *sort, L_fp select, char *
00023 sense, integer *n, doublereal *a, integer *lda, integer *sdim,
00024 doublereal *wr, doublereal *wi, doublereal *vs, integer *ldvs,
00025 doublereal *rconde, doublereal *rcondv, doublereal *work, integer *
00026 lwork, integer *iwork, integer *liwork, logical *bwork, integer *info)
00027 {
00028
00029 integer a_dim1, a_offset, vs_dim1, vs_offset, i__1, i__2, i__3;
00030
00031
00032 double sqrt(doublereal);
00033
00034
00035 integer i__, i1, i2, ip, ihi, ilo;
00036 doublereal dum[1], eps;
00037 integer ibal;
00038 doublereal anrm;
00039 integer ierr, itau, iwrk, lwrk, inxt, icond, ieval;
00040 extern logical lsame_(char *, char *);
00041 extern int dcopy_(integer *, doublereal *, integer *,
00042 doublereal *, integer *), dswap_(integer *, doublereal *, integer
00043 *, doublereal *, integer *);
00044 logical cursl;
00045 integer liwrk;
00046 extern int dlabad_(doublereal *, doublereal *), dgebak_(
00047 char *, char *, integer *, integer *, integer *, doublereal *,
00048 integer *, doublereal *, integer *, integer *),
00049 dgebal_(char *, integer *, doublereal *, integer *, integer *,
00050 integer *, doublereal *, integer *);
00051 logical lst2sl, scalea;
00052 extern doublereal dlamch_(char *);
00053 doublereal cscale;
00054 extern doublereal dlange_(char *, integer *, integer *, doublereal *,
00055 integer *, doublereal *);
00056 extern int dgehrd_(integer *, integer *, integer *,
00057 doublereal *, integer *, doublereal *, doublereal *, integer *,
00058 integer *), dlascl_(char *, integer *, integer *, doublereal *,
00059 doublereal *, integer *, integer *, doublereal *, integer *,
00060 integer *), dlacpy_(char *, integer *, integer *,
00061 doublereal *, integer *, doublereal *, integer *),
00062 xerbla_(char *, integer *);
00063 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
00064 integer *, integer *);
00065 doublereal bignum;
00066 extern int dorghr_(integer *, integer *, integer *,
00067 doublereal *, integer *, doublereal *, doublereal *, integer *,
00068 integer *), dhseqr_(char *, char *, integer *, integer *, integer
00069 *, doublereal *, integer *, doublereal *, doublereal *,
00070 doublereal *, integer *, doublereal *, integer *, integer *);
00071 logical wantsb;
00072 extern int dtrsen_(char *, char *, logical *, integer *,
00073 doublereal *, integer *, doublereal *, integer *, doublereal *,
00074 doublereal *, integer *, doublereal *, doublereal *, doublereal *,
00075 integer *, integer *, integer *, integer *);
00076 logical wantse, lastsl;
00077 integer minwrk, maxwrk;
00078 logical wantsn;
00079 doublereal smlnum;
00080 integer hswork;
00081 logical wantst, lquery, wantsv, wantvs;
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
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275 a_dim1 = *lda;
00276 a_offset = 1 + a_dim1;
00277 a -= a_offset;
00278 --wr;
00279 --wi;
00280 vs_dim1 = *ldvs;
00281 vs_offset = 1 + vs_dim1;
00282 vs -= vs_offset;
00283 --work;
00284 --iwork;
00285 --bwork;
00286
00287
00288 *info = 0;
00289 wantvs = lsame_(jobvs, "V");
00290 wantst = lsame_(sort, "S");
00291 wantsn = lsame_(sense, "N");
00292 wantse = lsame_(sense, "E");
00293 wantsv = lsame_(sense, "V");
00294 wantsb = lsame_(sense, "B");
00295 lquery = *lwork == -1 || *liwork == -1;
00296 if (! wantvs && ! lsame_(jobvs, "N")) {
00297 *info = -1;
00298 } else if (! wantst && ! lsame_(sort, "N")) {
00299 *info = -2;
00300 } else if (! (wantsn || wantse || wantsv || wantsb) || ! wantst && !
00301 wantsn) {
00302 *info = -4;
00303 } else if (*n < 0) {
00304 *info = -5;
00305 } else if (*lda < max(1,*n)) {
00306 *info = -7;
00307 } else if (*ldvs < 1 || wantvs && *ldvs < *n) {
00308 *info = -12;
00309 }
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325 if (*info == 0) {
00326 liwrk = 1;
00327 if (*n == 0) {
00328 minwrk = 1;
00329 lwrk = 1;
00330 } else {
00331 maxwrk = (*n << 1) + *n * ilaenv_(&c__1, "DGEHRD", " ", n, &c__1,
00332 n, &c__0);
00333 minwrk = *n * 3;
00334
00335 dhseqr_("S", jobvs, n, &c__1, n, &a[a_offset], lda, &wr[1], &wi[1]
00336 , &vs[vs_offset], ldvs, &work[1], &c_n1, &ieval);
00337 hswork = (integer) work[1];
00338
00339 if (! wantvs) {
00340
00341 i__1 = maxwrk, i__2 = *n + hswork;
00342 maxwrk = max(i__1,i__2);
00343 } else {
00344
00345 i__1 = maxwrk, i__2 = (*n << 1) + (*n - 1) * ilaenv_(&c__1,
00346 "DORGHR", " ", n, &c__1, n, &c_n1);
00347 maxwrk = max(i__1,i__2);
00348
00349 i__1 = maxwrk, i__2 = *n + hswork;
00350 maxwrk = max(i__1,i__2);
00351 }
00352 lwrk = maxwrk;
00353 if (! wantsn) {
00354
00355 i__1 = lwrk, i__2 = *n + *n * *n / 2;
00356 lwrk = max(i__1,i__2);
00357 }
00358 if (wantsv || wantsb) {
00359 liwrk = *n * *n / 4;
00360 }
00361 }
00362 iwork[1] = liwrk;
00363 work[1] = (doublereal) lwrk;
00364
00365 if (*lwork < minwrk && ! lquery) {
00366 *info = -16;
00367 } else if (*liwork < 1 && ! lquery) {
00368 *info = -18;
00369 }
00370 }
00371
00372 if (*info != 0) {
00373 i__1 = -(*info);
00374 xerbla_("DGEESX", &i__1);
00375 return 0;
00376 }
00377
00378
00379
00380 if (*n == 0) {
00381 *sdim = 0;
00382 return 0;
00383 }
00384
00385
00386
00387 eps = dlamch_("P");
00388 smlnum = dlamch_("S");
00389 bignum = 1. / smlnum;
00390 dlabad_(&smlnum, &bignum);
00391 smlnum = sqrt(smlnum) / eps;
00392 bignum = 1. / smlnum;
00393
00394
00395
00396 anrm = dlange_("M", n, n, &a[a_offset], lda, dum);
00397 scalea = FALSE_;
00398 if (anrm > 0. && anrm < smlnum) {
00399 scalea = TRUE_;
00400 cscale = smlnum;
00401 } else if (anrm > bignum) {
00402 scalea = TRUE_;
00403 cscale = bignum;
00404 }
00405 if (scalea) {
00406 dlascl_("G", &c__0, &c__0, &anrm, &cscale, n, n, &a[a_offset], lda, &
00407 ierr);
00408 }
00409
00410
00411
00412
00413 ibal = 1;
00414 dgebal_("P", n, &a[a_offset], lda, &ilo, &ihi, &work[ibal], &ierr);
00415
00416
00417
00418
00419 itau = *n + ibal;
00420 iwrk = *n + itau;
00421 i__1 = *lwork - iwrk + 1;
00422 dgehrd_(n, &ilo, &ihi, &a[a_offset], lda, &work[itau], &work[iwrk], &i__1,
00423 &ierr);
00424
00425 if (wantvs) {
00426
00427
00428
00429 dlacpy_("L", n, n, &a[a_offset], lda, &vs[vs_offset], ldvs)
00430 ;
00431
00432
00433
00434
00435 i__1 = *lwork - iwrk + 1;
00436 dorghr_(n, &ilo, &ihi, &vs[vs_offset], ldvs, &work[itau], &work[iwrk],
00437 &i__1, &ierr);
00438 }
00439
00440 *sdim = 0;
00441
00442
00443
00444
00445 iwrk = itau;
00446 i__1 = *lwork - iwrk + 1;
00447 dhseqr_("S", jobvs, n, &ilo, &ihi, &a[a_offset], lda, &wr[1], &wi[1], &vs[
00448 vs_offset], ldvs, &work[iwrk], &i__1, &ieval);
00449 if (ieval > 0) {
00450 *info = ieval;
00451 }
00452
00453
00454
00455 if (wantst && *info == 0) {
00456 if (scalea) {
00457 dlascl_("G", &c__0, &c__0, &cscale, &anrm, n, &c__1, &wr[1], n, &
00458 ierr);
00459 dlascl_("G", &c__0, &c__0, &cscale, &anrm, n, &c__1, &wi[1], n, &
00460 ierr);
00461 }
00462 i__1 = *n;
00463 for (i__ = 1; i__ <= i__1; ++i__) {
00464 bwork[i__] = (*select)(&wr[i__], &wi[i__]);
00465
00466 }
00467
00468
00469
00470
00471
00472
00473
00474
00475 i__1 = *lwork - iwrk + 1;
00476 dtrsen_(sense, jobvs, &bwork[1], n, &a[a_offset], lda, &vs[vs_offset],
00477 ldvs, &wr[1], &wi[1], sdim, rconde, rcondv, &work[iwrk], &
00478 i__1, &iwork[1], liwork, &icond);
00479 if (! wantsn) {
00480
00481 i__1 = maxwrk, i__2 = *n + (*sdim << 1) * (*n - *sdim);
00482 maxwrk = max(i__1,i__2);
00483 }
00484 if (icond == -15) {
00485
00486
00487
00488 *info = -16;
00489 } else if (icond == -17) {
00490
00491
00492
00493 *info = -18;
00494 } else if (icond > 0) {
00495
00496
00497
00498 *info = icond + *n;
00499 }
00500 }
00501
00502 if (wantvs) {
00503
00504
00505
00506
00507 dgebak_("P", "R", n, &ilo, &ihi, &work[ibal], n, &vs[vs_offset], ldvs,
00508 &ierr);
00509 }
00510
00511 if (scalea) {
00512
00513
00514
00515 dlascl_("H", &c__0, &c__0, &cscale, &anrm, n, n, &a[a_offset], lda, &
00516 ierr);
00517 i__1 = *lda + 1;
00518 dcopy_(n, &a[a_offset], &i__1, &wr[1], &c__1);
00519 if ((wantsv || wantsb) && *info == 0) {
00520 dum[0] = *rcondv;
00521 dlascl_("G", &c__0, &c__0, &cscale, &anrm, &c__1, &c__1, dum, &
00522 c__1, &ierr);
00523 *rcondv = dum[0];
00524 }
00525 if (cscale == smlnum) {
00526
00527
00528
00529
00530
00531 if (ieval > 0) {
00532 i1 = ieval + 1;
00533 i2 = ihi - 1;
00534 i__1 = ilo - 1;
00535 dlascl_("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &wi[
00536 1], n, &ierr);
00537 } else if (wantst) {
00538 i1 = 1;
00539 i2 = *n - 1;
00540 } else {
00541 i1 = ilo;
00542 i2 = ihi - 1;
00543 }
00544 inxt = i1 - 1;
00545 i__1 = i2;
00546 for (i__ = i1; i__ <= i__1; ++i__) {
00547 if (i__ < inxt) {
00548 goto L20;
00549 }
00550 if (wi[i__] == 0.) {
00551 inxt = i__ + 1;
00552 } else {
00553 if (a[i__ + 1 + i__ * a_dim1] == 0.) {
00554 wi[i__] = 0.;
00555 wi[i__ + 1] = 0.;
00556 } else if (a[i__ + 1 + i__ * a_dim1] != 0. && a[i__ + (
00557 i__ + 1) * a_dim1] == 0.) {
00558 wi[i__] = 0.;
00559 wi[i__ + 1] = 0.;
00560 if (i__ > 1) {
00561 i__2 = i__ - 1;
00562 dswap_(&i__2, &a[i__ * a_dim1 + 1], &c__1, &a[(
00563 i__ + 1) * a_dim1 + 1], &c__1);
00564 }
00565 if (*n > i__ + 1) {
00566 i__2 = *n - i__ - 1;
00567 dswap_(&i__2, &a[i__ + (i__ + 2) * a_dim1], lda, &
00568 a[i__ + 1 + (i__ + 2) * a_dim1], lda);
00569 }
00570 dswap_(n, &vs[i__ * vs_dim1 + 1], &c__1, &vs[(i__ + 1)
00571 * vs_dim1 + 1], &c__1);
00572 a[i__ + (i__ + 1) * a_dim1] = a[i__ + 1 + i__ *
00573 a_dim1];
00574 a[i__ + 1 + i__ * a_dim1] = 0.;
00575 }
00576 inxt = i__ + 2;
00577 }
00578 L20:
00579 ;
00580 }
00581 }
00582 i__1 = *n - ieval;
00583
00584 i__3 = *n - ieval;
00585 i__2 = max(i__3,1);
00586 dlascl_("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &wi[ieval +
00587 1], &i__2, &ierr);
00588 }
00589
00590 if (wantst && *info == 0) {
00591
00592
00593
00594 lastsl = TRUE_;
00595 lst2sl = TRUE_;
00596 *sdim = 0;
00597 ip = 0;
00598 i__1 = *n;
00599 for (i__ = 1; i__ <= i__1; ++i__) {
00600 cursl = (*select)(&wr[i__], &wi[i__]);
00601 if (wi[i__] == 0.) {
00602 if (cursl) {
00603 ++(*sdim);
00604 }
00605 ip = 0;
00606 if (cursl && ! lastsl) {
00607 *info = *n + 2;
00608 }
00609 } else {
00610 if (ip == 1) {
00611
00612
00613
00614 cursl = cursl || lastsl;
00615 lastsl = cursl;
00616 if (cursl) {
00617 *sdim += 2;
00618 }
00619 ip = -1;
00620 if (cursl && ! lst2sl) {
00621 *info = *n + 2;
00622 }
00623 } else {
00624
00625
00626
00627 ip = 1;
00628 }
00629 }
00630 lst2sl = lastsl;
00631 lastsl = cursl;
00632
00633 }
00634 }
00635
00636 work[1] = (doublereal) maxwrk;
00637 if (wantsv || wantsb) {
00638
00639 i__1 = 1, i__2 = *sdim * (*n - *sdim);
00640 iwork[1] = max(i__1,i__2);
00641 } else {
00642 iwork[1] = 1;
00643 }
00644
00645 return 0;
00646
00647
00648
00649 }