1 SUBROUTINE DGGBAL( 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 A( LDA, * ), B( LDB, * ), LSCALE( * ), 15 $ RSCALE( * ), WORK( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* DGGBAL balances a pair of general real 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) 73* is the scaling factor 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) 84* is the scaling factor applied to column j, then 85* LSCALE(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* .. 111* .. Local Scalars .. 112 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1, 113 $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN, 114 $ M, NR, NRP2 115 DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2, 116 $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX, 117 $ SFMIN, SUM, T, TA, TB, TC 118* .. 119* .. External Functions .. 120 LOGICAL LSAME 121 INTEGER IDAMAX 122 DOUBLE PRECISION DDOT, DLAMCH 123 EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH 124* .. 125* .. External Subroutines .. 126 EXTERNAL DAXPY, DSCAL, DSWAP, XERBLA 127* .. 128* .. Intrinsic Functions .. 129 INTRINSIC ABS, DBLE, INT, LOG10, MAX, MIN, SIGN 130* .. 131* .. Executable Statements .. 132* 133* Test the input parameters 134* 135 INFO = 0 136 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. 137 $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN 138 INFO = -1 139 ELSE IF( N.LT.0 ) THEN 140 INFO = -2 141 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 142 INFO = -4 143 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 144 INFO = -6 145 END IF 146 IF( INFO.NE.0 ) THEN 147 CALL XERBLA( 'DGGBAL', -INFO ) 148 RETURN 149 END IF 150* 151 K = 1 152 L = N 153* 154* Quick return if possible 155* 156 IF( N.EQ.0 ) 157 $ RETURN 158* 159 IF( LSAME( JOB, 'N' ) ) THEN 160 ILO = 1 161 IHI = N 162 DO 10 I = 1, N 163 LSCALE( I ) = ONE 164 RSCALE( I ) = ONE 165 10 CONTINUE 166 RETURN 167 END IF 168* 169 IF( K.EQ.L ) THEN 170 ILO = 1 171 IHI = 1 172 LSCALE( 1 ) = ONE 173 RSCALE( 1 ) = ONE 174 RETURN 175 END IF 176* 177 IF( LSAME( JOB, 'S' ) ) 178 $ GO TO 190 179* 180 GO TO 30 181* 182* Permute the matrices A and B to isolate the eigenvalues. 183* 184* Find row with one nonzero in columns 1 through L 185* 186 20 CONTINUE 187 L = LM1 188 IF( L.NE.1 ) 189 $ GO TO 30 190* 191 RSCALE( 1 ) = ONE 192 LSCALE( 1 ) = ONE 193 GO TO 190 194* 195 30 CONTINUE 196 LM1 = L - 1 197 DO 80 I = L, 1, -1 198 DO 40 J = 1, LM1 199 JP1 = J + 1 200 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 201 $ GO TO 50 202 40 CONTINUE 203 J = L 204 GO TO 70 205* 206 50 CONTINUE 207 DO 60 J = JP1, L 208 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 209 $ GO TO 80 210 60 CONTINUE 211 J = JP1 - 1 212* 213 70 CONTINUE 214 M = L 215 IFLOW = 1 216 GO TO 160 217 80 CONTINUE 218 GO TO 100 219* 220* Find column with one nonzero in rows K through N 221* 222 90 CONTINUE 223 K = K + 1 224* 225 100 CONTINUE 226 DO 150 J = K, L 227 DO 110 I = K, LM1 228 IP1 = I + 1 229 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 230 $ GO TO 120 231 110 CONTINUE 232 I = L 233 GO TO 140 234 120 CONTINUE 235 DO 130 I = IP1, L 236 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 237 $ GO TO 150 238 130 CONTINUE 239 I = IP1 - 1 240 140 CONTINUE 241 M = K 242 IFLOW = 2 243 GO TO 160 244 150 CONTINUE 245 GO TO 190 246* 247* Permute rows M and I 248* 249 160 CONTINUE 250 LSCALE( M ) = DBLE( I ) 251 IF( I.EQ.M ) 252 $ GO TO 170 253 CALL DSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA ) 254 CALL DSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB ) 255* 256* Permute columns M and J 257* 258 170 CONTINUE 259 RSCALE( M ) = DBLE( J ) 260 IF( J.EQ.M ) 261 $ GO TO 180 262 CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 ) 263 CALL DSWAP( L, B( 1, J ), 1, B( 1, M ), 1 ) 264* 265 180 CONTINUE 266 GO TO ( 20, 90 )IFLOW 267* 268 190 CONTINUE 269 ILO = K 270 IHI = L 271* 272 IF( ILO.EQ.IHI ) 273 $ RETURN 274* 275 IF( LSAME( JOB, 'P' ) ) 276 $ RETURN 277* 278* Balance the submatrix in rows ILO to IHI. 279* 280 NR = IHI - ILO + 1 281 DO 200 I = ILO, IHI 282 RSCALE( I ) = ZERO 283 LSCALE( I ) = ZERO 284* 285 WORK( I ) = ZERO 286 WORK( I+N ) = ZERO 287 WORK( I+2*N ) = ZERO 288 WORK( I+3*N ) = ZERO 289 WORK( I+4*N ) = ZERO 290 WORK( I+5*N ) = ZERO 291 200 CONTINUE 292* 293* Compute right side vector in resulting linear equations 294* 295 BASL = LOG10( SCLFAC ) 296 DO 240 I = ILO, IHI 297 DO 230 J = ILO, IHI 298 TB = B( I, J ) 299 TA = A( I, J ) 300 IF( TA.EQ.ZERO ) 301 $ GO TO 210 302 TA = LOG10( ABS( TA ) ) / BASL 303 210 CONTINUE 304 IF( TB.EQ.ZERO ) 305 $ GO TO 220 306 TB = LOG10( ABS( TB ) ) / BASL 307 220 CONTINUE 308 WORK( I+4*N ) = WORK( I+4*N ) - TA - TB 309 WORK( J+5*N ) = WORK( J+5*N ) - TA - TB 310 230 CONTINUE 311 240 CONTINUE 312* 313 COEF = ONE / DBLE( 2*NR ) 314 COEF2 = COEF*COEF 315 COEF5 = HALF*COEF2 316 NRP2 = NR + 2 317 BETA = ZERO 318 IT = 1 319* 320* Start generalized conjugate gradient iteration 321* 322 250 CONTINUE 323* 324 GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) + 325 $ DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 ) 326* 327 EW = ZERO 328 EWC = ZERO 329 DO 260 I = ILO, IHI 330 EW = EW + WORK( I+4*N ) 331 EWC = EWC + WORK( I+5*N ) 332 260 CONTINUE 333* 334 GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2 335 IF( GAMMA.EQ.ZERO ) 336 $ GO TO 350 337 IF( IT.NE.1 ) 338 $ BETA = GAMMA / PGAMMA 339 T = COEF5*( EWC-THREE*EW ) 340 TC = COEF5*( EW-THREE*EWC ) 341* 342 CALL DSCAL( NR, BETA, WORK( ILO ), 1 ) 343 CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 ) 344* 345 CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 ) 346 CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 ) 347* 348 DO 270 I = ILO, IHI 349 WORK( I ) = WORK( I ) + TC 350 WORK( I+N ) = WORK( I+N ) + T 351 270 CONTINUE 352* 353* Apply matrix to vector 354* 355 DO 300 I = ILO, IHI 356 KOUNT = 0 357 SUM = ZERO 358 DO 290 J = ILO, IHI 359 IF( A( I, J ).EQ.ZERO ) 360 $ GO TO 280 361 KOUNT = KOUNT + 1 362 SUM = SUM + WORK( J ) 363 280 CONTINUE 364 IF( B( I, J ).EQ.ZERO ) 365 $ GO TO 290 366 KOUNT = KOUNT + 1 367 SUM = SUM + WORK( J ) 368 290 CONTINUE 369 WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM 370 300 CONTINUE 371* 372 DO 330 J = ILO, IHI 373 KOUNT = 0 374 SUM = ZERO 375 DO 320 I = ILO, IHI 376 IF( A( I, J ).EQ.ZERO ) 377 $ GO TO 310 378 KOUNT = KOUNT + 1 379 SUM = SUM + WORK( I+N ) 380 310 CONTINUE 381 IF( B( I, J ).EQ.ZERO ) 382 $ GO TO 320 383 KOUNT = KOUNT + 1 384 SUM = SUM + WORK( I+N ) 385 320 CONTINUE 386 WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM 387 330 CONTINUE 388* 389 SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) + 390 $ DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 ) 391 ALPHA = GAMMA / SUM 392* 393* Determine correction to current iteration 394* 395 CMAX = ZERO 396 DO 340 I = ILO, IHI 397 COR = ALPHA*WORK( I+N ) 398 IF( ABS( COR ).GT.CMAX ) 399 $ CMAX = ABS( COR ) 400 LSCALE( I ) = LSCALE( I ) + COR 401 COR = ALPHA*WORK( I ) 402 IF( ABS( COR ).GT.CMAX ) 403 $ CMAX = ABS( COR ) 404 RSCALE( I ) = RSCALE( I ) + COR 405 340 CONTINUE 406 IF( CMAX.LT.HALF ) 407 $ GO TO 350 408* 409 CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 ) 410 CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 ) 411* 412 PGAMMA = GAMMA 413 IT = IT + 1 414 IF( IT.LE.NRP2 ) 415 $ GO TO 250 416* 417* End generalized conjugate gradient iteration 418* 419 350 CONTINUE 420 SFMIN = DLAMCH( 'S' ) 421 SFMAX = ONE / SFMIN 422 LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE ) 423 LSFMAX = INT( LOG10( SFMAX ) / BASL ) 424 DO 360 I = ILO, IHI 425 IRAB = IDAMAX( N-ILO+1, A( I, ILO ), LDA ) 426 RAB = ABS( A( I, IRAB+ILO-1 ) ) 427 IRAB = IDAMAX( N-ILO+1, B( I, ILO ), LDB ) 428 RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) ) 429 LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE ) 430 IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) ) 431 IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB ) 432 LSCALE( I ) = SCLFAC**IR 433 ICAB = IDAMAX( IHI, A( 1, I ), 1 ) 434 CAB = ABS( A( ICAB, I ) ) 435 ICAB = IDAMAX( IHI, B( 1, I ), 1 ) 436 CAB = MAX( CAB, ABS( B( ICAB, I ) ) ) 437 LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE ) 438 JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) ) 439 JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB ) 440 RSCALE( I ) = SCLFAC**JC 441 360 CONTINUE 442* 443* Row scaling of matrices A and B 444* 445 DO 370 I = ILO, IHI 446 CALL DSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA ) 447 CALL DSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB ) 448 370 CONTINUE 449* 450* Column scaling of matrices A and B 451* 452 DO 380 J = ILO, IHI 453 CALL DSCAL( IHI, RSCALE( J ), A( 1, J ), 1 ) 454 CALL DSCAL( IHI, RSCALE( J ), B( 1, J ), 1 ) 455 380 CONTINUE 456* 457 RETURN 458* 459* End of DGGBAL 460* 461 END 462