001:       SUBROUTINE ZGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
002:      $                   RSCALE, WORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          JOB
010:       INTEGER            IHI, ILO, INFO, LDA, LDB, N
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   LSCALE( * ), RSCALE( * ), WORK( * )
014:       COMPLEX*16         A( LDA, * ), B( LDB, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  ZGGBAL balances a pair of general complex matrices (A,B).  This
021: *  involves, first, permuting A and B by similarity transformations to
022: *  isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
023: *  elements on the diagonal; and second, applying a diagonal similarity
024: *  transformation to rows and columns ILO to IHI to make the rows
025: *  and columns as close in norm as possible. Both steps are optional.
026: *
027: *  Balancing may reduce the 1-norm of the matrices, and improve the
028: *  accuracy of the computed eigenvalues and/or eigenvectors in the
029: *  generalized eigenvalue problem A*x = lambda*B*x.
030: *
031: *  Arguments
032: *  =========
033: *
034: *  JOB     (input) CHARACTER*1
035: *          Specifies the operations to be performed on A and B:
036: *          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
037: *                  and RSCALE(I) = 1.0 for i=1,...,N;
038: *          = 'P':  permute only;
039: *          = 'S':  scale only;
040: *          = 'B':  both permute and scale.
041: *
042: *  N       (input) INTEGER
043: *          The order of the matrices A and B.  N >= 0.
044: *
045: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
046: *          On entry, the input matrix A.
047: *          On exit, A is overwritten by the balanced matrix.
048: *          If JOB = 'N', A is not referenced.
049: *
050: *  LDA     (input) INTEGER
051: *          The leading dimension of the array A. LDA >= max(1,N).
052: *
053: *  B       (input/output) COMPLEX*16 array, dimension (LDB,N)
054: *          On entry, the input matrix B.
055: *          On exit, B is overwritten by the balanced matrix.
056: *          If JOB = 'N', B is not referenced.
057: *
058: *  LDB     (input) INTEGER
059: *          The leading dimension of the array B. LDB >= max(1,N).
060: *
061: *  ILO     (output) INTEGER
062: *  IHI     (output) INTEGER
063: *          ILO and IHI are set to integers such that on exit
064: *          A(i,j) = 0 and B(i,j) = 0 if i > j and
065: *          j = 1,...,ILO-1 or i = IHI+1,...,N.
066: *          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
067: *
068: *  LSCALE  (output) DOUBLE PRECISION array, dimension (N)
069: *          Details of the permutations and scaling factors applied
070: *          to the left side of A and B.  If P(j) is the index of the
071: *          row interchanged with row j, and D(j) is the scaling factor
072: *          applied to row j, then
073: *            LSCALE(j) = P(j)    for J = 1,...,ILO-1
074: *                      = D(j)    for J = ILO,...,IHI
075: *                      = P(j)    for J = IHI+1,...,N.
076: *          The order in which the interchanges are made is N to IHI+1,
077: *          then 1 to ILO-1.
078: *
079: *  RSCALE  (output) DOUBLE PRECISION array, dimension (N)
080: *          Details of the permutations and scaling factors applied
081: *          to the right side of A and B.  If P(j) is the index of the
082: *          column interchanged with column j, and D(j) is the scaling
083: *          factor applied to column j, then
084: *            RSCALE(j) = P(j)    for J = 1,...,ILO-1
085: *                      = D(j)    for J = ILO,...,IHI
086: *                      = P(j)    for J = IHI+1,...,N.
087: *          The order in which the interchanges are made is N to IHI+1,
088: *          then 1 to ILO-1.
089: *
090: *  WORK    (workspace) REAL array, dimension (lwork)
091: *          lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
092: *          at least 1 when JOB = 'N' or 'P'.
093: *
094: *  INFO    (output) INTEGER
095: *          = 0:  successful exit
096: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
097: *
098: *  Further Details
099: *  ===============
100: *
101: *  See R.C. WARD, Balancing the generalized eigenvalue problem,
102: *                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
103: *
104: *  =====================================================================
105: *
106: *     .. Parameters ..
107:       DOUBLE PRECISION   ZERO, HALF, ONE
108:       PARAMETER          ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
109:       DOUBLE PRECISION   THREE, SCLFAC
110:       PARAMETER          ( THREE = 3.0D+0, SCLFAC = 1.0D+1 )
111:       COMPLEX*16         CZERO
112:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
113: *     ..
114: *     .. Local Scalars ..
115:       INTEGER            I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
116:      $                   K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
117:      $                   M, NR, NRP2
118:       DOUBLE PRECISION   ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
119:      $                   COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
120:      $                   SFMIN, SUM, T, TA, TB, TC
121:       COMPLEX*16         CDUM
122: *     ..
123: *     .. External Functions ..
124:       LOGICAL            LSAME
125:       INTEGER            IZAMAX
126:       DOUBLE PRECISION   DDOT, DLAMCH
127:       EXTERNAL           LSAME, IZAMAX, DDOT, DLAMCH
128: *     ..
129: *     .. External Subroutines ..
130:       EXTERNAL           DAXPY, DSCAL, XERBLA, ZDSCAL, ZSWAP
131: *     ..
132: *     .. Intrinsic Functions ..
133:       INTRINSIC          ABS, DBLE, DIMAG, INT, LOG10, MAX, MIN, SIGN
134: *     ..
135: *     .. Statement Functions ..
136:       DOUBLE PRECISION   CABS1
137: *     ..
138: *     .. Statement Function definitions ..
139:       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
140: *     ..
141: *     .. Executable Statements ..
142: *
143: *     Test the input parameters
144: *
145:       INFO = 0
146:       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
147:      $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
148:          INFO = -1
149:       ELSE IF( N.LT.0 ) THEN
150:          INFO = -2
151:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
152:          INFO = -4
153:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
154:          INFO = -6
155:       END IF
156:       IF( INFO.NE.0 ) THEN
157:          CALL XERBLA( 'ZGGBAL', -INFO )
158:          RETURN
159:       END IF
160: *
161: *     Quick return if possible
162: *
163:       IF( N.EQ.0 ) THEN
164:          ILO = 1
165:          IHI = N
166:          RETURN
167:       END IF
168: *
169:       IF( N.EQ.1 ) THEN
170:          ILO = 1
171:          IHI = N
172:          LSCALE( 1 ) = ONE
173:          RSCALE( 1 ) = ONE
174:          RETURN
175:       END IF
176: *
177:       IF( LSAME( JOB, 'N' ) ) THEN
178:          ILO = 1
179:          IHI = N
180:          DO 10 I = 1, N
181:             LSCALE( I ) = ONE
182:             RSCALE( I ) = ONE
183:    10    CONTINUE
184:          RETURN
185:       END IF
186: *
187:       K = 1
188:       L = N
189:       IF( LSAME( JOB, 'S' ) )
190:      $   GO TO 190
191: *
192:       GO TO 30
193: *
194: *     Permute the matrices A and B to isolate the eigenvalues.
195: *
196: *     Find row with one nonzero in columns 1 through L
197: *
198:    20 CONTINUE
199:       L = LM1
200:       IF( L.NE.1 )
201:      $   GO TO 30
202: *
203:       RSCALE( 1 ) = 1
204:       LSCALE( 1 ) = 1
205:       GO TO 190
206: *
207:    30 CONTINUE
208:       LM1 = L - 1
209:       DO 80 I = L, 1, -1
210:          DO 40 J = 1, LM1
211:             JP1 = J + 1
212:             IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
213:      $         GO TO 50
214:    40    CONTINUE
215:          J = L
216:          GO TO 70
217: *
218:    50    CONTINUE
219:          DO 60 J = JP1, L
220:             IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
221:      $         GO TO 80
222:    60    CONTINUE
223:          J = JP1 - 1
224: *
225:    70    CONTINUE
226:          M = L
227:          IFLOW = 1
228:          GO TO 160
229:    80 CONTINUE
230:       GO TO 100
231: *
232: *     Find column with one nonzero in rows K through N
233: *
234:    90 CONTINUE
235:       K = K + 1
236: *
237:   100 CONTINUE
238:       DO 150 J = K, L
239:          DO 110 I = K, LM1
240:             IP1 = I + 1
241:             IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
242:      $         GO TO 120
243:   110    CONTINUE
244:          I = L
245:          GO TO 140
246:   120    CONTINUE
247:          DO 130 I = IP1, L
248:             IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
249:      $         GO TO 150
250:   130    CONTINUE
251:          I = IP1 - 1
252:   140    CONTINUE
253:          M = K
254:          IFLOW = 2
255:          GO TO 160
256:   150 CONTINUE
257:       GO TO 190
258: *
259: *     Permute rows M and I
260: *
261:   160 CONTINUE
262:       LSCALE( M ) = I
263:       IF( I.EQ.M )
264:      $   GO TO 170
265:       CALL ZSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
266:       CALL ZSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
267: *
268: *     Permute columns M and J
269: *
270:   170 CONTINUE
271:       RSCALE( M ) = J
272:       IF( J.EQ.M )
273:      $   GO TO 180
274:       CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
275:       CALL ZSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
276: *
277:   180 CONTINUE
278:       GO TO ( 20, 90 )IFLOW
279: *
280:   190 CONTINUE
281:       ILO = K
282:       IHI = L
283: *
284:       IF( LSAME( JOB, 'P' ) ) THEN
285:          DO 195 I = ILO, IHI
286:             LSCALE( I ) = ONE
287:             RSCALE( I ) = ONE
288:   195    CONTINUE
289:          RETURN
290:       END IF
291: *
292:       IF( ILO.EQ.IHI )
293:      $   RETURN
294: *
295: *     Balance the submatrix in rows ILO to IHI.
296: *
297:       NR = IHI - ILO + 1
298:       DO 200 I = ILO, IHI
299:          RSCALE( I ) = ZERO
300:          LSCALE( I ) = ZERO
301: *
302:          WORK( I ) = ZERO
303:          WORK( I+N ) = ZERO
304:          WORK( I+2*N ) = ZERO
305:          WORK( I+3*N ) = ZERO
306:          WORK( I+4*N ) = ZERO
307:          WORK( I+5*N ) = ZERO
308:   200 CONTINUE
309: *
310: *     Compute right side vector in resulting linear equations
311: *
312:       BASL = LOG10( SCLFAC )
313:       DO 240 I = ILO, IHI
314:          DO 230 J = ILO, IHI
315:             IF( A( I, J ).EQ.CZERO ) THEN
316:                TA = ZERO
317:                GO TO 210
318:             END IF
319:             TA = LOG10( CABS1( A( I, J ) ) ) / BASL
320: *
321:   210       CONTINUE
322:             IF( B( I, J ).EQ.CZERO ) THEN
323:                TB = ZERO
324:                GO TO 220
325:             END IF
326:             TB = LOG10( CABS1( B( I, J ) ) ) / BASL
327: *
328:   220       CONTINUE
329:             WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
330:             WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
331:   230    CONTINUE
332:   240 CONTINUE
333: *
334:       COEF = ONE / DBLE( 2*NR )
335:       COEF2 = COEF*COEF
336:       COEF5 = HALF*COEF2
337:       NRP2 = NR + 2
338:       BETA = ZERO
339:       IT = 1
340: *
341: *     Start generalized conjugate gradient iteration
342: *
343:   250 CONTINUE
344: *
345:       GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
346:      $        DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
347: *
348:       EW = ZERO
349:       EWC = ZERO
350:       DO 260 I = ILO, IHI
351:          EW = EW + WORK( I+4*N )
352:          EWC = EWC + WORK( I+5*N )
353:   260 CONTINUE
354: *
355:       GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
356:       IF( GAMMA.EQ.ZERO )
357:      $   GO TO 350
358:       IF( IT.NE.1 )
359:      $   BETA = GAMMA / PGAMMA
360:       T = COEF5*( EWC-THREE*EW )
361:       TC = COEF5*( EW-THREE*EWC )
362: *
363:       CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
364:       CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
365: *
366:       CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
367:       CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
368: *
369:       DO 270 I = ILO, IHI
370:          WORK( I ) = WORK( I ) + TC
371:          WORK( I+N ) = WORK( I+N ) + T
372:   270 CONTINUE
373: *
374: *     Apply matrix to vector
375: *
376:       DO 300 I = ILO, IHI
377:          KOUNT = 0
378:          SUM = ZERO
379:          DO 290 J = ILO, IHI
380:             IF( A( I, J ).EQ.CZERO )
381:      $         GO TO 280
382:             KOUNT = KOUNT + 1
383:             SUM = SUM + WORK( J )
384:   280       CONTINUE
385:             IF( B( I, J ).EQ.CZERO )
386:      $         GO TO 290
387:             KOUNT = KOUNT + 1
388:             SUM = SUM + WORK( J )
389:   290    CONTINUE
390:          WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
391:   300 CONTINUE
392: *
393:       DO 330 J = ILO, IHI
394:          KOUNT = 0
395:          SUM = ZERO
396:          DO 320 I = ILO, IHI
397:             IF( A( I, J ).EQ.CZERO )
398:      $         GO TO 310
399:             KOUNT = KOUNT + 1
400:             SUM = SUM + WORK( I+N )
401:   310       CONTINUE
402:             IF( B( I, J ).EQ.CZERO )
403:      $         GO TO 320
404:             KOUNT = KOUNT + 1
405:             SUM = SUM + WORK( I+N )
406:   320    CONTINUE
407:          WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
408:   330 CONTINUE
409: *
410:       SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
411:      $      DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
412:       ALPHA = GAMMA / SUM
413: *
414: *     Determine correction to current iteration
415: *
416:       CMAX = ZERO
417:       DO 340 I = ILO, IHI
418:          COR = ALPHA*WORK( I+N )
419:          IF( ABS( COR ).GT.CMAX )
420:      $      CMAX = ABS( COR )
421:          LSCALE( I ) = LSCALE( I ) + COR
422:          COR = ALPHA*WORK( I )
423:          IF( ABS( COR ).GT.CMAX )
424:      $      CMAX = ABS( COR )
425:          RSCALE( I ) = RSCALE( I ) + COR
426:   340 CONTINUE
427:       IF( CMAX.LT.HALF )
428:      $   GO TO 350
429: *
430:       CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
431:       CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
432: *
433:       PGAMMA = GAMMA
434:       IT = IT + 1
435:       IF( IT.LE.NRP2 )
436:      $   GO TO 250
437: *
438: *     End generalized conjugate gradient iteration
439: *
440:   350 CONTINUE
441:       SFMIN = DLAMCH( 'S' )
442:       SFMAX = ONE / SFMIN
443:       LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
444:       LSFMAX = INT( LOG10( SFMAX ) / BASL )
445:       DO 360 I = ILO, IHI
446:          IRAB = IZAMAX( N-ILO+1, A( I, ILO ), LDA )
447:          RAB = ABS( A( I, IRAB+ILO-1 ) )
448:          IRAB = IZAMAX( N-ILO+1, B( I, ILO ), LDB )
449:          RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
450:          LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
451:          IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )
452:          IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
453:          LSCALE( I ) = SCLFAC**IR
454:          ICAB = IZAMAX( IHI, A( 1, I ), 1 )
455:          CAB = ABS( A( ICAB, I ) )
456:          ICAB = IZAMAX( IHI, B( 1, I ), 1 )
457:          CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
458:          LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
459:          JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) )
460:          JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
461:          RSCALE( I ) = SCLFAC**JC
462:   360 CONTINUE
463: *
464: *     Row scaling of matrices A and B
465: *
466:       DO 370 I = ILO, IHI
467:          CALL ZDSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
468:          CALL ZDSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
469:   370 CONTINUE
470: *
471: *     Column scaling of matrices A and B
472: *
473:       DO 380 J = ILO, IHI
474:          CALL ZDSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
475:          CALL ZDSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
476:   380 CONTINUE
477: *
478:       RETURN
479: *
480: *     End of ZGGBAL
481: *
482:       END
483: