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