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