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