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