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__2 = 2;
00019 static integer c__1 = 1;
00020
00021 int ctgex2_(logical *wantq, logical *wantz, integer *n,
00022 complex *a, integer *lda, complex *b, integer *ldb, complex *q,
00023 integer *ldq, complex *z__, integer *ldz, integer *j1, integer *info)
00024 {
00025
00026 integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1,
00027 z_offset, i__1, i__2, i__3;
00028 real r__1;
00029 complex q__1, q__2, q__3;
00030
00031
00032 double sqrt(doublereal), c_abs(complex *);
00033 void r_cnjg(complex *, complex *);
00034
00035
00036 complex f, g;
00037 integer i__, m;
00038 complex s[4] , t[4] ;
00039 real cq, sa, sb, cz;
00040 complex sq;
00041 real ss, ws;
00042 complex sz;
00043 real eps, sum;
00044 logical weak;
00045 complex cdum;
00046 extern int crot_(integer *, complex *, integer *,
00047 complex *, integer *, real *, complex *);
00048 complex work[8];
00049 real scale;
00050 extern doublereal slamch_(char *);
00051 extern int clacpy_(char *, integer *, integer *, complex
00052 *, integer *, complex *, integer *), clartg_(complex *,
00053 complex *, real *, complex *, complex *), classq_(integer *,
00054 complex *, integer *, real *, real *);
00055 real thresh, smlnum;
00056 logical strong;
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 a_dim1 = *lda;
00183 a_offset = 1 + a_dim1;
00184 a -= a_offset;
00185 b_dim1 = *ldb;
00186 b_offset = 1 + b_dim1;
00187 b -= b_offset;
00188 q_dim1 = *ldq;
00189 q_offset = 1 + q_dim1;
00190 q -= q_offset;
00191 z_dim1 = *ldz;
00192 z_offset = 1 + z_dim1;
00193 z__ -= z_offset;
00194
00195
00196 *info = 0;
00197
00198
00199
00200 if (*n <= 1) {
00201 return 0;
00202 }
00203
00204 m = 2;
00205 weak = FALSE_;
00206 strong = FALSE_;
00207
00208
00209
00210 clacpy_("Full", &m, &m, &a[*j1 + *j1 * a_dim1], lda, s, &c__2);
00211 clacpy_("Full", &m, &m, &b[*j1 + *j1 * b_dim1], ldb, t, &c__2);
00212
00213
00214
00215 eps = slamch_("P");
00216 smlnum = slamch_("S") / eps;
00217 scale = 0.f;
00218 sum = 1.f;
00219 clacpy_("Full", &m, &m, s, &c__2, work, &m);
00220 clacpy_("Full", &m, &m, t, &c__2, &work[m * m], &m);
00221 i__1 = (m << 1) * m;
00222 classq_(&i__1, work, &c__1, &scale, &sum);
00223 sa = scale * sqrt(sum);
00224
00225 r__1 = eps * 10.f * sa;
00226 thresh = dmax(r__1,smlnum);
00227
00228
00229
00230
00231 q__2.r = s[3].r * t[0].r - s[3].i * t[0].i, q__2.i = s[3].r * t[0].i + s[
00232 3].i * t[0].r;
00233 q__3.r = t[3].r * s[0].r - t[3].i * s[0].i, q__3.i = t[3].r * s[0].i + t[
00234 3].i * s[0].r;
00235 q__1.r = q__2.r - q__3.r, q__1.i = q__2.i - q__3.i;
00236 f.r = q__1.r, f.i = q__1.i;
00237 q__2.r = s[3].r * t[2].r - s[3].i * t[2].i, q__2.i = s[3].r * t[2].i + s[
00238 3].i * t[2].r;
00239 q__3.r = t[3].r * s[2].r - t[3].i * s[2].i, q__3.i = t[3].r * s[2].i + t[
00240 3].i * s[2].r;
00241 q__1.r = q__2.r - q__3.r, q__1.i = q__2.i - q__3.i;
00242 g.r = q__1.r, g.i = q__1.i;
00243 sa = c_abs(&s[3]);
00244 sb = c_abs(&t[3]);
00245 clartg_(&g, &f, &cz, &sz, &cdum);
00246 q__1.r = -sz.r, q__1.i = -sz.i;
00247 sz.r = q__1.r, sz.i = q__1.i;
00248 r_cnjg(&q__1, &sz);
00249 crot_(&c__2, s, &c__1, &s[2], &c__1, &cz, &q__1);
00250 r_cnjg(&q__1, &sz);
00251 crot_(&c__2, t, &c__1, &t[2], &c__1, &cz, &q__1);
00252 if (sa >= sb) {
00253 clartg_(s, &s[1], &cq, &sq, &cdum);
00254 } else {
00255 clartg_(t, &t[1], &cq, &sq, &cdum);
00256 }
00257 crot_(&c__2, s, &c__2, &s[1], &c__2, &cq, &sq);
00258 crot_(&c__2, t, &c__2, &t[1], &c__2, &cq, &sq);
00259
00260
00261
00262 ws = c_abs(&s[1]) + c_abs(&t[1]);
00263 weak = ws <= thresh;
00264 if (! weak) {
00265 goto L20;
00266 }
00267
00268 if (TRUE_) {
00269
00270
00271
00272
00273 clacpy_("Full", &m, &m, s, &c__2, work, &m);
00274 clacpy_("Full", &m, &m, t, &c__2, &work[m * m], &m);
00275 r_cnjg(&q__2, &sz);
00276 q__1.r = -q__2.r, q__1.i = -q__2.i;
00277 crot_(&c__2, work, &c__1, &work[2], &c__1, &cz, &q__1);
00278 r_cnjg(&q__2, &sz);
00279 q__1.r = -q__2.r, q__1.i = -q__2.i;
00280 crot_(&c__2, &work[4], &c__1, &work[6], &c__1, &cz, &q__1);
00281 q__1.r = -sq.r, q__1.i = -sq.i;
00282 crot_(&c__2, work, &c__2, &work[1], &c__2, &cq, &q__1);
00283 q__1.r = -sq.r, q__1.i = -sq.i;
00284 crot_(&c__2, &work[4], &c__2, &work[5], &c__2, &cq, &q__1);
00285 for (i__ = 1; i__ <= 2; ++i__) {
00286 i__1 = i__ - 1;
00287 i__2 = i__ - 1;
00288 i__3 = *j1 + i__ - 1 + *j1 * a_dim1;
00289 q__1.r = work[i__2].r - a[i__3].r, q__1.i = work[i__2].i - a[i__3]
00290 .i;
00291 work[i__1].r = q__1.r, work[i__1].i = q__1.i;
00292 i__1 = i__ + 1;
00293 i__2 = i__ + 1;
00294 i__3 = *j1 + i__ - 1 + (*j1 + 1) * a_dim1;
00295 q__1.r = work[i__2].r - a[i__3].r, q__1.i = work[i__2].i - a[i__3]
00296 .i;
00297 work[i__1].r = q__1.r, work[i__1].i = q__1.i;
00298 i__1 = i__ + 3;
00299 i__2 = i__ + 3;
00300 i__3 = *j1 + i__ - 1 + *j1 * b_dim1;
00301 q__1.r = work[i__2].r - b[i__3].r, q__1.i = work[i__2].i - b[i__3]
00302 .i;
00303 work[i__1].r = q__1.r, work[i__1].i = q__1.i;
00304 i__1 = i__ + 5;
00305 i__2 = i__ + 5;
00306 i__3 = *j1 + i__ - 1 + (*j1 + 1) * b_dim1;
00307 q__1.r = work[i__2].r - b[i__3].r, q__1.i = work[i__2].i - b[i__3]
00308 .i;
00309 work[i__1].r = q__1.r, work[i__1].i = q__1.i;
00310
00311 }
00312 scale = 0.f;
00313 sum = 1.f;
00314 i__1 = (m << 1) * m;
00315 classq_(&i__1, work, &c__1, &scale, &sum);
00316 ss = scale * sqrt(sum);
00317 strong = ss <= thresh;
00318 if (! strong) {
00319 goto L20;
00320 }
00321 }
00322
00323
00324
00325
00326 i__1 = *j1 + 1;
00327 r_cnjg(&q__1, &sz);
00328 crot_(&i__1, &a[*j1 * a_dim1 + 1], &c__1, &a[(*j1 + 1) * a_dim1 + 1], &
00329 c__1, &cz, &q__1);
00330 i__1 = *j1 + 1;
00331 r_cnjg(&q__1, &sz);
00332 crot_(&i__1, &b[*j1 * b_dim1 + 1], &c__1, &b[(*j1 + 1) * b_dim1 + 1], &
00333 c__1, &cz, &q__1);
00334 i__1 = *n - *j1 + 1;
00335 crot_(&i__1, &a[*j1 + *j1 * a_dim1], lda, &a[*j1 + 1 + *j1 * a_dim1], lda,
00336 &cq, &sq);
00337 i__1 = *n - *j1 + 1;
00338 crot_(&i__1, &b[*j1 + *j1 * b_dim1], ldb, &b[*j1 + 1 + *j1 * b_dim1], ldb,
00339 &cq, &sq);
00340
00341
00342
00343 i__1 = *j1 + 1 + *j1 * a_dim1;
00344 a[i__1].r = 0.f, a[i__1].i = 0.f;
00345 i__1 = *j1 + 1 + *j1 * b_dim1;
00346 b[i__1].r = 0.f, b[i__1].i = 0.f;
00347
00348
00349
00350 if (*wantz) {
00351 r_cnjg(&q__1, &sz);
00352 crot_(n, &z__[*j1 * z_dim1 + 1], &c__1, &z__[(*j1 + 1) * z_dim1 + 1],
00353 &c__1, &cz, &q__1);
00354 }
00355 if (*wantq) {
00356 r_cnjg(&q__1, &sq);
00357 crot_(n, &q[*j1 * q_dim1 + 1], &c__1, &q[(*j1 + 1) * q_dim1 + 1], &
00358 c__1, &cq, &q__1);
00359 }
00360
00361
00362
00363 return 0;
00364
00365
00366
00367 L20:
00368 *info = 1;
00369 return 0;
00370
00371
00372
00373 }