1
SUBROUTINE ZGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
4
* -- LAPACK routine (version 3.1) --
5
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8
* .. Scalar Arguments ..
10
INTEGER IHI, ILO, INFO, LDA, LDB, N
12
* .. Array Arguments ..
13
DOUBLE PRECISION LSCALE( * ), RSCALE( * ), WORK( * )
14
COMPLEX*16 A( LDA, * ), B( LDB, * )
20
* ZGGBAL balances a pair of general complex matrices (A,B). This
21
* involves, first, permuting A and B by similarity transformations to
22
* isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
23
* elements on the diagonal; and second, applying a diagonal similarity
24
* transformation to rows and columns ILO to IHI to make the rows
25
* and columns as close in norm as possible. Both steps are optional.
27
* Balancing may reduce the 1-norm of the matrices, and improve the
28
* accuracy of the computed eigenvalues and/or eigenvectors in the
29
* generalized eigenvalue problem A*x = lambda*B*x.
34
* JOB (input) CHARACTER*1
35
* Specifies the operations to be performed on A and B:
36
* = 'N': none: simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
37
* and RSCALE(I) = 1.0 for i=1,...,N;
38
* = 'P': permute only;
40
* = 'B': both permute and scale.
43
* The order of the matrices A and B. N >= 0.
45
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
46
* On entry, the input matrix A.
47
* On exit, A is overwritten by the balanced matrix.
48
* If JOB = 'N', A is not referenced.
51
* The leading dimension of the array A. LDA >= max(1,N).
53
* B (input/output) COMPLEX*16 array, dimension (LDB,N)
54
* On entry, the input matrix B.
55
* On exit, B is overwritten by the balanced matrix.
56
* If JOB = 'N', B is not referenced.
59
* The leading dimension of the array B. LDB >= max(1,N).
61
* ILO (output) INTEGER
62
* IHI (output) INTEGER
63
* ILO and IHI are set to integers such that on exit
64
* A(i,j) = 0 and B(i,j) = 0 if i > j and
65
* j = 1,...,ILO-1 or i = IHI+1,...,N.
66
* If JOB = 'N' or 'S', ILO = 1 and IHI = N.
68
* LSCALE (output) DOUBLE PRECISION array, dimension (N)
69
* Details of the permutations and scaling factors applied
70
* to the left side of A and B. If P(j) is the index of the
71
* row interchanged with row j, and D(j) is the scaling factor
72
* applied to row j, then
73
* LSCALE(j) = P(j) for J = 1,...,ILO-1
74
* = D(j) for J = ILO,...,IHI
75
* = P(j) for J = IHI+1,...,N.
76
* The order in which the interchanges are made is N to IHI+1,
79
* RSCALE (output) DOUBLE PRECISION array, dimension (N)
80
* Details of the permutations and scaling factors applied
81
* to the right side of A and B. If P(j) is the index of the
82
* column interchanged with column j, and D(j) is the scaling
83
* factor applied to column j, then
84
* RSCALE(j) = P(j) for J = 1,...,ILO-1
85
* = D(j) for J = ILO,...,IHI
86
* = P(j) for J = IHI+1,...,N.
87
* The order in which the interchanges are made is N to IHI+1,
90
* WORK (workspace) REAL array, dimension (lwork)
91
* lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
92
* at least 1 when JOB = 'N' or 'P'.
94
* INFO (output) INTEGER
95
* = 0: successful exit
96
* < 0: if INFO = -i, the i-th argument had an illegal value.
101
* See R.C. WARD, Balancing the generalized eigenvalue problem,
102
* SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
104
* =====================================================================
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 )
112
PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
114
* .. Local Scalars ..
115
INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
116
$ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
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
123
* .. External Functions ..
126
DOUBLE PRECISION DDOT, DLAMCH
127
EXTERNAL LSAME, IZAMAX, DDOT, DLAMCH
129
* .. External Subroutines ..
130
EXTERNAL DAXPY, DSCAL, XERBLA, ZDSCAL, ZSWAP
132
* .. Intrinsic Functions ..
133
INTRINSIC ABS, DBLE, DIMAG, INT, LOG10, MAX, MIN, SIGN
135
* .. Statement Functions ..
136
DOUBLE PRECISION CABS1
138
* .. Statement Function definitions ..
139
CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
141
* .. Executable Statements ..
143
* Test the input parameters
146
IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
147
$ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
149
ELSE IF( N.LT.0 ) THEN
151
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
153
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
157
CALL XERBLA( 'ZGGBAL', -INFO )
161
* Quick return if possible
177
IF( LSAME( JOB, 'N' ) ) THEN
189
IF( LSAME( JOB, 'S' ) )
194
* Permute the matrices A and B to isolate the eigenvalues.
196
* Find row with one nonzero in columns 1 through L
212
IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
220
IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
232
* Find column with one nonzero in rows K through N
241
IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
248
IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
259
* Permute rows M and I
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 )
268
* Permute columns M and J
274
CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
275
CALL ZSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
278
GO TO ( 20, 90 )IFLOW
284
IF( LSAME( JOB, 'P' ) ) THEN
295
* Balance the submatrix in rows ILO to IHI.
310
* Compute right side vector in resulting linear equations
312
BASL = LOG10( SCLFAC )
315
IF( A( I, J ).EQ.CZERO ) THEN
319
TA = LOG10( CABS1( A( I, J ) ) ) / BASL
322
IF( B( I, J ).EQ.CZERO ) THEN
326
TB = LOG10( CABS1( B( I, J ) ) ) / BASL
329
WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
330
WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
334
COEF = ONE / DBLE( 2*NR )
341
* Start generalized conjugate gradient iteration
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 )
351
EW = EW + WORK( I+4*N )
352
EWC = EWC + WORK( I+5*N )
355
GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
359
$ BETA = GAMMA / PGAMMA
360
T = COEF5*( EWC-THREE*EW )
361
TC = COEF5*( EW-THREE*EWC )
363
CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
364
CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
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 )
370
WORK( I ) = WORK( I ) + TC
371
WORK( I+N ) = WORK( I+N ) + T
374
* Apply matrix to vector
380
IF( A( I, J ).EQ.CZERO )
383
SUM = SUM + WORK( J )
385
IF( B( I, J ).EQ.CZERO )
388
SUM = SUM + WORK( J )
390
WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
397
IF( A( I, J ).EQ.CZERO )
400
SUM = SUM + WORK( I+N )
402
IF( B( I, J ).EQ.CZERO )
405
SUM = SUM + WORK( I+N )
407
WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
410
SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
411
$ DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
414
* Determine correction to current iteration
418
COR = ALPHA*WORK( I+N )
419
IF( ABS( COR ).GT.CMAX )
421
LSCALE( I ) = LSCALE( I ) + COR
422
COR = ALPHA*WORK( I )
423
IF( ABS( COR ).GT.CMAX )
425
RSCALE( I ) = RSCALE( I ) + COR
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 )
438
* End generalized conjugate gradient iteration
441
SFMIN = DLAMCH( 'S' )
443
LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
444
LSFMAX = INT( LOG10( SFMAX ) / BASL )
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
464
* Row scaling of matrices A and B
467
CALL ZDSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
468
CALL ZDSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
471
* Column scaling of matrices A and B
474
CALL ZDSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
475
CALL ZDSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
added
removed
libcruft/lapack/zggbal.f
