35 #ifndef TEMPLATE_LAPACK_TGEVC_HEADER
36 #define TEMPLATE_LAPACK_TGEVC_HEADER
256 integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1,
257 vr_offset, i__1, i__2, i__3, i__4, i__5;
258 Treal d__1, d__2, d__3, d__4, d__5, d__6;
261 Treal dmin__, temp, suma[4] , sumb[4]
263 Treal cim2a, cim2b, cre2a, cre2b, temp2, bdiag[2];
280 Treal acoefa, bcoefa, cimaga, cimagb;
283 Treal bcoefi, ascale, bscale, creala;
288 Treal salfar, safmin;
289 Treal xscale, bignum;
295 #define suma_ref(a_1,a_2) suma[(a_2)*2 + a_1 - 3]
296 #define sumb_ref(a_1,a_2) sumb[(a_2)*2 + a_1 - 3]
297 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
298 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
299 #define vl_ref(a_1,a_2) vl[(a_2)*vl_dim1 + a_1]
300 #define vr_ref(a_1,a_2) vr[(a_2)*vr_dim1 + a_1]
301 #define sum_ref(a_1,a_2) sum[(a_2)*2 + a_1 - 3]
306 a_offset = 1 + a_dim1 * 1;
309 b_offset = 1 + b_dim1 * 1;
312 vl_offset = 1 + vl_dim1 * 1;
315 vr_offset = 1 + vr_dim1 * 1;
359 }
else if (ihwmny < 0) {
380 for (j = 1; j <= i__1; ++j) {
386 if (
a_ref(j + 1, j) != 0.) {
391 if (select[j] || select[j + 1]) {
411 for (j = 1; j <= i__1; ++j) {
412 if (
a_ref(j + 1, j) != 0.) {
418 if (
a_ref(j + 2, j + 1) != 0.) {
430 }
else if ( ( compl_AAAA && *ldvl < *n ) || *ldvl < 1) {
432 }
else if ( ( compr && *ldvr < *n ) || *ldvr < 1) {
434 }
else if (*mm < im) {
456 small = safmin * *n / ulp;
458 bignum = 1. / (safmin * *n);
474 for (j = 2; j <= i__1; ++j) {
477 if (
a_ref(j, j - 1) == 0.) {
483 for (i__ = 1; i__ <= i__2; ++i__) {
489 work[*n + j] = temp2;
493 for (i__ = iend + 1; i__ <= i__2; ++i__) {
503 ascale = 1. /
maxMACRO(anorm,safmin);
504 bscale = 1. /
maxMACRO(bnorm,safmin);
515 for (je = 1; je <= i__1; ++je) {
528 if (
a_ref(je + 1, je) != 0.) {
536 ilcomp = select[je] || select[je + 1];
548 if ((d__1 =
a_ref(je, je),
absMACRO(d__1)) <= safmin && (d__2 =
555 for (jr = 1; jr <= i__2; ++jr) {
567 for (jr = 1; jr <= i__2; ++jr) {
568 work[(*n << 1) + jr] = 0.;
581 d__3 = (d__1 =
a_ref(je, je),
absMACRO(d__1)) * ascale, d__4 = (
585 salfar = temp *
a_ref(je, je) * ascale;
586 sbeta = temp *
b_ref(je, je) * bscale;
587 acoef = sbeta * ascale;
588 bcoefr = salfar * bscale;
609 d__1 = scale, d__2 = 1. / (safmin *
maxMACRO(d__3,d__4));
612 acoef = ascale * (scale * sbeta);
614 acoef = scale * acoef;
617 bcoefr = bscale * (scale * salfar);
619 bcoefr = scale * bcoefr;
627 work[(*n << 1) + je] = 1.;
633 d__1 = safmin * 100.;
635 acoef, &temp, &bcoefr, &temp2, &bcoefi);
647 if (acoefa * ulp < safmin && acoefa >= safmin) {
648 scale = safmin / ulp / acoefa;
650 if (bcoefa * ulp < safmin && bcoefa >= safmin) {
652 d__1 = scale, d__2 = safmin / ulp / bcoefa;
655 if (safmin * acoefa > ascale) {
656 scale = ascale / (safmin * acoefa);
658 if (safmin * bcoefa > bscale) {
660 d__1 = scale, d__2 = bscale / (safmin * bcoefa);
664 acoef = scale * acoef;
666 bcoefr = scale * bcoefr;
667 bcoefi = scale * bcoefi;
673 temp = acoef *
a_ref(je + 1, je);
674 temp2r = acoef * a_ref(je, je) - bcoefr *
b_ref(je, je);
675 temp2i = -bcoefi * b_ref(je, je);
677 work[(*n << 1) + je] = 1.;
678 work[*n * 3 + je] = 0.;
679 work[(*n << 1) + je + 1] = -temp2r / temp;
680 work[*n * 3 + je + 1] = -temp2i / temp;
682 work[(*n << 1) + je + 1] = 1.;
683 work[*n * 3 + je + 1] = 0.;
684 temp = acoef * a_ref(je, je + 1);
685 work[(*n << 1) + je] = (bcoefr * b_ref(je + 1, je + 1) -
686 acoef * a_ref(je + 1, je + 1)) / temp;
687 work[*n * 3 + je] = bcoefi * b_ref(je + 1, je + 1) / temp;
690 d__5 = (d__1 = work[(*n << 1) + je],
absMACRO(d__1)) + (d__2 =
691 work[*n * 3 + je],
absMACRO(d__2)), d__6 = (d__3 = work[(*
692 n << 1) + je + 1],
absMACRO(d__3)) + (d__4 = work[*n * 3 +
698 d__1 = ulp * acoefa * anorm, d__2 = ulp * bcoefa * bnorm, d__1 =
711 for (j = je + nw; j <= i__2; ++j) {
718 bdiag[0] =
b_ref(j, j);
720 if (
a_ref(j + 1, j) != 0.) {
722 bdiag[1] =
b_ref(j + 1, j + 1);
731 d__1 = work[j], d__2 = work[*n + j], d__1 =
maxMACRO(d__1,d__2),
732 d__2 = acoefa * work[j] + bcoefa * work[*n + j];
736 d__1 = temp, d__2 = work[j + 1], d__1 =
maxMACRO(d__1,d__2),
737 d__2 = work[*n + j + 1], d__1 =
maxMACRO(d__1,d__2),
738 d__2 = acoefa * work[j + 1] + bcoefa * work[*n +
742 if (temp > bignum * xscale) {
744 for (jw = 0; jw <= i__3; ++jw) {
746 for (jr = je; jr <= i__4; ++jr) {
747 work[(jw + 2) * *n + jr] = xscale * work[(jw + 2)
788 for (jw = 1; jw <= i__3; ++jw) {
803 for (ja = 1; ja <= i__4; ++ja) {
808 for (jr = je; jr <= i__5; ++jr) {
810 + ja - 1) * work[(jw + 1) * *n + jr];
812 + ja - 1) * work[(jw + 1) * *n + jr];
833 for (ja = 1; ja <= i__3; ++ja) {
838 sumb_ref(ja, 2) + bcoefi * sumb_ref(ja, 1);
851 bdiag, &bdiag[1], sum, &c__2, &bcoefr, &bcoefi, &
852 work[(*n << 1) + j], n, &scale, &temp, &iinfo);
855 for (jw = 0; jw <= i__3; ++jw) {
857 for (jr = je; jr <= i__4; ++jr) {
858 work[(jw + 2) * *n + jr] = scale * work[(jw + 2) *
877 for (jw = 0; jw <= i__2; ++jw) {
880 (jw + 2) * *n + je], &c__1, &c_b37, &work[(jw + 4)
898 for (j = ibeg; j <= i__2; ++j) {
907 for (j = ibeg; j <= i__2; ++j) {
919 for (jw = 0; jw <= i__2; ++jw) {
921 for (jr = ibeg; jr <= i__3; ++jr) {
929 ieig = ieig + nw - 1;
944 for (je = *n; je >= 1; --je) {
960 if (
a_ref(je, je - 1) != 0.) {
968 ilcomp = select[je] || select[je - 1];
980 if ((d__1 =
a_ref(je, je),
absMACRO(d__1)) <= safmin && (d__2 =
987 for (jr = 1; jr <= i__1; ++jr) {
999 for (jw = 0; jw <= i__1; ++jw) {
1001 for (jr = 1; jr <= i__2; ++jr) {
1002 work[(jw + 2) * *n + jr] = 0.;
1017 d__3 = (d__1 =
a_ref(je, je),
absMACRO(d__1)) * ascale, d__4 = (
1021 salfar = temp *
a_ref(je, je) * ascale;
1022 sbeta = temp *
b_ref(je, je) * bscale;
1023 acoef = sbeta * ascale;
1024 bcoefr = salfar * bscale;
1045 d__1 = scale, d__2 = 1. / (safmin *
maxMACRO(d__3,d__4));
1048 acoef = ascale * (scale * sbeta);
1050 acoef = scale * acoef;
1053 bcoefr = bscale * (scale * salfar);
1055 bcoefr = scale * bcoefr;
1063 work[(*n << 1) + je] = 1.;
1070 for (jr = 1; jr <= i__1; ++jr) {
1071 work[(*n << 1) + jr] = bcoefr *
b_ref(jr, je) - acoef *
1079 d__1 = safmin * 100.;
1081 ldb, &d__1, &acoef, &temp, &bcoefr, &temp2, &bcoefi);
1092 if (acoefa * ulp < safmin && acoefa >= safmin) {
1093 scale = safmin / ulp / acoefa;
1095 if (bcoefa * ulp < safmin && bcoefa >= safmin) {
1097 d__1 = scale, d__2 = safmin / ulp / bcoefa;
1100 if (safmin * acoefa > ascale) {
1101 scale = ascale / (safmin * acoefa);
1103 if (safmin * bcoefa > bscale) {
1105 d__1 = scale, d__2 = bscale / (safmin * bcoefa);
1109 acoef = scale * acoef;
1111 bcoefr = scale * bcoefr;
1112 bcoefi = scale * bcoefi;
1119 temp = acoef *
a_ref(je, je - 1);
1120 temp2r = acoef * a_ref(je, je) - bcoefr *
b_ref(je, je);
1121 temp2i = -bcoefi * b_ref(je, je);
1123 work[(*n << 1) + je] = 1.;
1124 work[*n * 3 + je] = 0.;
1125 work[(*n << 1) + je - 1] = -temp2r / temp;
1126 work[*n * 3 + je - 1] = -temp2i / temp;
1128 work[(*n << 1) + je - 1] = 1.;
1129 work[*n * 3 + je - 1] = 0.;
1130 temp = acoef * a_ref(je - 1, je);
1131 work[(*n << 1) + je] = (bcoefr * b_ref(je - 1, je - 1) -
1132 acoef * a_ref(je - 1, je - 1)) / temp;
1133 work[*n * 3 + je] = bcoefi * b_ref(je - 1, je - 1) / temp;
1137 d__5 = (d__1 = work[(*n << 1) + je],
absMACRO(d__1)) + (d__2 =
1138 work[*n * 3 + je],
absMACRO(d__2)), d__6 = (d__3 = work[(*
1139 n << 1) + je - 1],
absMACRO(d__3)) + (d__4 = work[*n * 3 +
1146 creala = acoef * work[(*n << 1) + je - 1];
1147 cimaga = acoef * work[*n * 3 + je - 1];
1148 crealb = bcoefr * work[(*n << 1) + je - 1] - bcoefi * work[*n
1150 cimagb = bcoefi * work[(*n << 1) + je - 1] + bcoefr * work[*n
1152 cre2a = acoef * work[(*n << 1) + je];
1153 cim2a = acoef * work[*n * 3 + je];
1154 cre2b = bcoefr * work[(*n << 1) + je] - bcoefi * work[*n * 3
1156 cim2b = bcoefi * work[(*n << 1) + je] + bcoefr * work[*n * 3
1159 for (jr = 1; jr <= i__1; ++jr) {
1160 work[(*n << 1) + jr] = -creala * a_ref(jr, je - 1) +
1161 crealb * b_ref(jr, je - 1) - cre2a * a_ref(jr, je)
1162 + cre2b * b_ref(jr, je);
1163 work[*n * 3 + jr] = -cimaga * a_ref(jr, je - 1) + cimagb *
1164 b_ref(jr, je - 1) - cim2a * a_ref(jr, je) +
1165 cim2b * b_ref(jr, je);
1171 d__1 = ulp * acoefa * anorm, d__2 = ulp * bcoefa * bnorm, d__1 =
1178 for (j = je - nw; j >= 1; --j) {
1183 if (! il2by2 && j > 1) {
1184 if (
a_ref(j, j - 1) != 0.) {
1189 bdiag[0] =
b_ref(j, j);
1192 bdiag[1] =
b_ref(j + 1, j + 1);
1200 lda, bdiag, &bdiag[1], &work[(*n << 1) + j], n, &
1201 bcoefr, &bcoefi, sum, &c__2, &scale, &temp, &iinfo);
1205 for (jw = 0; jw <= i__1; ++jw) {
1207 for (jr = 1; jr <= i__2; ++jr) {
1208 work[(jw + 2) * *n + jr] = scale * work[(jw + 2) *
1216 d__1 = scale * xmax;
1220 for (jw = 1; jw <= i__1; ++jw) {
1222 for (ja = 1; ja <= i__2; ++ja) {
1223 work[(jw + 1) * *n + j + ja - 1] =
sum_ref(ja, jw);
1236 temp = acoefa * work[j] + bcoefa * work[*n + j];
1239 d__1 = temp, d__2 = acoefa * work[j + 1] + bcoefa *
1246 if (temp > bignum * xscale) {
1249 for (jw = 0; jw <= i__1; ++jw) {
1251 for (jr = 1; jr <= i__2; ++jr) {
1252 work[(jw + 2) * *n + jr] = xscale * work[(jw
1267 for (ja = 1; ja <= i__1; ++ja) {
1269 creala = acoef * work[(*n << 1) + j + ja - 1];
1270 cimaga = acoef * work[*n * 3 + j + ja - 1];
1271 crealb = bcoefr * work[(*n << 1) + j + ja - 1] -
1272 bcoefi * work[*n * 3 + j + ja - 1];
1273 cimagb = bcoefi * work[(*n << 1) + j + ja - 1] +
1274 bcoefr * work[*n * 3 + j + ja - 1];
1276 for (jr = 1; jr <= i__2; ++jr) {
1277 work[(*n << 1) + jr] = work[(*n << 1) + jr] -
1278 creala *
a_ref(jr, j + ja - 1) +
1279 crealb *
b_ref(jr, j + ja - 1);
1280 work[*n * 3 + jr] = work[*n * 3 + jr] -
1281 cimaga *
a_ref(jr, j + ja - 1) +
1282 cimagb * b_ref(jr, j + ja - 1);
1286 creala = acoef * work[(*n << 1) + j + ja - 1];
1287 crealb = bcoefr * work[(*n << 1) + j + ja - 1];
1289 for (jr = 1; jr <= i__2; ++jr) {
1290 work[(*n << 1) + jr] = work[(*n << 1) + jr] -
1291 creala *
a_ref(jr, j + ja - 1) +
1292 crealb *
b_ref(jr, j + ja - 1);
1312 for (jw = 0; jw <= i__1; ++jw) {
1314 for (jr = 1; jr <= i__2; ++jr) {
1315 work[(jw + 4) * *n + jr] = work[(jw + 2) * *n + 1] *
1325 for (jc = 2; jc <= i__2; ++jc) {
1327 for (jr = 1; jr <= i__3; ++jr) {
1328 work[(jw + 4) * *n + jr] += work[(jw + 2) * *n +
1338 for (jw = 0; jw <= i__1; ++jw) {
1340 for (jr = 1; jr <= i__2; ++jr) {
1341 vr_ref(jr, ieig + jw) = work[(jw + 4) * *n + jr];
1350 for (jw = 0; jw <= i__1; ++jw) {
1352 for (jr = 1; jr <= i__2; ++jr) {
1353 vr_ref(jr, ieig + jw) = work[(jw + 2) * *n + jr];
1367 for (j = 1; j <= i__1; ++j) {
1376 for (j = 1; j <= i__1; ++j) {
1384 if (xmax > safmin) {
1387 for (jw = 0; jw <= i__1; ++jw) {
1389 for (jr = 1; jr <= i__2; ++jr) {