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_GGEV_HEADER
00036 #define TEMPLATE_LAPACK_GGEV_HEADER
00037
00038
00039 template<class Treal>
00040 int template_lapack_ggev(const char *jobvl, const char *jobvr, const integer *n, Treal *
00041 a, const integer *lda, Treal *b, const integer *ldb, Treal *alphar,
00042 Treal *alphai, Treal *beta, Treal *vl, const integer *ldvl,
00043 Treal *vr, const integer *ldvr, Treal *work, const integer *lwork,
00044 integer *info)
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 integer c__1 = 1;
00182 integer c__0 = 0;
00183 Treal c_b26 = 0.;
00184 Treal c_b27 = 1.;
00185
00186
00187 integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1,
00188 vr_offset, i__1, i__2;
00189 Treal d__1, d__2, d__3, d__4;
00190
00191 Treal anrm, bnrm;
00192 integer ierr, itau;
00193 Treal temp;
00194 logical ilvl, ilvr;
00195 integer iwrk;
00196 integer ileft, icols, irows;
00197 integer jc;
00198 integer in;
00199 integer jr;
00200 logical ilascl, ilbscl;
00201 logical ldumma[1];
00202 char chtemp[1];
00203 Treal bignum;
00204 integer ijobvl, iright, ijobvr;
00205 Treal anrmto, bnrmto;
00206 integer minwrk, maxwrk;
00207 Treal smlnum;
00208 logical lquery;
00209 integer ihi, ilo;
00210 Treal eps;
00211 logical ilv;
00212 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00213 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
00214 #define vl_ref(a_1,a_2) vl[(a_2)*vl_dim1 + a_1]
00215 #define vr_ref(a_1,a_2) vr[(a_2)*vr_dim1 + a_1]
00216
00217
00218 a_dim1 = *lda;
00219 a_offset = 1 + a_dim1 * 1;
00220 a -= a_offset;
00221 b_dim1 = *ldb;
00222 b_offset = 1 + b_dim1 * 1;
00223 b -= b_offset;
00224 --alphar;
00225 --alphai;
00226 --beta;
00227 vl_dim1 = *ldvl;
00228 vl_offset = 1 + vl_dim1 * 1;
00229 vl -= vl_offset;
00230 vr_dim1 = *ldvr;
00231 vr_offset = 1 + vr_dim1 * 1;
00232 vr -= vr_offset;
00233 --work;
00234
00235
00236 maxwrk = 0;
00237
00238 if (template_blas_lsame(jobvl, "N")) {
00239 ijobvl = 1;
00240 ilvl = FALSE_;
00241 } else if (template_blas_lsame(jobvl, "V")) {
00242 ijobvl = 2;
00243 ilvl = TRUE_;
00244 } else {
00245 ijobvl = -1;
00246 ilvl = FALSE_;
00247 }
00248
00249 if (template_blas_lsame(jobvr, "N")) {
00250 ijobvr = 1;
00251 ilvr = FALSE_;
00252 } else if (template_blas_lsame(jobvr, "V")) {
00253 ijobvr = 2;
00254 ilvr = TRUE_;
00255 } else {
00256 ijobvr = -1;
00257 ilvr = FALSE_;
00258 }
00259 ilv = ilvl || ilvr;
00260
00261
00262
00263 *info = 0;
00264 lquery = *lwork == -1;
00265 if (ijobvl <= 0) {
00266 *info = -1;
00267 } else if (ijobvr <= 0) {
00268 *info = -2;
00269 } else if (*n < 0) {
00270 *info = -3;
00271 } else if (*lda < maxMACRO(1,*n)) {
00272 *info = -5;
00273 } else if (*ldb < maxMACRO(1,*n)) {
00274 *info = -7;
00275 } else if (*ldvl < 1 || ( ilvl && *ldvl < *n ) ) {
00276 *info = -12;
00277 } else if (*ldvr < 1 || ( ilvr && *ldvr < *n ) ) {
00278 *info = -14;
00279 }
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289 minwrk = 1;
00290 if (*info == 0 && (*lwork >= 1 || lquery)) {
00291 maxwrk = *n * 7 + *n * template_lapack_ilaenv(&c__1, "DGEQRF", " ", n, &c__1, n, &
00292 c__0, (ftnlen)6, (ftnlen)1);
00293
00294 i__1 = 1, i__2 = *n << 3;
00295 minwrk = maxMACRO(i__1,i__2);
00296 work[1] = (Treal) maxwrk;
00297 }
00298
00299 if (*lwork < minwrk && ! lquery) {
00300 *info = -16;
00301 }
00302
00303 if (*info != 0) {
00304 i__1 = -(*info);
00305 template_blas_erbla("GGEV ", &i__1);
00306 return 0;
00307 } else if (lquery) {
00308 return 0;
00309 }
00310
00311
00312
00313 if (*n == 0) {
00314 return 0;
00315 }
00316
00317
00318
00319 eps = template_lapack_lamch("P", (Treal)0);
00320 smlnum = template_lapack_lamch("S", (Treal)0);
00321 bignum = 1. / smlnum;
00322 template_lapack_labad(&smlnum, &bignum);
00323 smlnum = template_blas_sqrt(smlnum) / eps;
00324 bignum = 1. / smlnum;
00325
00326
00327
00328 anrm = template_lapack_lange("M", n, n, &a[a_offset], lda, &work[1]);
00329 ilascl = FALSE_;
00330 if (anrm > 0. && anrm < smlnum) {
00331 anrmto = smlnum;
00332 ilascl = TRUE_;
00333 } else if (anrm > bignum) {
00334 anrmto = bignum;
00335 ilascl = TRUE_;
00336 }
00337 if (ilascl) {
00338 template_lapack_lascl("G", &c__0, &c__0, &anrm, &anrmto, n, n, &a[a_offset], lda, &
00339 ierr);
00340 }
00341
00342
00343
00344 bnrm = template_lapack_lange("M", n, n, &b[b_offset], ldb, &work[1]);
00345 ilbscl = FALSE_;
00346 if (bnrm > 0. && bnrm < smlnum) {
00347 bnrmto = smlnum;
00348 ilbscl = TRUE_;
00349 } else if (bnrm > bignum) {
00350 bnrmto = bignum;
00351 ilbscl = TRUE_;
00352 }
00353 if (ilbscl) {
00354 template_lapack_lascl("G", &c__0, &c__0, &bnrm, &bnrmto, n, n, &b[b_offset], ldb, &
00355 ierr);
00356 }
00357
00358
00359
00360
00361 ileft = 1;
00362 iright = *n + 1;
00363 iwrk = iright + *n;
00364 template_lapack_ggbal("P", n, &a[a_offset], lda, &b[b_offset], ldb, &ilo, &ihi, &work[
00365 ileft], &work[iright], &work[iwrk], &ierr);
00366
00367
00368
00369
00370 irows = ihi + 1 - ilo;
00371 if (ilv) {
00372 icols = *n + 1 - ilo;
00373 } else {
00374 icols = irows;
00375 }
00376 itau = iwrk;
00377 iwrk = itau + irows;
00378 i__1 = *lwork + 1 - iwrk;
00379 template_lapack_geqrf(&irows, &icols, &b_ref(ilo, ilo), ldb, &work[itau], &work[iwrk], &
00380 i__1, &ierr);
00381
00382
00383
00384
00385 i__1 = *lwork + 1 - iwrk;
00386
00387 char str_L[] = {'L', 0};
00388 char str_T[] = {'T', 0};
00389 template_lapack_ormqr(str_L, str_T, &irows, &icols, &irows, &b_ref(ilo, ilo), ldb, &work[
00390 itau], &a_ref(ilo, ilo), lda, &work[iwrk], &i__1, &ierr);
00391
00392
00393
00394
00395 if (ilvl) {
00396 template_lapack_laset("Full", n, n, &c_b26, &c_b27, &vl[vl_offset], ldvl)
00397 ;
00398 i__1 = irows - 1;
00399 i__2 = irows - 1;
00400 template_lapack_lacpy("L", &i__1, &i__2, &b_ref(ilo + 1, ilo), ldb, &vl_ref(ilo + 1,
00401 ilo), ldvl);
00402 i__1 = *lwork + 1 - iwrk;
00403 template_lapack_orgqr(&irows, &irows, &irows, &vl_ref(ilo, ilo), ldvl, &work[itau],
00404 &work[iwrk], &i__1, &ierr);
00405 }
00406
00407
00408
00409 if (ilvr) {
00410 template_lapack_laset("Full", n, n, &c_b26, &c_b27, &vr[vr_offset], ldvr)
00411 ;
00412 }
00413
00414
00415
00416
00417 if (ilv) {
00418
00419
00420
00421 template_lapack_gghrd(jobvl, jobvr, n, &ilo, &ihi, &a[a_offset], lda, &b[b_offset],
00422 ldb, &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, &ierr);
00423 } else {
00424 template_lapack_gghrd("N", "N", &irows, &c__1, &irows, &a_ref(ilo, ilo), lda, &
00425 b_ref(ilo, ilo), ldb, &vl[vl_offset], ldvl, &vr[vr_offset],
00426 ldvr, &ierr);
00427 }
00428
00429
00430
00431
00432
00433 iwrk = itau;
00434 if (ilv) {
00435 *(unsigned char *)chtemp = 'S';
00436 } else {
00437 *(unsigned char *)chtemp = 'E';
00438 }
00439 i__1 = *lwork + 1 - iwrk;
00440 template_lapack_hgeqz(chtemp, jobvl, jobvr, n, &ilo, &ihi, &a[a_offset], lda, &b[
00441 b_offset], ldb, &alphar[1], &alphai[1], &beta[1], &vl[vl_offset],
00442 ldvl, &vr[vr_offset], ldvr, &work[iwrk], &i__1, &ierr);
00443 if (ierr != 0) {
00444 if (ierr > 0 && ierr <= *n) {
00445 *info = ierr;
00446 } else if (ierr > *n && ierr <= *n << 1) {
00447 *info = ierr - *n;
00448 } else {
00449 *info = *n + 1;
00450 }
00451 goto L110;
00452 }
00453
00454
00455
00456
00457 if (ilv) {
00458 if (ilvl) {
00459 if (ilvr) {
00460 *(unsigned char *)chtemp = 'B';
00461 } else {
00462 *(unsigned char *)chtemp = 'L';
00463 }
00464 } else {
00465 *(unsigned char *)chtemp = 'R';
00466 }
00467 template_lapack_tgevc(chtemp, "B", ldumma, n, &a[a_offset], lda, &b[b_offset], ldb,
00468 &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, n, &in, &work[
00469 iwrk], &ierr);
00470 if (ierr != 0) {
00471 *info = *n + 2;
00472 goto L110;
00473 }
00474
00475
00476
00477
00478 if (ilvl) {
00479 template_lapack_ggbak("P", "L", n, &ilo, &ihi, &work[ileft], &work[iright], n, &
00480 vl[vl_offset], ldvl, &ierr);
00481 i__1 = *n;
00482 for (jc = 1; jc <= i__1; ++jc) {
00483 if (alphai[jc] < 0.) {
00484 goto L50;
00485 }
00486 temp = 0.;
00487 if (alphai[jc] == 0.) {
00488 i__2 = *n;
00489 for (jr = 1; jr <= i__2; ++jr) {
00490
00491 d__2 = temp, d__3 = (d__1 = vl_ref(jr, jc), absMACRO(d__1))
00492 ;
00493 temp = maxMACRO(d__2,d__3);
00494
00495 }
00496 } else {
00497 i__2 = *n;
00498 for (jr = 1; jr <= i__2; ++jr) {
00499
00500 d__3 = temp, d__4 = (d__1 = vl_ref(jr, jc), absMACRO(d__1))
00501 + (d__2 = vl_ref(jr, jc + 1), absMACRO(d__2));
00502 temp = maxMACRO(d__3,d__4);
00503
00504 }
00505 }
00506 if (temp < smlnum) {
00507 goto L50;
00508 }
00509 temp = 1. / temp;
00510 if (alphai[jc] == 0.) {
00511 i__2 = *n;
00512 for (jr = 1; jr <= i__2; ++jr) {
00513 vl_ref(jr, jc) = vl_ref(jr, jc) * temp;
00514
00515 }
00516 } else {
00517 i__2 = *n;
00518 for (jr = 1; jr <= i__2; ++jr) {
00519 vl_ref(jr, jc) = vl_ref(jr, jc) * temp;
00520 vl_ref(jr, jc + 1) = vl_ref(jr, jc + 1) * temp;
00521
00522 }
00523 }
00524 L50:
00525 ;
00526 }
00527 }
00528 if (ilvr) {
00529 template_lapack_ggbak("P", "R", n, &ilo, &ihi, &work[ileft], &work[iright], n, &
00530 vr[vr_offset], ldvr, &ierr);
00531 i__1 = *n;
00532 for (jc = 1; jc <= i__1; ++jc) {
00533 if (alphai[jc] < 0.) {
00534 goto L100;
00535 }
00536 temp = 0.;
00537 if (alphai[jc] == 0.) {
00538 i__2 = *n;
00539 for (jr = 1; jr <= i__2; ++jr) {
00540
00541 d__2 = temp, d__3 = (d__1 = vr_ref(jr, jc), absMACRO(d__1))
00542 ;
00543 temp = maxMACRO(d__2,d__3);
00544
00545 }
00546 } else {
00547 i__2 = *n;
00548 for (jr = 1; jr <= i__2; ++jr) {
00549
00550 d__3 = temp, d__4 = (d__1 = vr_ref(jr, jc), absMACRO(d__1))
00551 + (d__2 = vr_ref(jr, jc + 1), absMACRO(d__2));
00552 temp = maxMACRO(d__3,d__4);
00553
00554 }
00555 }
00556 if (temp < smlnum) {
00557 goto L100;
00558 }
00559 temp = 1. / temp;
00560 if (alphai[jc] == 0.) {
00561 i__2 = *n;
00562 for (jr = 1; jr <= i__2; ++jr) {
00563 vr_ref(jr, jc) = vr_ref(jr, jc) * temp;
00564
00565 }
00566 } else {
00567 i__2 = *n;
00568 for (jr = 1; jr <= i__2; ++jr) {
00569 vr_ref(jr, jc) = vr_ref(jr, jc) * temp;
00570 vr_ref(jr, jc + 1) = vr_ref(jr, jc + 1) * temp;
00571
00572 }
00573 }
00574 L100:
00575 ;
00576 }
00577 }
00578
00579
00580
00581 }
00582
00583
00584
00585 if (ilascl) {
00586 template_lapack_lascl("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alphar[1], n, &
00587 ierr);
00588 template_lapack_lascl("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alphai[1], n, &
00589 ierr);
00590 }
00591
00592 if (ilbscl) {
00593 template_lapack_lascl("G", &c__0, &c__0, &bnrmto, &bnrm, n, &c__1, &beta[1], n, &
00594 ierr);
00595 }
00596
00597 L110:
00598
00599 work[1] = (Treal) maxwrk;
00600
00601 return 0;
00602
00603
00604
00605 }
00606
00607 #undef vr_ref
00608 #undef vl_ref
00609 #undef b_ref
00610 #undef a_ref
00611
00612
00613 #endif