1 /* ./src_f77/zggbal.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11 static doublereal c_b35 = 10.;
12 static doublereal c_b71 = .5;
13
zggbal_(char * job,integer * n,doublecomplex * a,integer * lda,doublecomplex * b,integer * ldb,integer * ilo,integer * ihi,doublereal * lscale,doublereal * rscale,doublereal * work,integer * info,ftnlen job_len)14 /* Subroutine */ int zggbal_(char *job, integer *n, doublecomplex *a, integer
15 *lda, doublecomplex *b, integer *ldb, integer *ilo, integer *ihi,
16 doublereal *lscale, doublereal *rscale, doublereal *work, integer *
17 info, ftnlen job_len)
18 {
19 /* System generated locals */
20 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4;
21 doublereal d__1, d__2, d__3;
22
23 /* Builtin functions */
24 double d_lg10(doublereal *), d_imag(doublecomplex *), z_abs(doublecomplex
25 *), d_sign(doublereal *, doublereal *), pow_di(doublereal *,
26 integer *);
27
28 /* Local variables */
29 static integer i__, j, k, l, m;
30 static doublereal t;
31 static integer jc;
32 static doublereal ta, tb, tc;
33 static integer ir;
34 static doublereal ew;
35 static integer it, nr, ip1, jp1, lm1;
36 static doublereal cab, rab, ewc, cor, sum;
37 static integer nrp2, icab, lcab;
38 static doublereal beta, coef;
39 static integer irab, lrab;
40 static doublereal basl, cmax;
41 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
42 integer *);
43 static doublereal coef2, coef5, gamma, alpha;
44 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
45 integer *);
46 extern logical lsame_(char *, char *, ftnlen, ftnlen);
47 static doublereal sfmin, sfmax;
48 static integer iflow;
49 extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
50 integer *, doublereal *, integer *);
51 static integer kount;
52 extern /* Subroutine */ int zswap_(integer *, doublecomplex *, integer *,
53 doublecomplex *, integer *);
54 extern doublereal dlamch_(char *, ftnlen);
55 static doublereal pgamma;
56 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), zdscal_(
57 integer *, doublereal *, doublecomplex *, integer *);
58 static integer lsfmin;
59 extern integer izamax_(integer *, doublecomplex *, integer *);
60 static integer lsfmax;
61
62
63 /* -- LAPACK routine (version 3.0) -- */
64 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
65 /* Courant Institute, Argonne National Lab, and Rice University */
66 /* September 30, 1994 */
67
68 /* .. Scalar Arguments .. */
69 /* .. */
70 /* .. Array Arguments .. */
71 /* .. */
72
73 /* Purpose */
74 /* ======= */
75
76 /* ZGGBAL balances a pair of general complex matrices (A,B). This */
77 /* involves, first, permuting A and B by similarity transformations to */
78 /* isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N */
79 /* elements on the diagonal; and second, applying a diagonal similarity */
80 /* transformation to rows and columns ILO to IHI to make the rows */
81 /* and columns as close in norm as possible. Both steps are optional. */
82
83 /* Balancing may reduce the 1-norm of the matrices, and improve the */
84 /* accuracy of the computed eigenvalues and/or eigenvectors in the */
85 /* generalized eigenvalue problem A*x = lambda*B*x. */
86
87 /* Arguments */
88 /* ========= */
89
90 /* JOB (input) CHARACTER*1 */
91 /* Specifies the operations to be performed on A and B: */
92 /* = 'N': none: simply set ILO = 1, IHI = N, LSCALE(I) = 1.0 */
93 /* and RSCALE(I) = 1.0 for i=1,...,N; */
94 /* = 'P': permute only; */
95 /* = 'S': scale only; */
96 /* = 'B': both permute and scale. */
97
98 /* N (input) INTEGER */
99 /* The order of the matrices A and B. N >= 0. */
100
101 /* A (input/output) COMPLEX*16 array, dimension (LDA,N) */
102 /* On entry, the input matrix A. */
103 /* On exit, A is overwritten by the balanced matrix. */
104 /* If JOB = 'N', A is not referenced. */
105
106 /* LDA (input) INTEGER */
107 /* The leading dimension of the array A. LDA >= max(1,N). */
108
109 /* B (input/output) COMPLEX*16 array, dimension (LDB,N) */
110 /* On entry, the input matrix B. */
111 /* On exit, B is overwritten by the balanced matrix. */
112 /* If JOB = 'N', B is not referenced. */
113
114 /* LDB (input) INTEGER */
115 /* The leading dimension of the array B. LDB >= max(1,N). */
116
117 /* ILO (output) INTEGER */
118 /* IHI (output) INTEGER */
119 /* ILO and IHI are set to integers such that on exit */
120 /* A(i,j) = 0 and B(i,j) = 0 if i > j and */
121 /* j = 1,...,ILO-1 or i = IHI+1,...,N. */
122 /* If JOB = 'N' or 'S', ILO = 1 and IHI = N. */
123
124 /* LSCALE (output) DOUBLE PRECISION array, dimension (N) */
125 /* Details of the permutations and scaling factors applied */
126 /* to the left side of A and B. If P(j) is the index of the */
127 /* row interchanged with row j, and D(j) is the scaling factor */
128 /* applied to row j, then */
129 /* LSCALE(j) = P(j) for J = 1,...,ILO-1 */
130 /* = D(j) for J = ILO,...,IHI */
131 /* = P(j) for J = IHI+1,...,N. */
132 /* The order in which the interchanges are made is N to IHI+1, */
133 /* then 1 to ILO-1. */
134
135 /* RSCALE (output) DOUBLE PRECISION array, dimension (N) */
136 /* Details of the permutations and scaling factors applied */
137 /* to the right side of A and B. If P(j) is the index of the */
138 /* column interchanged with column j, and D(j) is the scaling */
139 /* factor applied to column j, then */
140 /* RSCALE(j) = P(j) for J = 1,...,ILO-1 */
141 /* = D(j) for J = ILO,...,IHI */
142 /* = P(j) for J = IHI+1,...,N. */
143 /* The order in which the interchanges are made is N to IHI+1, */
144 /* then 1 to ILO-1. */
145
146 /* WORK (workspace) DOUBLE PRECISION array, dimension (6*N) */
147
148 /* INFO (output) INTEGER */
149 /* = 0: successful exit */
150 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
151
152 /* Further Details */
153 /* =============== */
154
155 /* See R.C. WARD, Balancing the generalized eigenvalue problem, */
156 /* SIAM J. Sci. Stat. Comp. 2 (1981), 141-152. */
157
158 /* ===================================================================== */
159
160 /* .. Parameters .. */
161 /* .. */
162 /* .. Local Scalars .. */
163 /* .. */
164 /* .. External Functions .. */
165 /* .. */
166 /* .. External Subroutines .. */
167 /* .. */
168 /* .. Intrinsic Functions .. */
169 /* .. */
170 /* .. Statement Functions .. */
171 /* .. */
172 /* .. Statement Function definitions .. */
173 /* .. */
174 /* .. Executable Statements .. */
175
176 /* Test the input parameters */
177
178 /* Parameter adjustments */
179 a_dim1 = *lda;
180 a_offset = 1 + a_dim1;
181 a -= a_offset;
182 b_dim1 = *ldb;
183 b_offset = 1 + b_dim1;
184 b -= b_offset;
185 --lscale;
186 --rscale;
187 --work;
188
189 /* Function Body */
190 *info = 0;
191 if (! lsame_(job, "N", (ftnlen)1, (ftnlen)1) && ! lsame_(job, "P", (
192 ftnlen)1, (ftnlen)1) && ! lsame_(job, "S", (ftnlen)1, (ftnlen)1)
193 && ! lsame_(job, "B", (ftnlen)1, (ftnlen)1)) {
194 *info = -1;
195 } else if (*n < 0) {
196 *info = -2;
197 } else if (*lda < max(1,*n)) {
198 *info = -4;
199 } else if (*ldb < max(1,*n)) {
200 *info = -5;
201 }
202 if (*info != 0) {
203 i__1 = -(*info);
204 xerbla_("ZGGBAL", &i__1, (ftnlen)6);
205 return 0;
206 }
207
208 k = 1;
209 l = *n;
210
211 /* Quick return if possible */
212
213 if (*n == 0) {
214 return 0;
215 }
216
217 if (lsame_(job, "N", (ftnlen)1, (ftnlen)1)) {
218 *ilo = 1;
219 *ihi = *n;
220 i__1 = *n;
221 for (i__ = 1; i__ <= i__1; ++i__) {
222 lscale[i__] = 1.;
223 rscale[i__] = 1.;
224 /* L10: */
225 }
226 return 0;
227 }
228
229 if (k == l) {
230 *ilo = 1;
231 *ihi = 1;
232 lscale[1] = 1.;
233 rscale[1] = 1.;
234 return 0;
235 }
236
237 if (lsame_(job, "S", (ftnlen)1, (ftnlen)1)) {
238 goto L190;
239 }
240
241 goto L30;
242
243 /* Permute the matrices A and B to isolate the eigenvalues. */
244
245 /* Find row with one nonzero in columns 1 through L */
246
247 L20:
248 l = lm1;
249 if (l != 1) {
250 goto L30;
251 }
252
253 rscale[1] = 1.;
254 lscale[1] = 1.;
255 goto L190;
256
257 L30:
258 lm1 = l - 1;
259 for (i__ = l; i__ >= 1; --i__) {
260 i__1 = lm1;
261 for (j = 1; j <= i__1; ++j) {
262 jp1 = j + 1;
263 i__2 = i__ + j * a_dim1;
264 i__3 = i__ + j * b_dim1;
265 if (a[i__2].r != 0. || a[i__2].i != 0. || (b[i__3].r != 0. || b[
266 i__3].i != 0.)) {
267 goto L50;
268 }
269 /* L40: */
270 }
271 j = l;
272 goto L70;
273
274 L50:
275 i__1 = l;
276 for (j = jp1; j <= i__1; ++j) {
277 i__2 = i__ + j * a_dim1;
278 i__3 = i__ + j * b_dim1;
279 if (a[i__2].r != 0. || a[i__2].i != 0. || (b[i__3].r != 0. || b[
280 i__3].i != 0.)) {
281 goto L80;
282 }
283 /* L60: */
284 }
285 j = jp1 - 1;
286
287 L70:
288 m = l;
289 iflow = 1;
290 goto L160;
291 L80:
292 ;
293 }
294 goto L100;
295
296 /* Find column with one nonzero in rows K through N */
297
298 L90:
299 ++k;
300
301 L100:
302 i__1 = l;
303 for (j = k; j <= i__1; ++j) {
304 i__2 = lm1;
305 for (i__ = k; i__ <= i__2; ++i__) {
306 ip1 = i__ + 1;
307 i__3 = i__ + j * a_dim1;
308 i__4 = i__ + j * b_dim1;
309 if (a[i__3].r != 0. || a[i__3].i != 0. || (b[i__4].r != 0. || b[
310 i__4].i != 0.)) {
311 goto L120;
312 }
313 /* L110: */
314 }
315 i__ = l;
316 goto L140;
317 L120:
318 i__2 = l;
319 for (i__ = ip1; i__ <= i__2; ++i__) {
320 i__3 = i__ + j * a_dim1;
321 i__4 = i__ + j * b_dim1;
322 if (a[i__3].r != 0. || a[i__3].i != 0. || (b[i__4].r != 0. || b[
323 i__4].i != 0.)) {
324 goto L150;
325 }
326 /* L130: */
327 }
328 i__ = ip1 - 1;
329 L140:
330 m = k;
331 iflow = 2;
332 goto L160;
333 L150:
334 ;
335 }
336 goto L190;
337
338 /* Permute rows M and I */
339
340 L160:
341 lscale[m] = (doublereal) i__;
342 if (i__ == m) {
343 goto L170;
344 }
345 i__1 = *n - k + 1;
346 zswap_(&i__1, &a[i__ + k * a_dim1], lda, &a[m + k * a_dim1], lda);
347 i__1 = *n - k + 1;
348 zswap_(&i__1, &b[i__ + k * b_dim1], ldb, &b[m + k * b_dim1], ldb);
349
350 /* Permute columns M and J */
351
352 L170:
353 rscale[m] = (doublereal) j;
354 if (j == m) {
355 goto L180;
356 }
357 zswap_(&l, &a[j * a_dim1 + 1], &c__1, &a[m * a_dim1 + 1], &c__1);
358 zswap_(&l, &b[j * b_dim1 + 1], &c__1, &b[m * b_dim1 + 1], &c__1);
359
360 L180:
361 switch (iflow) {
362 case 1: goto L20;
363 case 2: goto L90;
364 }
365
366 L190:
367 *ilo = k;
368 *ihi = l;
369
370 if (*ilo == *ihi) {
371 return 0;
372 }
373
374 if (lsame_(job, "P", (ftnlen)1, (ftnlen)1)) {
375 return 0;
376 }
377
378 /* Balance the submatrix in rows ILO to IHI. */
379
380 nr = *ihi - *ilo + 1;
381 i__1 = *ihi;
382 for (i__ = *ilo; i__ <= i__1; ++i__) {
383 rscale[i__] = 0.;
384 lscale[i__] = 0.;
385
386 work[i__] = 0.;
387 work[i__ + *n] = 0.;
388 work[i__ + (*n << 1)] = 0.;
389 work[i__ + *n * 3] = 0.;
390 work[i__ + (*n << 2)] = 0.;
391 work[i__ + *n * 5] = 0.;
392 /* L200: */
393 }
394
395 /* Compute right side vector in resulting linear equations */
396
397 basl = d_lg10(&c_b35);
398 i__1 = *ihi;
399 for (i__ = *ilo; i__ <= i__1; ++i__) {
400 i__2 = *ihi;
401 for (j = *ilo; j <= i__2; ++j) {
402 i__3 = i__ + j * a_dim1;
403 if (a[i__3].r == 0. && a[i__3].i == 0.) {
404 ta = 0.;
405 goto L210;
406 }
407 i__3 = i__ + j * a_dim1;
408 d__3 = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[i__ + j *
409 a_dim1]), abs(d__2));
410 ta = d_lg10(&d__3) / basl;
411
412 L210:
413 i__3 = i__ + j * b_dim1;
414 if (b[i__3].r == 0. && b[i__3].i == 0.) {
415 tb = 0.;
416 goto L220;
417 }
418 i__3 = i__ + j * b_dim1;
419 d__3 = (d__1 = b[i__3].r, abs(d__1)) + (d__2 = d_imag(&b[i__ + j *
420 b_dim1]), abs(d__2));
421 tb = d_lg10(&d__3) / basl;
422
423 L220:
424 work[i__ + (*n << 2)] = work[i__ + (*n << 2)] - ta - tb;
425 work[j + *n * 5] = work[j + *n * 5] - ta - tb;
426 /* L230: */
427 }
428 /* L240: */
429 }
430
431 coef = 1. / (doublereal) (nr << 1);
432 coef2 = coef * coef;
433 coef5 = coef2 * .5;
434 nrp2 = nr + 2;
435 beta = 0.;
436 it = 1;
437
438 /* Start generalized conjugate gradient iteration */
439
440 L250:
441
442 gamma = ddot_(&nr, &work[*ilo + (*n << 2)], &c__1, &work[*ilo + (*n << 2)]
443 , &c__1) + ddot_(&nr, &work[*ilo + *n * 5], &c__1, &work[*ilo + *
444 n * 5], &c__1);
445
446 ew = 0.;
447 ewc = 0.;
448 i__1 = *ihi;
449 for (i__ = *ilo; i__ <= i__1; ++i__) {
450 ew += work[i__ + (*n << 2)];
451 ewc += work[i__ + *n * 5];
452 /* L260: */
453 }
454
455 /* Computing 2nd power */
456 d__1 = ew;
457 /* Computing 2nd power */
458 d__2 = ewc;
459 /* Computing 2nd power */
460 d__3 = ew - ewc;
461 gamma = coef * gamma - coef2 * (d__1 * d__1 + d__2 * d__2) - coef5 * (
462 d__3 * d__3);
463 if (gamma == 0.) {
464 goto L350;
465 }
466 if (it != 1) {
467 beta = gamma / pgamma;
468 }
469 t = coef5 * (ewc - ew * 3.);
470 tc = coef5 * (ew - ewc * 3.);
471
472 dscal_(&nr, &beta, &work[*ilo], &c__1);
473 dscal_(&nr, &beta, &work[*ilo + *n], &c__1);
474
475 daxpy_(&nr, &coef, &work[*ilo + (*n << 2)], &c__1, &work[*ilo + *n], &
476 c__1);
477 daxpy_(&nr, &coef, &work[*ilo + *n * 5], &c__1, &work[*ilo], &c__1);
478
479 i__1 = *ihi;
480 for (i__ = *ilo; i__ <= i__1; ++i__) {
481 work[i__] += tc;
482 work[i__ + *n] += t;
483 /* L270: */
484 }
485
486 /* Apply matrix to vector */
487
488 i__1 = *ihi;
489 for (i__ = *ilo; i__ <= i__1; ++i__) {
490 kount = 0;
491 sum = 0.;
492 i__2 = *ihi;
493 for (j = *ilo; j <= i__2; ++j) {
494 i__3 = i__ + j * a_dim1;
495 if (a[i__3].r == 0. && a[i__3].i == 0.) {
496 goto L280;
497 }
498 ++kount;
499 sum += work[j];
500 L280:
501 i__3 = i__ + j * b_dim1;
502 if (b[i__3].r == 0. && b[i__3].i == 0.) {
503 goto L290;
504 }
505 ++kount;
506 sum += work[j];
507 L290:
508 ;
509 }
510 work[i__ + (*n << 1)] = (doublereal) kount * work[i__ + *n] + sum;
511 /* L300: */
512 }
513
514 i__1 = *ihi;
515 for (j = *ilo; j <= i__1; ++j) {
516 kount = 0;
517 sum = 0.;
518 i__2 = *ihi;
519 for (i__ = *ilo; i__ <= i__2; ++i__) {
520 i__3 = i__ + j * a_dim1;
521 if (a[i__3].r == 0. && a[i__3].i == 0.) {
522 goto L310;
523 }
524 ++kount;
525 sum += work[i__ + *n];
526 L310:
527 i__3 = i__ + j * b_dim1;
528 if (b[i__3].r == 0. && b[i__3].i == 0.) {
529 goto L320;
530 }
531 ++kount;
532 sum += work[i__ + *n];
533 L320:
534 ;
535 }
536 work[j + *n * 3] = (doublereal) kount * work[j] + sum;
537 /* L330: */
538 }
539
540 sum = ddot_(&nr, &work[*ilo + *n], &c__1, &work[*ilo + (*n << 1)], &c__1)
541 + ddot_(&nr, &work[*ilo], &c__1, &work[*ilo + *n * 3], &c__1);
542 alpha = gamma / sum;
543
544 /* Determine correction to current iteration */
545
546 cmax = 0.;
547 i__1 = *ihi;
548 for (i__ = *ilo; i__ <= i__1; ++i__) {
549 cor = alpha * work[i__ + *n];
550 if (abs(cor) > cmax) {
551 cmax = abs(cor);
552 }
553 lscale[i__] += cor;
554 cor = alpha * work[i__];
555 if (abs(cor) > cmax) {
556 cmax = abs(cor);
557 }
558 rscale[i__] += cor;
559 /* L340: */
560 }
561 if (cmax < .5) {
562 goto L350;
563 }
564
565 d__1 = -alpha;
566 daxpy_(&nr, &d__1, &work[*ilo + (*n << 1)], &c__1, &work[*ilo + (*n << 2)]
567 , &c__1);
568 d__1 = -alpha;
569 daxpy_(&nr, &d__1, &work[*ilo + *n * 3], &c__1, &work[*ilo + *n * 5], &
570 c__1);
571
572 pgamma = gamma;
573 ++it;
574 if (it <= nrp2) {
575 goto L250;
576 }
577
578 /* End generalized conjugate gradient iteration */
579
580 L350:
581 sfmin = dlamch_("S", (ftnlen)1);
582 sfmax = 1. / sfmin;
583 lsfmin = (integer) (d_lg10(&sfmin) / basl + 1.);
584 lsfmax = (integer) (d_lg10(&sfmax) / basl);
585 i__1 = *ihi;
586 for (i__ = *ilo; i__ <= i__1; ++i__) {
587 i__2 = *n - *ilo + 1;
588 irab = izamax_(&i__2, &a[i__ + *ilo * a_dim1], lda);
589 rab = z_abs(&a[i__ + (irab + *ilo - 1) * a_dim1]);
590 i__2 = *n - *ilo + 1;
591 irab = izamax_(&i__2, &b[i__ + *ilo * b_dim1], lda);
592 /* Computing MAX */
593 d__1 = rab, d__2 = z_abs(&b[i__ + (irab + *ilo - 1) * b_dim1]);
594 rab = max(d__1,d__2);
595 d__1 = rab + sfmin;
596 lrab = (integer) (d_lg10(&d__1) / basl + 1.);
597 ir = (integer) (lscale[i__] + d_sign(&c_b71, &lscale[i__]));
598 /* Computing MIN */
599 i__2 = max(ir,lsfmin), i__2 = min(i__2,lsfmax), i__3 = lsfmax - lrab;
600 ir = min(i__2,i__3);
601 lscale[i__] = pow_di(&c_b35, &ir);
602 icab = izamax_(ihi, &a[i__ * a_dim1 + 1], &c__1);
603 cab = z_abs(&a[icab + i__ * a_dim1]);
604 icab = izamax_(ihi, &b[i__ * b_dim1 + 1], &c__1);
605 /* Computing MAX */
606 d__1 = cab, d__2 = z_abs(&b[icab + i__ * b_dim1]);
607 cab = max(d__1,d__2);
608 d__1 = cab + sfmin;
609 lcab = (integer) (d_lg10(&d__1) / basl + 1.);
610 jc = (integer) (rscale[i__] + d_sign(&c_b71, &rscale[i__]));
611 /* Computing MIN */
612 i__2 = max(jc,lsfmin), i__2 = min(i__2,lsfmax), i__3 = lsfmax - lcab;
613 jc = min(i__2,i__3);
614 rscale[i__] = pow_di(&c_b35, &jc);
615 /* L360: */
616 }
617
618 /* Row scaling of matrices A and B */
619
620 i__1 = *ihi;
621 for (i__ = *ilo; i__ <= i__1; ++i__) {
622 i__2 = *n - *ilo + 1;
623 zdscal_(&i__2, &lscale[i__], &a[i__ + *ilo * a_dim1], lda);
624 i__2 = *n - *ilo + 1;
625 zdscal_(&i__2, &lscale[i__], &b[i__ + *ilo * b_dim1], ldb);
626 /* L370: */
627 }
628
629 /* Column scaling of matrices A and B */
630
631 i__1 = *ihi;
632 for (j = *ilo; j <= i__1; ++j) {
633 zdscal_(ihi, &rscale[j], &a[j * a_dim1 + 1], &c__1);
634 zdscal_(ihi, &rscale[j], &b[j * b_dim1 + 1], &c__1);
635 /* L380: */
636 }
637
638 return 0;
639
640 /* End of ZGGBAL */
641
642 } /* zggbal_ */
643
644