1 /* lapack/complex16/zggbal.f -- translated by f2c (version 20090411).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #ifdef __cplusplus
14 extern "C" {
15 #endif
16 #include "v3p_netlib.h"
17
18 /* Table of constant values */
19
20 static integer c__1 = 1;
21 static doublereal c_b36 = 10.;
22 static doublereal c_b72 = .5;
23
24 /*< >*/
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)25 /* Subroutine */ int zggbal_(char *job, integer *n, doublecomplex *a, integer
26 *lda, doublecomplex *b, integer *ldb, integer *ilo, integer *ihi,
27 doublereal *lscale, doublereal *rscale, doublereal *work, integer *
28 info, ftnlen job_len)
29 {
30 /* System generated locals */
31 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4;
32 doublereal d__1, d__2, d__3;
33
34 /* Builtin functions */
35 double d_lg10(doublereal *), d_imag(doublecomplex *), z_abs(doublecomplex
36 *), d_sign(doublereal *, doublereal *), pow_di(doublereal *,
37 integer *);
38
39 /* Local variables */
40 integer i__, j, k, l, m;
41 doublereal t;
42 integer jc;
43 doublereal ta, tb, tc;
44 integer ir;
45 doublereal ew;
46 integer it, nr, ip1, jp1, lm1;
47 doublereal cab, rab, ewc, cor, sum;
48 integer nrp2, icab, lcab;
49 doublereal beta, coef;
50 integer irab, lrab;
51 doublereal basl, cmax;
52 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
53 integer *);
54 doublereal coef2, coef5, gamma, alpha;
55 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
56 integer *);
57 extern logical lsame_(const char *, const char *, ftnlen, ftnlen);
58 doublereal sfmin, sfmax;
59 integer iflow;
60 extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
61 integer *, doublereal *, integer *);
62 integer kount;
63 extern /* Subroutine */ int zswap_(integer *, doublecomplex *, integer *,
64 doublecomplex *, integer *);
65 extern doublereal dlamch_(char *, ftnlen);
66 doublereal pgamma;
67 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), zdscal_(
68 integer *, doublereal *, doublecomplex *, integer *);
69 integer lsfmin;
70 extern integer izamax_(integer *, doublecomplex *, integer *);
71 integer lsfmax;
72 (void)job_len;
73
74 /* -- LAPACK routine (version 3.2) -- */
75 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
76 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
77 /* November 2006 */
78
79 /* .. Scalar Arguments .. */
80 /*< CHARACTER JOB >*/
81 /*< INTEGER IHI, ILO, INFO, LDA, LDB, N >*/
82 /* .. */
83 /* .. Array Arguments .. */
84 /*< DOUBLE PRECISION LSCALE( * ), RSCALE( * ), WORK( * ) >*/
85 /*< COMPLEX*16 A( LDA, * ), B( LDB, * ) >*/
86 /* .. */
87
88 /* Purpose */
89 /* ======= */
90
91 /* ZGGBAL balances a pair of general complex matrices (A,B). This */
92 /* involves, first, permuting A and B by similarity transformations to */
93 /* isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N */
94 /* elements on the diagonal; and second, applying a diagonal similarity */
95 /* transformation to rows and columns ILO to IHI to make the rows */
96 /* and columns as close in norm as possible. Both steps are optional. */
97
98 /* Balancing may reduce the 1-norm of the matrices, and improve the */
99 /* accuracy of the computed eigenvalues and/or eigenvectors in the */
100 /* generalized eigenvalue problem A*x = lambda*B*x. */
101
102 /* Arguments */
103 /* ========= */
104
105 /* JOB (input) CHARACTER*1 */
106 /* Specifies the operations to be performed on A and B: */
107 /* = 'N': none: simply set ILO = 1, IHI = N, LSCALE(I) = 1.0 */
108 /* and RSCALE(I) = 1.0 for i=1,...,N; */
109 /* = 'P': permute only; */
110 /* = 'S': scale only; */
111 /* = 'B': both permute and scale. */
112
113 /* N (input) INTEGER */
114 /* The order of the matrices A and B. N >= 0. */
115
116 /* A (input/output) COMPLEX*16 array, dimension (LDA,N) */
117 /* On entry, the input matrix A. */
118 /* On exit, A is overwritten by the balanced matrix. */
119 /* If JOB = 'N', A is not referenced. */
120
121 /* LDA (input) INTEGER */
122 /* The leading dimension of the array A. LDA >= max(1,N). */
123
124 /* B (input/output) COMPLEX*16 array, dimension (LDB,N) */
125 /* On entry, the input matrix B. */
126 /* On exit, B is overwritten by the balanced matrix. */
127 /* If JOB = 'N', B is not referenced. */
128
129 /* LDB (input) INTEGER */
130 /* The leading dimension of the array B. LDB >= max(1,N). */
131
132 /* ILO (output) INTEGER */
133 /* IHI (output) INTEGER */
134 /* ILO and IHI are set to integers such that on exit */
135 /* A(i,j) = 0 and B(i,j) = 0 if i > j and */
136 /* j = 1,...,ILO-1 or i = IHI+1,...,N. */
137 /* If JOB = 'N' or 'S', ILO = 1 and IHI = N. */
138
139 /* LSCALE (output) DOUBLE PRECISION array, dimension (N) */
140 /* Details of the permutations and scaling factors applied */
141 /* to the left side of A and B. If P(j) is the index of the */
142 /* row interchanged with row j, and D(j) is the scaling factor */
143 /* applied to row j, then */
144 /* LSCALE(j) = P(j) for J = 1,...,ILO-1 */
145 /* = D(j) for J = ILO,...,IHI */
146 /* = P(j) for J = IHI+1,...,N. */
147 /* The order in which the interchanges are made is N to IHI+1, */
148 /* then 1 to ILO-1. */
149
150 /* RSCALE (output) DOUBLE PRECISION array, dimension (N) */
151 /* Details of the permutations and scaling factors applied */
152 /* to the right side of A and B. If P(j) is the index of the */
153 /* column interchanged with column j, and D(j) is the scaling */
154 /* factor applied to column j, then */
155 /* RSCALE(j) = P(j) for J = 1,...,ILO-1 */
156 /* = D(j) for J = ILO,...,IHI */
157 /* = P(j) for J = IHI+1,...,N. */
158 /* The order in which the interchanges are made is N to IHI+1, */
159 /* then 1 to ILO-1. */
160
161 /* WORK (workspace) REAL array, dimension (lwork) */
162 /* lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and */
163 /* at least 1 when JOB = 'N' or 'P'. */
164
165 /* INFO (output) INTEGER */
166 /* = 0: successful exit */
167 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
168
169 /* Further Details */
170 /* =============== */
171
172 /* See R.C. WARD, Balancing the generalized eigenvalue problem, */
173 /* SIAM J. Sci. Stat. Comp. 2 (1981), 141-152. */
174
175 /* ===================================================================== */
176
177 /* .. Parameters .. */
178 /*< DOUBLE PRECISION ZERO, HALF, ONE >*/
179 /*< PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) >*/
180 /*< DOUBLE PRECISION THREE, SCLFAC >*/
181 /*< PARAMETER ( THREE = 3.0D+0, SCLFAC = 1.0D+1 ) >*/
182 /*< COMPLEX*16 CZERO >*/
183 /*< PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) >*/
184 /* .. */
185 /* .. Local Scalars .. */
186 /*< >*/
187 /*< >*/
188 /*< COMPLEX*16 CDUM >*/
189 /* .. */
190 /* .. External Functions .. */
191 /*< LOGICAL LSAME >*/
192 /*< INTEGER IZAMAX >*/
193 /*< DOUBLE PRECISION DDOT, DLAMCH >*/
194 /*< EXTERNAL LSAME, IZAMAX, DDOT, DLAMCH >*/
195 /* .. */
196 /* .. External Subroutines .. */
197 /*< EXTERNAL DAXPY, DSCAL, XERBLA, ZDSCAL, ZSWAP >*/
198 /* .. */
199 /* .. Intrinsic Functions .. */
200 /*< INTRINSIC ABS, DBLE, DIMAG, INT, LOG10, MAX, MIN, SIGN >*/
201 /* .. */
202 /* .. Statement Functions .. */
203 /*< DOUBLE PRECISION CABS1 >*/
204 /* .. */
205 /* .. Statement Function definitions .. */
206 /*< CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) >*/
207 /* .. */
208 /* .. Executable Statements .. */
209
210 /* Test the input parameters */
211
212 /*< INFO = 0 >*/
213 /* Parameter adjustments */
214 a_dim1 = *lda;
215 a_offset = 1 + a_dim1;
216 a -= a_offset;
217 b_dim1 = *ldb;
218 b_offset = 1 + b_dim1;
219 b -= b_offset;
220 --lscale;
221 --rscale;
222 --work;
223
224 /* Function Body */
225 *info = 0;
226 /*< >*/
227 if (! lsame_(job, "N", (ftnlen)1, (ftnlen)1) && ! lsame_(job, "P", (
228 ftnlen)1, (ftnlen)1) && ! lsame_(job, "S", (ftnlen)1, (ftnlen)1)
229 && ! lsame_(job, "B", (ftnlen)1, (ftnlen)1)) {
230 /*< INFO = -1 >*/
231 *info = -1;
232 /*< ELSE IF( N.LT.0 ) THEN >*/
233 } else if (*n < 0) {
234 /*< INFO = -2 >*/
235 *info = -2;
236 /*< ELSE IF( LDA.LT.MAX( 1, N ) ) THEN >*/
237 } else if (*lda < max(1,*n)) {
238 /*< INFO = -4 >*/
239 *info = -4;
240 /*< ELSE IF( LDB.LT.MAX( 1, N ) ) THEN >*/
241 } else if (*ldb < max(1,*n)) {
242 /*< INFO = -6 >*/
243 *info = -6;
244 /*< END IF >*/
245 }
246 /*< IF( INFO.NE.0 ) THEN >*/
247 if (*info != 0) {
248 /*< CALL XERBLA( 'ZGGBAL', -INFO ) >*/
249 i__1 = -(*info);
250 xerbla_("ZGGBAL", &i__1, (ftnlen)6);
251 /*< RETURN >*/
252 return 0;
253 /*< END IF >*/
254 }
255
256 /* Quick return if possible */
257
258 /*< IF( N.EQ.0 ) THEN >*/
259 if (*n == 0) {
260 /*< ILO = 1 >*/
261 *ilo = 1;
262 /*< IHI = N >*/
263 *ihi = *n;
264 /*< RETURN >*/
265 return 0;
266 /*< END IF >*/
267 }
268
269 /*< IF( N.EQ.1 ) THEN >*/
270 if (*n == 1) {
271 /*< ILO = 1 >*/
272 *ilo = 1;
273 /*< IHI = N >*/
274 *ihi = *n;
275 /*< LSCALE( 1 ) = ONE >*/
276 lscale[1] = 1.;
277 /*< RSCALE( 1 ) = ONE >*/
278 rscale[1] = 1.;
279 /*< RETURN >*/
280 return 0;
281 /*< END IF >*/
282 }
283
284 /*< IF( LSAME( JOB, 'N' ) ) THEN >*/
285 if (lsame_(job, "N", (ftnlen)1, (ftnlen)1)) {
286 /*< ILO = 1 >*/
287 *ilo = 1;
288 /*< IHI = N >*/
289 *ihi = *n;
290 /*< DO 10 I = 1, N >*/
291 i__1 = *n;
292 for (i__ = 1; i__ <= i__1; ++i__) {
293 /*< LSCALE( I ) = ONE >*/
294 lscale[i__] = 1.;
295 /*< RSCALE( I ) = ONE >*/
296 rscale[i__] = 1.;
297 /*< 10 CONTINUE >*/
298 /* L10: */
299 }
300 /*< RETURN >*/
301 return 0;
302 /*< END IF >*/
303 }
304
305 /*< K = 1 >*/
306 k = 1;
307 /*< L = N >*/
308 l = *n;
309 /*< >*/
310 if (lsame_(job, "S", (ftnlen)1, (ftnlen)1)) {
311 goto L190;
312 }
313
314 /*< GO TO 30 >*/
315 goto L30;
316
317 /* Permute the matrices A and B to isolate the eigenvalues. */
318
319 /* Find row with one nonzero in columns 1 through L */
320
321 /*< 20 CONTINUE >*/
322 L20:
323 /*< L = LM1 >*/
324 l = lm1;
325 /*< >*/
326 if (l != 1) {
327 goto L30;
328 }
329
330 /*< RSCALE( 1 ) = 1 >*/
331 rscale[1] = 1.;
332 /*< LSCALE( 1 ) = 1 >*/
333 lscale[1] = 1.;
334 /*< GO TO 190 >*/
335 goto L190;
336
337 /*< 30 CONTINUE >*/
338 L30:
339 /*< LM1 = L - 1 >*/
340 lm1 = l - 1;
341 /*< DO 80 I = L, 1, -1 >*/
342 for (i__ = l; i__ >= 1; --i__) {
343 /*< DO 40 J = 1, LM1 >*/
344 i__1 = lm1;
345 for (j = 1; j <= i__1; ++j) {
346 /*< JP1 = J + 1 >*/
347 jp1 = j + 1;
348 /*< >*/
349 i__2 = i__ + j * a_dim1;
350 i__3 = i__ + j * b_dim1;
351 if (a[i__2].r != 0. || a[i__2].i != 0. || (b[i__3].r != 0. || b[
352 i__3].i != 0.)) {
353 goto L50;
354 }
355 /*< 40 CONTINUE >*/
356 /* L40: */
357 }
358 /*< J = L >*/
359 j = l;
360 /*< GO TO 70 >*/
361 goto L70;
362
363 /*< 50 CONTINUE >*/
364 L50:
365 /*< DO 60 J = JP1, L >*/
366 i__1 = l;
367 for (j = jp1; j <= i__1; ++j) {
368 /*< >*/
369 i__2 = i__ + j * a_dim1;
370 i__3 = i__ + j * b_dim1;
371 if (a[i__2].r != 0. || a[i__2].i != 0. || (b[i__3].r != 0. || b[
372 i__3].i != 0.)) {
373 goto L80;
374 }
375 /*< 60 CONTINUE >*/
376 /* L60: */
377 }
378 /*< J = JP1 - 1 >*/
379 j = jp1 - 1;
380
381 /*< 70 CONTINUE >*/
382 L70:
383 /*< M = L >*/
384 m = l;
385 /*< IFLOW = 1 >*/
386 iflow = 1;
387 /*< GO TO 160 >*/
388 goto L160;
389 /*< 80 CONTINUE >*/
390 L80:
391 ;
392 }
393 /*< GO TO 100 >*/
394 goto L100;
395
396 /* Find column with one nonzero in rows K through N */
397
398 /*< 90 CONTINUE >*/
399 L90:
400 /*< K = K + 1 >*/
401 ++k;
402
403 /*< 100 CONTINUE >*/
404 L100:
405 /*< DO 150 J = K, L >*/
406 i__1 = l;
407 for (j = k; j <= i__1; ++j) {
408 /*< DO 110 I = K, LM1 >*/
409 i__2 = lm1;
410 for (i__ = k; i__ <= i__2; ++i__) {
411 /*< IP1 = I + 1 >*/
412 ip1 = i__ + 1;
413 /*< >*/
414 i__3 = i__ + j * a_dim1;
415 i__4 = i__ + j * b_dim1;
416 if (a[i__3].r != 0. || a[i__3].i != 0. || (b[i__4].r != 0. || b[
417 i__4].i != 0.)) {
418 goto L120;
419 }
420 /*< 110 CONTINUE >*/
421 /* L110: */
422 }
423 /*< I = L >*/
424 i__ = l;
425 /*< GO TO 140 >*/
426 goto L140;
427 /*< 120 CONTINUE >*/
428 L120:
429 /*< DO 130 I = IP1, L >*/
430 i__2 = l;
431 for (i__ = ip1; i__ <= i__2; ++i__) {
432 /*< >*/
433 i__3 = i__ + j * a_dim1;
434 i__4 = i__ + j * b_dim1;
435 if (a[i__3].r != 0. || a[i__3].i != 0. || (b[i__4].r != 0. || b[
436 i__4].i != 0.)) {
437 goto L150;
438 }
439 /*< 130 CONTINUE >*/
440 /* L130: */
441 }
442 /*< I = IP1 - 1 >*/
443 i__ = ip1 - 1;
444 /*< 140 CONTINUE >*/
445 L140:
446 /*< M = K >*/
447 m = k;
448 /*< IFLOW = 2 >*/
449 iflow = 2;
450 /*< GO TO 160 >*/
451 goto L160;
452 /*< 150 CONTINUE >*/
453 L150:
454 ;
455 }
456 /*< GO TO 190 >*/
457 goto L190;
458
459 /* Permute rows M and I */
460
461 /*< 160 CONTINUE >*/
462 L160:
463 /*< LSCALE( M ) = I >*/
464 lscale[m] = (doublereal) i__;
465 /*< >*/
466 if (i__ == m) {
467 goto L170;
468 }
469 /*< CALL ZSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA ) >*/
470 i__1 = *n - k + 1;
471 zswap_(&i__1, &a[i__ + k * a_dim1], lda, &a[m + k * a_dim1], lda);
472 /*< CALL ZSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB ) >*/
473 i__1 = *n - k + 1;
474 zswap_(&i__1, &b[i__ + k * b_dim1], ldb, &b[m + k * b_dim1], ldb);
475
476 /* Permute columns M and J */
477
478 /*< 170 CONTINUE >*/
479 L170:
480 /*< RSCALE( M ) = J >*/
481 rscale[m] = (doublereal) j;
482 /*< >*/
483 if (j == m) {
484 goto L180;
485 }
486 /*< CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 ) >*/
487 zswap_(&l, &a[j * a_dim1 + 1], &c__1, &a[m * a_dim1 + 1], &c__1);
488 /*< CALL ZSWAP( L, B( 1, J ), 1, B( 1, M ), 1 ) >*/
489 zswap_(&l, &b[j * b_dim1 + 1], &c__1, &b[m * b_dim1 + 1], &c__1);
490
491 /*< 180 CONTINUE >*/
492 L180:
493 /*< GO TO ( 20, 90 )IFLOW >*/
494 switch (iflow) {
495 case 1: goto L20;
496 case 2: goto L90;
497 }
498
499 /*< 190 CONTINUE >*/
500 L190:
501 /*< ILO = K >*/
502 *ilo = k;
503 /*< IHI = L >*/
504 *ihi = l;
505
506 /*< IF( LSAME( JOB, 'P' ) ) THEN >*/
507 if (lsame_(job, "P", (ftnlen)1, (ftnlen)1)) {
508 /*< DO 195 I = ILO, IHI >*/
509 i__1 = *ihi;
510 for (i__ = *ilo; i__ <= i__1; ++i__) {
511 /*< LSCALE( I ) = ONE >*/
512 lscale[i__] = 1.;
513 /*< RSCALE( I ) = ONE >*/
514 rscale[i__] = 1.;
515 /*< 195 CONTINUE >*/
516 /* L195: */
517 }
518 /*< RETURN >*/
519 return 0;
520 /*< END IF >*/
521 }
522
523 /*< >*/
524 if (*ilo == *ihi) {
525 return 0;
526 }
527
528 /* Balance the submatrix in rows ILO to IHI. */
529
530 /*< NR = IHI - ILO + 1 >*/
531 nr = *ihi - *ilo + 1;
532 /*< DO 200 I = ILO, IHI >*/
533 i__1 = *ihi;
534 for (i__ = *ilo; i__ <= i__1; ++i__) {
535 /*< RSCALE( I ) = ZERO >*/
536 rscale[i__] = 0.;
537 /*< LSCALE( I ) = ZERO >*/
538 lscale[i__] = 0.;
539
540 /*< WORK( I ) = ZERO >*/
541 work[i__] = 0.;
542 /*< WORK( I+N ) = ZERO >*/
543 work[i__ + *n] = 0.;
544 /*< WORK( I+2*N ) = ZERO >*/
545 work[i__ + (*n << 1)] = 0.;
546 /*< WORK( I+3*N ) = ZERO >*/
547 work[i__ + *n * 3] = 0.;
548 /*< WORK( I+4*N ) = ZERO >*/
549 work[i__ + (*n << 2)] = 0.;
550 /*< WORK( I+5*N ) = ZERO >*/
551 work[i__ + *n * 5] = 0.;
552 /*< 200 CONTINUE >*/
553 /* L200: */
554 }
555
556 /* Compute right side vector in resulting linear equations */
557
558 /*< BASL = LOG10( SCLFAC ) >*/
559 basl = d_lg10(&c_b36);
560 /*< DO 240 I = ILO, IHI >*/
561 i__1 = *ihi;
562 for (i__ = *ilo; i__ <= i__1; ++i__) {
563 /*< DO 230 J = ILO, IHI >*/
564 i__2 = *ihi;
565 for (j = *ilo; j <= i__2; ++j) {
566 /*< IF( A( I, J ).EQ.CZERO ) THEN >*/
567 i__3 = i__ + j * a_dim1;
568 if (a[i__3].r == 0. && a[i__3].i == 0.) {
569 /*< TA = ZERO >*/
570 ta = 0.;
571 /*< GO TO 210 >*/
572 goto L210;
573 /*< END IF >*/
574 }
575 /*< TA = LOG10( CABS1( A( I, J ) ) ) / BASL >*/
576 i__3 = i__ + j * a_dim1;
577 d__3 = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[i__ + j *
578 a_dim1]), abs(d__2));
579 ta = d_lg10(&d__3) / basl;
580
581 /*< 210 CONTINUE >*/
582 L210:
583 /*< IF( B( I, J ).EQ.CZERO ) THEN >*/
584 i__3 = i__ + j * b_dim1;
585 if (b[i__3].r == 0. && b[i__3].i == 0.) {
586 /*< TB = ZERO >*/
587 tb = 0.;
588 /*< GO TO 220 >*/
589 goto L220;
590 /*< END IF >*/
591 }
592 /*< TB = LOG10( CABS1( B( I, J ) ) ) / BASL >*/
593 i__3 = i__ + j * b_dim1;
594 d__3 = (d__1 = b[i__3].r, abs(d__1)) + (d__2 = d_imag(&b[i__ + j *
595 b_dim1]), abs(d__2));
596 tb = d_lg10(&d__3) / basl;
597
598 /*< 220 CONTINUE >*/
599 L220:
600 /*< WORK( I+4*N ) = WORK( I+4*N ) - TA - TB >*/
601 work[i__ + (*n << 2)] = work[i__ + (*n << 2)] - ta - tb;
602 /*< WORK( J+5*N ) = WORK( J+5*N ) - TA - TB >*/
603 work[j + *n * 5] = work[j + *n * 5] - ta - tb;
604 /*< 230 CONTINUE >*/
605 /* L230: */
606 }
607 /*< 240 CONTINUE >*/
608 /* L240: */
609 }
610
611 /*< COEF = ONE / DBLE( 2*NR ) >*/
612 coef = 1. / (doublereal) (nr << 1);
613 /*< COEF2 = COEF*COEF >*/
614 coef2 = coef * coef;
615 /*< COEF5 = HALF*COEF2 >*/
616 coef5 = coef2 * .5;
617 /*< NRP2 = NR + 2 >*/
618 nrp2 = nr + 2;
619 /*< BETA = ZERO >*/
620 beta = 0.;
621 /*< IT = 1 >*/
622 it = 1;
623
624 /* Start generalized conjugate gradient iteration */
625
626 /*< 250 CONTINUE >*/
627 L250:
628
629 /*< >*/
630 gamma = ddot_(&nr, &work[*ilo + (*n << 2)], &c__1, &work[*ilo + (*n << 2)]
631 , &c__1) + ddot_(&nr, &work[*ilo + *n * 5], &c__1, &work[*ilo + *
632 n * 5], &c__1);
633
634 /*< EW = ZERO >*/
635 ew = 0.;
636 /*< EWC = ZERO >*/
637 ewc = 0.;
638 /*< DO 260 I = ILO, IHI >*/
639 i__1 = *ihi;
640 for (i__ = *ilo; i__ <= i__1; ++i__) {
641 /*< EW = EW + WORK( I+4*N ) >*/
642 ew += work[i__ + (*n << 2)];
643 /*< EWC = EWC + WORK( I+5*N ) >*/
644 ewc += work[i__ + *n * 5];
645 /*< 260 CONTINUE >*/
646 /* L260: */
647 }
648
649 /*< GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2 >*/
650 /* Computing 2nd power */
651 d__1 = ew;
652 /* Computing 2nd power */
653 d__2 = ewc;
654 /* Computing 2nd power */
655 d__3 = ew - ewc;
656 gamma = coef * gamma - coef2 * (d__1 * d__1 + d__2 * d__2) - coef5 * (
657 d__3 * d__3);
658 /*< >*/
659 if (gamma == 0.) {
660 goto L350;
661 }
662 /*< >*/
663 if (it != 1) {
664 beta = gamma / pgamma;
665 }
666 /*< T = COEF5*( EWC-THREE*EW ) >*/
667 t = coef5 * (ewc - ew * 3.);
668 /*< TC = COEF5*( EW-THREE*EWC ) >*/
669 tc = coef5 * (ew - ewc * 3.);
670
671 /*< CALL DSCAL( NR, BETA, WORK( ILO ), 1 ) >*/
672 dscal_(&nr, &beta, &work[*ilo], &c__1);
673 /*< CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 ) >*/
674 dscal_(&nr, &beta, &work[*ilo + *n], &c__1);
675
676 /*< CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 ) >*/
677 daxpy_(&nr, &coef, &work[*ilo + (*n << 2)], &c__1, &work[*ilo + *n], &
678 c__1);
679 /*< CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 ) >*/
680 daxpy_(&nr, &coef, &work[*ilo + *n * 5], &c__1, &work[*ilo], &c__1);
681
682 /*< DO 270 I = ILO, IHI >*/
683 i__1 = *ihi;
684 for (i__ = *ilo; i__ <= i__1; ++i__) {
685 /*< WORK( I ) = WORK( I ) + TC >*/
686 work[i__] += tc;
687 /*< WORK( I+N ) = WORK( I+N ) + T >*/
688 work[i__ + *n] += t;
689 /*< 270 CONTINUE >*/
690 /* L270: */
691 }
692
693 /* Apply matrix to vector */
694
695 /*< DO 300 I = ILO, IHI >*/
696 i__1 = *ihi;
697 for (i__ = *ilo; i__ <= i__1; ++i__) {
698 /*< KOUNT = 0 >*/
699 kount = 0;
700 /*< SUM = ZERO >*/
701 sum = 0.;
702 /*< DO 290 J = ILO, IHI >*/
703 i__2 = *ihi;
704 for (j = *ilo; j <= i__2; ++j) {
705 /*< >*/
706 i__3 = i__ + j * a_dim1;
707 if (a[i__3].r == 0. && a[i__3].i == 0.) {
708 goto L280;
709 }
710 /*< KOUNT = KOUNT + 1 >*/
711 ++kount;
712 /*< SUM = SUM + WORK( J ) >*/
713 sum += work[j];
714 /*< 280 CONTINUE >*/
715 L280:
716 /*< >*/
717 i__3 = i__ + j * b_dim1;
718 if (b[i__3].r == 0. && b[i__3].i == 0.) {
719 goto L290;
720 }
721 /*< KOUNT = KOUNT + 1 >*/
722 ++kount;
723 /*< SUM = SUM + WORK( J ) >*/
724 sum += work[j];
725 /*< 290 CONTINUE >*/
726 L290:
727 ;
728 }
729 /*< WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM >*/
730 work[i__ + (*n << 1)] = (doublereal) kount * work[i__ + *n] + sum;
731 /*< 300 CONTINUE >*/
732 /* L300: */
733 }
734
735 /*< DO 330 J = ILO, IHI >*/
736 i__1 = *ihi;
737 for (j = *ilo; j <= i__1; ++j) {
738 /*< KOUNT = 0 >*/
739 kount = 0;
740 /*< SUM = ZERO >*/
741 sum = 0.;
742 /*< DO 320 I = ILO, IHI >*/
743 i__2 = *ihi;
744 for (i__ = *ilo; i__ <= i__2; ++i__) {
745 /*< >*/
746 i__3 = i__ + j * a_dim1;
747 if (a[i__3].r == 0. && a[i__3].i == 0.) {
748 goto L310;
749 }
750 /*< KOUNT = KOUNT + 1 >*/
751 ++kount;
752 /*< SUM = SUM + WORK( I+N ) >*/
753 sum += work[i__ + *n];
754 /*< 310 CONTINUE >*/
755 L310:
756 /*< >*/
757 i__3 = i__ + j * b_dim1;
758 if (b[i__3].r == 0. && b[i__3].i == 0.) {
759 goto L320;
760 }
761 /*< KOUNT = KOUNT + 1 >*/
762 ++kount;
763 /*< SUM = SUM + WORK( I+N ) >*/
764 sum += work[i__ + *n];
765 /*< 320 CONTINUE >*/
766 L320:
767 ;
768 }
769 /*< WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM >*/
770 work[j + *n * 3] = (doublereal) kount * work[j] + sum;
771 /*< 330 CONTINUE >*/
772 /* L330: */
773 }
774
775 /*< >*/
776 sum = ddot_(&nr, &work[*ilo + *n], &c__1, &work[*ilo + (*n << 1)], &c__1)
777 + ddot_(&nr, &work[*ilo], &c__1, &work[*ilo + *n * 3], &c__1);
778 /*< ALPHA = GAMMA / SUM >*/
779 alpha = gamma / sum;
780
781 /* Determine correction to current iteration */
782
783 /*< CMAX = ZERO >*/
784 cmax = 0.;
785 /*< DO 340 I = ILO, IHI >*/
786 i__1 = *ihi;
787 for (i__ = *ilo; i__ <= i__1; ++i__) {
788 /*< COR = ALPHA*WORK( I+N ) >*/
789 cor = alpha * work[i__ + *n];
790 /*< >*/
791 if (abs(cor) > cmax) {
792 cmax = abs(cor);
793 }
794 /*< LSCALE( I ) = LSCALE( I ) + COR >*/
795 lscale[i__] += cor;
796 /*< COR = ALPHA*WORK( I ) >*/
797 cor = alpha * work[i__];
798 /*< >*/
799 if (abs(cor) > cmax) {
800 cmax = abs(cor);
801 }
802 /*< RSCALE( I ) = RSCALE( I ) + COR >*/
803 rscale[i__] += cor;
804 /*< 340 CONTINUE >*/
805 /* L340: */
806 }
807 /*< >*/
808 if (cmax < .5) {
809 goto L350;
810 }
811
812 /*< CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 ) >*/
813 d__1 = -alpha;
814 daxpy_(&nr, &d__1, &work[*ilo + (*n << 1)], &c__1, &work[*ilo + (*n << 2)]
815 , &c__1);
816 /*< CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 ) >*/
817 d__1 = -alpha;
818 daxpy_(&nr, &d__1, &work[*ilo + *n * 3], &c__1, &work[*ilo + *n * 5], &
819 c__1);
820
821 /*< PGAMMA = GAMMA >*/
822 pgamma = gamma;
823 /*< IT = IT + 1 >*/
824 ++it;
825 /*< >*/
826 if (it <= nrp2) {
827 goto L250;
828 }
829
830 /* End generalized conjugate gradient iteration */
831
832 /*< 350 CONTINUE >*/
833 L350:
834 /*< SFMIN = DLAMCH( 'S' ) >*/
835 sfmin = dlamch_("S", (ftnlen)1);
836 /*< SFMAX = ONE / SFMIN >*/
837 sfmax = 1. / sfmin;
838 /*< LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE ) >*/
839 lsfmin = (integer) (d_lg10(&sfmin) / basl + 1.);
840 /*< LSFMAX = INT( LOG10( SFMAX ) / BASL ) >*/
841 lsfmax = (integer) (d_lg10(&sfmax) / basl);
842 /*< DO 360 I = ILO, IHI >*/
843 i__1 = *ihi;
844 for (i__ = *ilo; i__ <= i__1; ++i__) {
845 /*< IRAB = IZAMAX( N-ILO+1, A( I, ILO ), LDA ) >*/
846 i__2 = *n - *ilo + 1;
847 irab = izamax_(&i__2, &a[i__ + *ilo * a_dim1], lda);
848 /*< RAB = ABS( A( I, IRAB+ILO-1 ) ) >*/
849 rab = z_abs(&a[i__ + (irab + *ilo - 1) * a_dim1]);
850 /*< IRAB = IZAMAX( N-ILO+1, B( I, ILO ), LDB ) >*/
851 i__2 = *n - *ilo + 1;
852 irab = izamax_(&i__2, &b[i__ + *ilo * b_dim1], ldb);
853 /*< RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) ) >*/
854 /* Computing MAX */
855 d__1 = rab, d__2 = z_abs(&b[i__ + (irab + *ilo - 1) * b_dim1]);
856 rab = max(d__1,d__2);
857 /*< LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE ) >*/
858 d__1 = rab + sfmin;
859 lrab = (integer) (d_lg10(&d__1) / basl + 1.);
860 /*< IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) ) >*/
861 ir = (integer) (lscale[i__] + d_sign(&c_b72, &lscale[i__]));
862 /*< IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB ) >*/
863 /* Computing MIN */
864 i__2 = max(ir,lsfmin), i__2 = min(i__2,lsfmax), i__3 = lsfmax - lrab;
865 ir = min(i__2,i__3);
866 /*< LSCALE( I ) = SCLFAC**IR >*/
867 lscale[i__] = pow_di(&c_b36, &ir);
868 /*< ICAB = IZAMAX( IHI, A( 1, I ), 1 ) >*/
869 icab = izamax_(ihi, &a[i__ * a_dim1 + 1], &c__1);
870 /*< CAB = ABS( A( ICAB, I ) ) >*/
871 cab = z_abs(&a[icab + i__ * a_dim1]);
872 /*< ICAB = IZAMAX( IHI, B( 1, I ), 1 ) >*/
873 icab = izamax_(ihi, &b[i__ * b_dim1 + 1], &c__1);
874 /*< CAB = MAX( CAB, ABS( B( ICAB, I ) ) ) >*/
875 /* Computing MAX */
876 d__1 = cab, d__2 = z_abs(&b[icab + i__ * b_dim1]);
877 cab = max(d__1,d__2);
878 /*< LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE ) >*/
879 d__1 = cab + sfmin;
880 lcab = (integer) (d_lg10(&d__1) / basl + 1.);
881 /*< JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) ) >*/
882 jc = (integer) (rscale[i__] + d_sign(&c_b72, &rscale[i__]));
883 /*< JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB ) >*/
884 /* Computing MIN */
885 i__2 = max(jc,lsfmin), i__2 = min(i__2,lsfmax), i__3 = lsfmax - lcab;
886 jc = min(i__2,i__3);
887 /*< RSCALE( I ) = SCLFAC**JC >*/
888 rscale[i__] = pow_di(&c_b36, &jc);
889 /*< 360 CONTINUE >*/
890 /* L360: */
891 }
892
893 /* Row scaling of matrices A and B */
894
895 /*< DO 370 I = ILO, IHI >*/
896 i__1 = *ihi;
897 for (i__ = *ilo; i__ <= i__1; ++i__) {
898 /*< CALL ZDSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA ) >*/
899 i__2 = *n - *ilo + 1;
900 zdscal_(&i__2, &lscale[i__], &a[i__ + *ilo * a_dim1], lda);
901 /*< CALL ZDSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB ) >*/
902 i__2 = *n - *ilo + 1;
903 zdscal_(&i__2, &lscale[i__], &b[i__ + *ilo * b_dim1], ldb);
904 /*< 370 CONTINUE >*/
905 /* L370: */
906 }
907
908 /* Column scaling of matrices A and B */
909
910 /*< DO 380 J = ILO, IHI >*/
911 i__1 = *ihi;
912 for (j = *ilo; j <= i__1; ++j) {
913 /*< CALL ZDSCAL( IHI, RSCALE( J ), A( 1, J ), 1 ) >*/
914 zdscal_(ihi, &rscale[j], &a[j * a_dim1 + 1], &c__1);
915 /*< CALL ZDSCAL( IHI, RSCALE( J ), B( 1, J ), 1 ) >*/
916 zdscal_(ihi, &rscale[j], &b[j * b_dim1 + 1], &c__1);
917 /*< 380 CONTINUE >*/
918 /* L380: */
919 }
920
921 /*< RETURN >*/
922 return 0;
923
924 /* End of ZGGBAL */
925
926 /*< END >*/
927 } /* zggbal_ */
928
929 #ifdef __cplusplus
930 }
931 #endif
932