1 SUBROUTINE ZGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, 2 $ RSCALE, WORK, INFO ) 3* 4* -- LAPACK routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* April 12, 2001 8* 9* .. Scalar Arguments .. 10 CHARACTER JOB 11 INTEGER IHI, ILO, INFO, LDA, LDB, N 12* .. 13* .. Array Arguments .. 14 DOUBLE PRECISION LSCALE( * ), RSCALE( * ), WORK( * ) 15 COMPLEX*16 A( LDA, * ), B( LDB, * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* ZGGBAL balances a pair of general complex matrices (A,B). This 22* involves, first, permuting A and B by similarity transformations to 23* isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N 24* elements on the diagonal; and second, applying a diagonal similarity 25* transformation to rows and columns ILO to IHI to make the rows 26* and columns as close in norm as possible. Both steps are optional. 27* 28* Balancing may reduce the 1-norm of the matrices, and improve the 29* accuracy of the computed eigenvalues and/or eigenvectors in the 30* generalized eigenvalue problem A*x = lambda*B*x. 31* 32* Arguments 33* ========= 34* 35* JOB (input) CHARACTER*1 36* Specifies the operations to be performed on A and B: 37* = 'N': none: simply set ILO = 1, IHI = N, LSCALE(I) = 1.0 38* and RSCALE(I) = 1.0 for i=1,...,N; 39* = 'P': permute only; 40* = 'S': scale only; 41* = 'B': both permute and scale. 42* 43* N (input) INTEGER 44* The order of the matrices A and B. N >= 0. 45* 46* A (input/output) COMPLEX*16 array, dimension (LDA,N) 47* On entry, the input matrix A. 48* On exit, A is overwritten by the balanced matrix. 49* If JOB = 'N', A is not referenced. 50* 51* LDA (input) INTEGER 52* The leading dimension of the array A. LDA >= max(1,N). 53* 54* B (input/output) COMPLEX*16 array, dimension (LDB,N) 55* On entry, the input matrix B. 56* On exit, B is overwritten by the balanced matrix. 57* If JOB = 'N', B is not referenced. 58* 59* LDB (input) INTEGER 60* The leading dimension of the array B. LDB >= max(1,N). 61* 62* ILO (output) INTEGER 63* IHI (output) INTEGER 64* ILO and IHI are set to integers such that on exit 65* A(i,j) = 0 and B(i,j) = 0 if i > j and 66* j = 1,...,ILO-1 or i = IHI+1,...,N. 67* If JOB = 'N' or 'S', ILO = 1 and IHI = N. 68* 69* LSCALE (output) DOUBLE PRECISION array, dimension (N) 70* Details of the permutations and scaling factors applied 71* to the left side of A and B. If P(j) is the index of the 72* row interchanged with row j, and D(j) is the scaling factor 73* applied to row j, then 74* LSCALE(j) = P(j) for J = 1,...,ILO-1 75* = D(j) for J = ILO,...,IHI 76* = P(j) for J = IHI+1,...,N. 77* The order in which the interchanges are made is N to IHI+1, 78* then 1 to ILO-1. 79* 80* RSCALE (output) DOUBLE PRECISION array, dimension (N) 81* Details of the permutations and scaling factors applied 82* to the right side of A and B. If P(j) is the index of the 83* column interchanged with column j, and D(j) is the scaling 84* factor applied to column j, then 85* RSCALE(j) = P(j) for J = 1,...,ILO-1 86* = D(j) for J = ILO,...,IHI 87* = P(j) for J = IHI+1,...,N. 88* The order in which the interchanges are made is N to IHI+1, 89* then 1 to ILO-1. 90* 91* WORK (workspace) DOUBLE PRECISION array, dimension (6*N) 92* 93* INFO (output) INTEGER 94* = 0: successful exit 95* < 0: if INFO = -i, the i-th argument had an illegal value. 96* 97* Further Details 98* =============== 99* 100* See R.C. WARD, Balancing the generalized eigenvalue problem, 101* SIAM J. Sci. Stat. Comp. 2 (1981), 141-152. 102* 103* ===================================================================== 104* 105* .. Parameters .. 106 DOUBLE PRECISION ZERO, HALF, ONE 107 PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) 108 DOUBLE PRECISION THREE, SCLFAC 109 PARAMETER ( THREE = 3.0D+0, SCLFAC = 1.0D+1 ) 110 COMPLEX*16 CZERO 111 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) 112* .. 113* .. Local Scalars .. 114 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1, 115 $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN, 116 $ M, NR, NRP2 117 DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2, 118 $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX, 119 $ SFMIN, SUM, T, TA, TB, TC 120 COMPLEX*16 CDUM 121* .. 122* .. External Functions .. 123 LOGICAL LSAME 124 INTEGER IZAMAX 125 DOUBLE PRECISION DDOT, DLAMCH 126 EXTERNAL LSAME, IZAMAX, DDOT, DLAMCH 127* .. 128* .. External Subroutines .. 129 EXTERNAL DAXPY, DSCAL, XERBLA, ZDSCAL, ZSWAP 130* .. 131* .. Intrinsic Functions .. 132 INTRINSIC ABS, DBLE, DIMAG, INT, LOG10, MAX, MIN, SIGN 133* .. 134* .. Statement Functions .. 135 DOUBLE PRECISION CABS1 136* .. 137* .. Statement Function definitions .. 138 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 139* .. 140* .. Executable Statements .. 141* 142* Test the input parameters 143* 144 INFO = 0 145 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. 146 $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN 147 INFO = -1 148 ELSE IF( N.LT.0 ) THEN 149 INFO = -2 150 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 151 INFO = -4 152 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 153 INFO = -6 154 END IF 155 IF( INFO.NE.0 ) THEN 156 CALL XERBLA( 'ZGGBAL', -INFO ) 157 RETURN 158 END IF 159* 160 K = 1 161 L = N 162* 163* Quick return if possible 164* 165 IF( N.EQ.0 ) 166 $ RETURN 167* 168 IF( LSAME( JOB, 'N' ) ) THEN 169 ILO = 1 170 IHI = N 171 DO 10 I = 1, N 172 LSCALE( I ) = ONE 173 RSCALE( I ) = ONE 174 10 CONTINUE 175 RETURN 176 END IF 177* 178 IF( K.EQ.L ) THEN 179 ILO = 1 180 IHI = 1 181 LSCALE( 1 ) = ONE 182 RSCALE( 1 ) = ONE 183 RETURN 184 END IF 185* 186 IF( LSAME( JOB, 'S' ) ) 187 $ GO TO 190 188* 189 GO TO 30 190* 191* Permute the matrices A and B to isolate the eigenvalues. 192* 193* Find row with one nonzero in columns 1 through L 194* 195 20 CONTINUE 196 L = LM1 197 IF( L.NE.1 ) 198 $ GO TO 30 199* 200 RSCALE( 1 ) = ONE 201 LSCALE( 1 ) = ONE 202 GO TO 190 203* 204 30 CONTINUE 205 LM1 = L - 1 206 DO 80 I = L, 1, -1 207 DO 40 J = 1, LM1 208 JP1 = J + 1 209 IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO ) 210 $ GO TO 50 211 40 CONTINUE 212 J = L 213 GO TO 70 214* 215 50 CONTINUE 216 DO 60 J = JP1, L 217 IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO ) 218 $ GO TO 80 219 60 CONTINUE 220 J = JP1 - 1 221* 222 70 CONTINUE 223 M = L 224 IFLOW = 1 225 GO TO 160 226 80 CONTINUE 227 GO TO 100 228* 229* Find column with one nonzero in rows K through N 230* 231 90 CONTINUE 232 K = K + 1 233* 234 100 CONTINUE 235 DO 150 J = K, L 236 DO 110 I = K, LM1 237 IP1 = I + 1 238 IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO ) 239 $ GO TO 120 240 110 CONTINUE 241 I = L 242 GO TO 140 243 120 CONTINUE 244 DO 130 I = IP1, L 245 IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO ) 246 $ GO TO 150 247 130 CONTINUE 248 I = IP1 - 1 249 140 CONTINUE 250 M = K 251 IFLOW = 2 252 GO TO 160 253 150 CONTINUE 254 GO TO 190 255* 256* Permute rows M and I 257* 258 160 CONTINUE 259 LSCALE( M ) = DBLE( I ) 260 IF( I.EQ.M ) 261 $ GO TO 170 262 CALL ZSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA ) 263 CALL ZSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB ) 264* 265* Permute columns M and J 266* 267 170 CONTINUE 268 RSCALE( M ) = DBLE( J ) 269 IF( J.EQ.M ) 270 $ GO TO 180 271 CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 ) 272 CALL ZSWAP( L, B( 1, J ), 1, B( 1, M ), 1 ) 273* 274 180 CONTINUE 275 GO TO ( 20, 90 )IFLOW 276* 277 190 CONTINUE 278 ILO = K 279 IHI = L 280* 281 IF( ILO.EQ.IHI ) 282 $ RETURN 283* 284 IF( LSAME( JOB, 'P' ) ) 285 $ RETURN 286* 287* Balance the submatrix in rows ILO to IHI. 288* 289 NR = IHI - ILO + 1 290 DO 200 I = ILO, IHI 291 RSCALE( I ) = ZERO 292 LSCALE( I ) = ZERO 293* 294 WORK( I ) = ZERO 295 WORK( I+N ) = ZERO 296 WORK( I+2*N ) = ZERO 297 WORK( I+3*N ) = ZERO 298 WORK( I+4*N ) = ZERO 299 WORK( I+5*N ) = ZERO 300 200 CONTINUE 301* 302* Compute right side vector in resulting linear equations 303* 304 BASL = LOG10( SCLFAC ) 305 DO 240 I = ILO, IHI 306 DO 230 J = ILO, IHI 307 IF( A( I, J ).EQ.CZERO ) THEN 308 TA = ZERO 309 GO TO 210 310 END IF 311 TA = LOG10( CABS1( A( I, J ) ) ) / BASL 312* 313 210 CONTINUE 314 IF( B( I, J ).EQ.CZERO ) THEN 315 TB = ZERO 316 GO TO 220 317 END IF 318 TB = LOG10( CABS1( B( I, J ) ) ) / BASL 319* 320 220 CONTINUE 321 WORK( I+4*N ) = WORK( I+4*N ) - TA - TB 322 WORK( J+5*N ) = WORK( J+5*N ) - TA - TB 323 230 CONTINUE 324 240 CONTINUE 325* 326 COEF = ONE / DBLE( 2*NR ) 327 COEF2 = COEF*COEF 328 COEF5 = HALF*COEF2 329 NRP2 = NR + 2 330 BETA = ZERO 331 IT = 1 332* 333* Start generalized conjugate gradient iteration 334* 335 250 CONTINUE 336* 337 GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) + 338 $ DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 ) 339* 340 EW = ZERO 341 EWC = ZERO 342 DO 260 I = ILO, IHI 343 EW = EW + WORK( I+4*N ) 344 EWC = EWC + WORK( I+5*N ) 345 260 CONTINUE 346* 347 GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2 348 IF( GAMMA.EQ.ZERO ) 349 $ GO TO 350 350 IF( IT.NE.1 ) 351 $ BETA = GAMMA / PGAMMA 352 T = COEF5*( EWC-THREE*EW ) 353 TC = COEF5*( EW-THREE*EWC ) 354* 355 CALL DSCAL( NR, BETA, WORK( ILO ), 1 ) 356 CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 ) 357* 358 CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 ) 359 CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 ) 360* 361 DO 270 I = ILO, IHI 362 WORK( I ) = WORK( I ) + TC 363 WORK( I+N ) = WORK( I+N ) + T 364 270 CONTINUE 365* 366* Apply matrix to vector 367* 368 DO 300 I = ILO, IHI 369 KOUNT = 0 370 SUM = ZERO 371 DO 290 J = ILO, IHI 372 IF( A( I, J ).EQ.CZERO ) 373 $ GO TO 280 374 KOUNT = KOUNT + 1 375 SUM = SUM + WORK( J ) 376 280 CONTINUE 377 IF( B( I, J ).EQ.CZERO ) 378 $ GO TO 290 379 KOUNT = KOUNT + 1 380 SUM = SUM + WORK( J ) 381 290 CONTINUE 382 WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM 383 300 CONTINUE 384* 385 DO 330 J = ILO, IHI 386 KOUNT = 0 387 SUM = ZERO 388 DO 320 I = ILO, IHI 389 IF( A( I, J ).EQ.CZERO ) 390 $ GO TO 310 391 KOUNT = KOUNT + 1 392 SUM = SUM + WORK( I+N ) 393 310 CONTINUE 394 IF( B( I, J ).EQ.CZERO ) 395 $ GO TO 320 396 KOUNT = KOUNT + 1 397 SUM = SUM + WORK( I+N ) 398 320 CONTINUE 399 WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM 400 330 CONTINUE 401* 402 SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) + 403 $ DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 ) 404 ALPHA = GAMMA / SUM 405* 406* Determine correction to current iteration 407* 408 CMAX = ZERO 409 DO 340 I = ILO, IHI 410 COR = ALPHA*WORK( I+N ) 411 IF( ABS( COR ).GT.CMAX ) 412 $ CMAX = ABS( COR ) 413 LSCALE( I ) = LSCALE( I ) + COR 414 COR = ALPHA*WORK( I ) 415 IF( ABS( COR ).GT.CMAX ) 416 $ CMAX = ABS( COR ) 417 RSCALE( I ) = RSCALE( I ) + COR 418 340 CONTINUE 419 IF( CMAX.LT.HALF ) 420 $ GO TO 350 421* 422 CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 ) 423 CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 ) 424* 425 PGAMMA = GAMMA 426 IT = IT + 1 427 IF( IT.LE.NRP2 ) 428 $ GO TO 250 429* 430* End generalized conjugate gradient iteration 431* 432 350 CONTINUE 433 SFMIN = DLAMCH( 'S' ) 434 SFMAX = ONE / SFMIN 435 LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE ) 436 LSFMAX = INT( LOG10( SFMAX ) / BASL ) 437 DO 360 I = ILO, IHI 438 IRAB = IZAMAX( N-ILO+1, A( I, ILO ), LDA ) 439 RAB = ABS( A( I, IRAB+ILO-1 ) ) 440 IRAB = IZAMAX( N-ILO+1, B( I, ILO ), LDB ) 441 RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) ) 442 LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE ) 443 IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) ) 444 IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB ) 445 LSCALE( I ) = SCLFAC**IR 446 ICAB = IZAMAX( IHI, A( 1, I ), 1 ) 447 CAB = ABS( A( ICAB, I ) ) 448 ICAB = IZAMAX( IHI, B( 1, I ), 1 ) 449 CAB = MAX( CAB, ABS( B( ICAB, I ) ) ) 450 LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE ) 451 JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) ) 452 JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB ) 453 RSCALE( I ) = SCLFAC**JC 454 360 CONTINUE 455* 456* Row scaling of matrices A and B 457* 458 DO 370 I = ILO, IHI 459 CALL ZDSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA ) 460 CALL ZDSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB ) 461 370 CONTINUE 462* 463* Column scaling of matrices A and B 464* 465 DO 380 J = ILO, IHI 466 CALL ZDSCAL( IHI, RSCALE( J ), A( 1, J ), 1 ) 467 CALL ZDSCAL( IHI, RSCALE( J ), B( 1, J ), 1 ) 468 380 CONTINUE 469* 470 RETURN 471* 472* End of ZGGBAL 473* 474 END 475