35 #ifndef TEMPLATE_LAPACK_HGEQZ_HEADER
36 #define TEMPLATE_LAPACK_HGEQZ_HEADER
42 b,
const integer *ldb, Treal *alphar, Treal *alphai, Treal *
43 beta, Treal *q,
const integer *ldq, Treal *z__,
const integer *ldz,
254 integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1,
255 z_offset, i__1, i__2, i__3, i__4;
256 Treal d__1, d__2, d__3, d__4;
258 Treal ad11l, ad12l, ad21l, ad22l, ad32l, wabs, atol, btol,
260 Treal temp2, s1inv, c__;
262 Treal s, t, v[3], scale;
266 Treal tempi, tempr, s1, s2, u1, u2;
268 Treal a11, a12, a21, a22, b11, b22, c12, c21;
270 Treal an, bn, cl, cq,
cr;
272 Treal ascale, bscale, u12, w11;
274 Treal cz, sl, w12, w21, w22, wi;
295 Treal wr2, ad11, ad12, ad21, ad22, c11i, c22i;
297 Treal c11r, c22r, u12l;
301 Treal ulp, sqr, szi, szr;
302 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
303 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
304 #define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1]
305 #define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
309 a_offset = 1 + a_dim1 * 1;
312 b_offset = 1 + b_dim1 * 1;
318 q_offset = 1 + q_dim1 * 1;
321 z_offset = 1 + z_dim1 * 1;
326 ilschr = ilq = ilz = 0;
368 lquery = *lwork == -1;
371 }
else if (icompq == 0) {
373 }
else if (icompz == 0) {
377 }
else if (*ilo < 1) {
379 }
else if (*ihi > *n || *ihi < *ilo - 1) {
381 }
else if (*lda < *n) {
383 }
else if (*ldb < *n) {
385 }
else if (*ldq < 1 || ( ilq && *ldq < *n ) ) {
387 }
else if (*ldz < 1 || ( ilz && *ldz < *n ) ) {
389 }
else if (*lwork <
maxMACRO(1,*n) && ! lquery) {
418 in = *ihi + 1 - *ilo;
420 safmax = 1. / safmin;
422 anorm =
dlanhs_(
"F", &in, &
a_ref(*ilo, *ilo), lda, &work[1]);
423 bnorm =
dlanhs_(
"F", &in, &
b_ref(*ilo, *ilo), ldb, &work[1]);
425 d__1 = safmin, d__2 = ulp * anorm;
428 d__1 = safmin, d__2 = ulp * bnorm;
430 ascale = 1. /
maxMACRO(safmin,anorm);
431 bscale = 1. /
maxMACRO(safmin,bnorm);
436 for (j = *ihi + 1; j <= i__1; ++j) {
437 if (
b_ref(j, j) < 0.) {
440 for (jr = 1; jr <= i__2; ++jr) {
451 for (jr = 1; jr <= i__2; ++jr) {
457 alphar[j] =
a_ref(j, j);
459 beta[j] =
b_ref(j, j);
494 maxit = (*ihi - *ilo + 1) * 30;
497 for (jiter = 1; jiter <= i__1; ++jiter) {
511 if ((d__1 =
a_ref(ilast, ilast - 1),
absMACRO(d__1)) <= atol) {
512 a_ref(ilast, ilast - 1) = 0.;
518 b_ref(ilast, ilast) = 0.;
525 for (j = ilast - 1; j >= i__2; --j) {
533 a_ref(j, j - 1) = 0.;
552 if (tempr < 1. && tempr != 0.) {
556 if (temp * (ascale * (d__1 =
a_ref(j + 1, j),
absMACRO(d__1)))
557 <= temp2 * (ascale * atol)) {
568 if (ilazro || ilazr2) {
570 for (jch = j; jch <= i__3; ++jch) {
571 temp =
a_ref(jch, jch);
574 a_ref(jch + 1, jch) = 0.;
577 1, jch + 1), lda, &c__, &s);
580 1, jch + 1), ldb, &c__, &s);
586 a_ref(jch, jch - 1) =
a_ref(jch, jch - 1) * c__;
591 if (jch + 1 >= ilast) {
598 b_ref(jch + 1, jch + 1) = 0.;
608 for (jch = j; jch <= i__3; ++jch) {
609 temp =
b_ref(jch, jch + 1);
611 b_ref(jch, jch + 1));
612 b_ref(jch + 1, jch + 1) = 0.;
613 if (jch < ilastm - 1) {
614 i__4 = ilastm - jch - 1;
616 jch + 1, jch + 2), ldb, &c__, &s);
618 i__4 = ilastm - jch + 2;
620 1, jch - 1), lda, &c__, &s);
625 temp =
a_ref(jch + 1, jch);
627 a_ref(jch + 1, jch));
628 a_ref(jch + 1, jch - 1) = 0.;
629 i__4 = jch + 1 - ifrstm;
631 ifrstm, jch - 1), &c__1, &c__, &s);
634 ifrstm, jch - 1), &c__1, &c__, &s);
637 - 1), &c__1, &c__, &s);
665 temp =
a_ref(ilast, ilast);
668 a_ref(ilast, ilast - 1) = 0.;
669 i__2 = ilast - ifrstm;
672 i__2 = ilast - ifrstm;
684 if (
b_ref(ilast, ilast) < 0.) {
687 for (j = ifrstm; j <= i__2; ++j) {
698 for (j = 1; j <= i__2; ++j) {
704 alphar[ilast] =
a_ref(ilast, ilast);
706 beta[ilast] =
b_ref(ilast, ilast);
721 if (ifrstm > ilast) {
744 if (iiter / 10 * 10 == iiter) {
749 if ((Treal) maxit * safmin * (d__1 =
a_ref(ilast - 1, ilast),
752 eshift +=
a_ref(ilast - 1, ilast) /
b_ref(ilast - 1, ilast -
755 eshift += 1. / (safmin * (Treal) maxit);
766 d__1 = safmin * 100.;
768 - 1), ldb, &d__1, &s1, &s2, &wr, &wr2, &wi);
773 d__1 = s1, d__2 = safmin *
maxMACRO(d__3,d__4);
782 temp =
minMACRO(ascale,1.) * (safmax * .5);
789 temp =
minMACRO(bscale,1.) * (safmax * .5);
792 d__1 = scale, d__2 = temp /
absMACRO(wr);
801 for (j = ilast - 1; j >= i__2; --j) {
806 if (tempr < 1. && tempr != 0.) {
810 if ((d__1 = ascale *
a_ref(j + 1, j) * temp,
absMACRO(d__1)) <= ascale
824 temp = s1 *
a_ref(istart, istart) - wr *
b_ref(istart, istart);
825 temp2 = s1 *
a_ref(istart + 1, istart);
831 for (j = istart; j <= i__2; ++j) {
833 temp =
a_ref(j, j - 1);
836 a_ref(j + 1, j - 1) = 0.;
840 for (jc = j; jc <= i__3; ++jc) {
841 temp = c__ *
a_ref(j, jc) + s *
a_ref(j + 1, jc);
844 temp2 = c__ *
b_ref(j, jc) + s *
b_ref(j + 1, jc);
846 b_ref(j, jc) = temp2;
851 for (jr = 1; jr <= i__3; ++jr) {
852 temp = c__ *
q_ref(jr, j) + s *
q_ref(jr, j + 1);
860 temp =
b_ref(j + 1, j + 1);
862 b_ref(j + 1, j) = 0.;
867 for (jr = ifrstm; jr <= i__3; ++jr) {
868 temp = c__ *
a_ref(jr, j + 1) + s *
a_ref(jr, j);
870 a_ref(jr, j + 1) = temp;
874 for (jr = ifrstm; jr <= i__3; ++jr) {
875 temp = c__ *
b_ref(jr, j + 1) + s *
b_ref(jr, j);
877 b_ref(jr, j + 1) = temp;
882 for (jr = 1; jr <= i__3; ++jr) {
903 if (ifirst + 1 == ilast) {
914 b_ref(ilast, ilast), &b22, &b11, &sr, &cr, &sl, &cl);
923 i__2 = ilastm + 1 - ifirst;
925 ilast - 1), lda, &cl, &sl);
926 i__2 = ilast + 1 - ifrstm;
928 ilast), &c__1, &cr, &sr);
930 if (ilast < ilastm) {
931 i__2 = ilastm - ilast;
933 ilast + 1), lda, &cl, &sl);
935 if (ifrstm < ilast - 1) {
936 i__2 = ifirst - ifrstm;
938 ilast), &c__1, &cr, &sr);
950 b_ref(ilast - 1, ilast - 1) = b11;
951 b_ref(ilast - 1, ilast) = 0.;
952 b_ref(ilast, ilast - 1) = 0.;
953 b_ref(ilast, ilast) = b22;
959 for (j = ifrstm; j <= i__2; ++j) {
967 for (j = 1; j <= i__2; ++j) {
978 d__1 = safmin * 100.;
980 - 1), ldb, &d__1, &s1, &temp, &wr, &temp2, &wi);
992 a11 =
a_ref(ilast - 1, ilast - 1);
993 a21 =
a_ref(ilast, ilast - 1);
994 a12 =
a_ref(ilast - 1, ilast);
995 a22 =
a_ref(ilast, ilast);
1003 c11r = s1 * a11 - wr * b11;
1007 c22r = s1 * a22 - wr * b22;
1027 szr = -c21 * tempr / t;
1028 szi = c21 * tempi / t;
1041 if (s1 * an > wabs * bn) {
1046 a1r = cz * a11 + szr * a12;
1048 a2r = cz * a21 + szr * a22;
1058 sqr = tempr * a2r + tempi * a2i;
1059 sqi = tempi * a2r - tempr * a2i;
1069 tempr = sqr * szr - sqi * szi;
1070 tempi = sqr * szi + sqi * szr;
1071 b1r = cq * cz * b11 + tempr * b22;
1074 b2r = cq * cz * b22 + tempr * b11;
1080 beta[ilast - 1] = b1a;
1082 alphar[ilast - 1] = wr * b1a * s1inv;
1083 alphai[ilast - 1] = wi * b1a * s1inv;
1084 alphar[ilast] = wr * b2a * s1inv;
1085 alphai[ilast] = -(wi * b2a) * s1inv;
1100 if (ifrstm > ilast) {
1119 ad11 = ascale *
a_ref(ilast - 1, ilast - 1) / (bscale *
b_ref(
1120 ilast - 1, ilast - 1));
1121 ad21 = ascale *
a_ref(ilast, ilast - 1) / (bscale *
b_ref(ilast -
1123 ad12 = ascale *
a_ref(ilast - 1, ilast) / (bscale *
b_ref(ilast,
1125 ad22 = ascale *
a_ref(ilast, ilast) / (bscale *
b_ref(ilast,
1127 u12 =
b_ref(ilast - 1, ilast) /
b_ref(ilast, ilast);
1128 ad11l = ascale *
a_ref(ifirst, ifirst) / (bscale *
b_ref(ifirst,
1130 ad21l = ascale *
a_ref(ifirst + 1, ifirst) / (bscale *
b_ref(
1132 ad12l = ascale *
a_ref(ifirst, ifirst + 1) / (bscale *
b_ref(
1133 ifirst + 1, ifirst + 1));
1134 ad22l = ascale *
a_ref(ifirst + 1, ifirst + 1) / (bscale *
b_ref(
1135 ifirst + 1, ifirst + 1));
1136 ad32l = ascale *
a_ref(ifirst + 2, ifirst + 1) / (bscale *
b_ref(
1137 ifirst + 1, ifirst + 1));
1138 u12l =
b_ref(ifirst, ifirst + 1) /
b_ref(ifirst + 1, ifirst + 1);
1140 v[0] = (ad11 - ad11l) * (ad22 - ad11l) - ad12 * ad21 + ad21 * u12
1141 * ad11l + (ad12l - ad11l * u12l) * ad21l;
1142 v[1] = (ad22l - ad11l - ad21l * u12l - (ad11 - ad11l) - (ad22 -
1143 ad11l) + ad21 * u12) * ad21l;
1144 v[2] = ad32l * ad21l;
1154 for (j = istart; j <= i__2; ++j) {
1161 v[0] =
a_ref(j, j - 1);
1162 v[1] =
a_ref(j + 1, j - 1);
1163 v[2] =
a_ref(j + 2, j - 1);
1167 a_ref(j + 1, j - 1) = 0.;
1168 a_ref(j + 2, j - 1) = 0.;
1172 for (jc = j; jc <= i__3; ++jc) {
1173 temp = tau * (
a_ref(j, jc) + v[1] *
a_ref(j + 1, jc) + v[
1174 2] *
a_ref(j + 2, jc));
1176 a_ref(j + 1, jc) =
a_ref(j + 1, jc) - temp * v[1];
1177 a_ref(j + 2, jc) =
a_ref(j + 2, jc) - temp * v[2];
1178 temp2 = tau * (
b_ref(j, jc) + v[1] *
b_ref(j + 1, jc) + v[
1179 2] *
b_ref(j + 2, jc));
1181 b_ref(j + 1, jc) =
b_ref(j + 1, jc) - temp2 * v[1];
1182 b_ref(j + 2, jc) =
b_ref(j + 2, jc) - temp2 * v[2];
1187 for (jr = 1; jr <= i__3; ++jr) {
1188 temp = tau * (
q_ref(jr, j) + v[1] *
q_ref(jr, j + 1)
1189 + v[2] *
q_ref(jr, j + 2));
1191 q_ref(jr, j + 1) =
q_ref(jr, j + 1) - temp * v[1];
1192 q_ref(jr, j + 2) =
q_ref(jr, j + 2) - temp * v[2];
1203 d__3 = (d__1 =
b_ref(j + 1, j + 1),
absMACRO(d__1)), d__4 = (d__2 =
1207 d__3 = (d__1 =
b_ref(j + 2, j + 1),
absMACRO(d__1)), d__4 = (d__2 =
1210 if (
maxMACRO(temp,temp2) < safmin) {
1215 }
else if (temp >= temp2) {
1216 w11 =
b_ref(j + 1, j + 1);
1217 w21 =
b_ref(j + 2, j + 1);
1218 w12 =
b_ref(j + 1, j + 2);
1219 w22 =
b_ref(j + 2, j + 2);
1220 u1 =
b_ref(j + 1, j);
1221 u2 =
b_ref(j + 2, j);
1223 w21 =
b_ref(j + 1, j + 1);
1224 w11 =
b_ref(j + 2, j + 1);
1225 w22 =
b_ref(j + 1, j + 2);
1226 w12 =
b_ref(j + 2, j + 2);
1227 u2 =
b_ref(j + 1, j);
1228 u1 =
b_ref(j + 2, j);
1260 scale = (d__1 = w22 / u2,
absMACRO(d__1));
1264 d__2 = scale, d__3 = (d__1 = w11 / u1,
absMACRO(d__1));
1270 u2 = scale * u2 / w22;
1271 u1 = (scale * u1 - w12 * u2) / w11;
1289 tau = scale / t + 1.;
1290 vs = -1. / (scale + t);
1300 for (jr = ifrstm; jr <= i__3; ++jr) {
1301 temp = tau * (
a_ref(jr, j) + v[1] *
a_ref(jr, j + 1) + v[
1302 2] *
a_ref(jr, j + 2));
1304 a_ref(jr, j + 1) =
a_ref(jr, j + 1) - temp * v[1];
1305 a_ref(jr, j + 2) =
a_ref(jr, j + 2) - temp * v[2];
1309 for (jr = ifrstm; jr <= i__3; ++jr) {
1310 temp = tau * (
b_ref(jr, j) + v[1] *
b_ref(jr, j + 1) + v[
1311 2] *
b_ref(jr, j + 2));
1313 b_ref(jr, j + 1) =
b_ref(jr, j + 1) - temp * v[1];
1314 b_ref(jr, j + 2) =
b_ref(jr, j + 2) - temp * v[2];
1319 for (jr = 1; jr <= i__3; ++jr) {
1321 1) + v[2] *
z___ref(jr, j + 2));
1328 b_ref(j + 1, j) = 0.;
1329 b_ref(j + 2, j) = 0.;
1338 temp =
a_ref(j, j - 1);
1340 a_ref(j + 1, j - 1) = 0.;
1343 for (jc = j; jc <= i__2; ++jc) {
1344 temp = c__ *
a_ref(j, jc) + s *
a_ref(j + 1, jc);
1346 a_ref(j, jc) = temp;
1347 temp2 = c__ *
b_ref(j, jc) + s *
b_ref(j + 1, jc);
1349 b_ref(j, jc) = temp2;
1354 for (jr = 1; jr <= i__2; ++jr) {
1355 temp = c__ *
q_ref(jr, j) + s *
q_ref(jr, j + 1);
1358 q_ref(jr, j) = temp;
1365 temp =
b_ref(j + 1, j + 1);
1367 b_ref(j + 1, j) = 0.;
1370 for (jr = ifrstm; jr <= i__2; ++jr) {
1371 temp = c__ *
a_ref(jr, j + 1) + s *
a_ref(jr, j);
1373 a_ref(jr, j + 1) = temp;
1377 for (jr = ifrstm; jr <= i__2; ++jr) {
1378 temp = c__ *
b_ref(jr, j + 1) + s *
b_ref(jr, j);
1380 b_ref(jr, j + 1) = temp;
1385 for (jr = 1; jr <= i__2; ++jr) {
1420 for (j = 1; j <= i__1; ++j) {
1421 if (
b_ref(j, j) < 0.) {
1424 for (jr = 1; jr <= i__2; ++jr) {
1435 for (jr = 1; jr <= i__2; ++jr) {
1441 alphar[j] =
a_ref(j, j);
1443 beta[j] =
b_ref(j, j);
1454 work[1] = (Treal) (*n);
#define absMACRO(x)
Definition: template_blas_common.h:45
int template_lapack_laset(const char *uplo, const integer *m, const integer *n, const Treal *alpha, const Treal *beta, Treal *a, const integer *lda)
Definition: template_lapack_laset.h:40
Treal template_lapack_lapy3(Treal *x, Treal *y, Treal *z__)
Definition: template_lapack_lapy3.h:40
int integer
Definition: template_blas_common.h:38
Treal template_lapack_lapy2(Treal *x, Treal *y)
Definition: template_lapack_lapy2.h:40
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
int template_blas_rot(const integer *n, Treal *dx, const integer *incx, Treal *dy, const integer *incy, const Treal *c__, const Treal *s)
Definition: template_blas_rot.h:40
#define minMACRO(a, b)
Definition: template_blas_common.h:44
int template_lapack_lag2(const Treal *a, const integer *lda, const Treal *b, const integer *ldb, const Treal *safmin, Treal *scale1, Treal *scale2, Treal *wr1, Treal *wr2, Treal *wi)
Definition: template_lapack_lag2.h:40
#define z___ref(a_1, a_2)
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
int template_lapack_lartg(const Treal *f, const Treal *g, Treal *cs, Treal *sn, Treal *r__)
Definition: template_lapack_lartg.h:40
int template_lapack_hgeqz(const char *job, const char *compq, const char *compz, const integer *n, const integer *ilo, const integer *ihi, Treal *a, const integer *lda, Treal *b, const integer *ldb, Treal *alphar, Treal *alphai, Treal *beta, Treal *q, const integer *ldq, Treal *z__, const integer *ldz, Treal *work, const integer *lwork, integer *info)
Definition: template_lapack_hgeqz.h:40
int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x, const integer *incx, Treal *tau)
Definition: template_lapack_larfg.h:41
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:199
bool logical
Definition: template_blas_common.h:39
#define TRUE_
Definition: template_lapack_common.h:40
#define FALSE_
Definition: template_lapack_common.h:41
int template_lapack_lasv2(const Treal *f, const Treal *g, const Treal *h__, Treal *ssmin, Treal *ssmax, Treal *snr, Treal *csr, Treal *snl, Treal *csl)
Definition: template_lapack_lasv2.h:40
Treal template_blas_sqrt(Treal x)
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44
Treal dlanhs_(const char *norm, const integer *n, const Treal *a, const integer *lda, Treal *work)
Definition: template_lapack_lanhs.h:40