1*> \brief \b SGGBAL 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SGGBAL + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggbal.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggbal.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggbal.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SGGBAL( 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 A( LDA, * ), B( LDB, * ), LSCALE( * ), 30* $ RSCALE( * ), WORK( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> SGGBAL balances a pair of general real 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 REAL 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 REAL 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) 119*> is the scaling factor 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) 133*> is the scaling factor applied to column j, then 134*> LSCALE(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 realGBcomputational 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 SGGBAL( 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 A( LDA, * ), B( LDB, * ), LSCALE( * ), 188 $ RSCALE( * ), WORK( * ) 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* .. 199* .. Local Scalars .. 200 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1, 201 $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN, 202 $ M, NR, NRP2 203 REAL ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2, 204 $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX, 205 $ SFMIN, SUM, T, TA, TB, TC 206* .. 207* .. External Functions .. 208 LOGICAL LSAME 209 INTEGER ISAMAX 210 REAL SDOT, SLAMCH 211 EXTERNAL LSAME, ISAMAX, SDOT, SLAMCH 212* .. 213* .. External Subroutines .. 214 EXTERNAL SAXPY, SSCAL, SSWAP, XERBLA 215* .. 216* .. Intrinsic Functions .. 217 INTRINSIC ABS, INT, LOG10, MAX, MIN, REAL, SIGN 218* .. 219* .. Executable Statements .. 220* 221* Test the input parameters 222* 223 INFO = 0 224 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. 225 $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN 226 INFO = -1 227 ELSE IF( N.LT.0 ) THEN 228 INFO = -2 229 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 230 INFO = -4 231 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 232 INFO = -6 233 END IF 234 IF( INFO.NE.0 ) THEN 235 CALL XERBLA( 'SGGBAL', -INFO ) 236 RETURN 237 END IF 238* 239* Quick return if possible 240* 241 IF( N.EQ.0 ) THEN 242 ILO = 1 243 IHI = N 244 RETURN 245 END IF 246* 247 IF( N.EQ.1 ) THEN 248 ILO = 1 249 IHI = N 250 LSCALE( 1 ) = ONE 251 RSCALE( 1 ) = ONE 252 RETURN 253 END IF 254* 255 IF( LSAME( JOB, 'N' ) ) THEN 256 ILO = 1 257 IHI = N 258 DO 10 I = 1, N 259 LSCALE( I ) = ONE 260 RSCALE( I ) = ONE 261 10 CONTINUE 262 RETURN 263 END IF 264* 265 K = 1 266 L = N 267 IF( LSAME( JOB, 'S' ) ) 268 $ GO TO 190 269* 270 GO TO 30 271* 272* Permute the matrices A and B to isolate the eigenvalues. 273* 274* Find row with one nonzero in columns 1 through L 275* 276 20 CONTINUE 277 L = LM1 278 IF( L.NE.1 ) 279 $ GO TO 30 280* 281 RSCALE( 1 ) = ONE 282 LSCALE( 1 ) = ONE 283 GO TO 190 284* 285 30 CONTINUE 286 LM1 = L - 1 287 DO 80 I = L, 1, -1 288 DO 40 J = 1, LM1 289 JP1 = J + 1 290 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 291 $ GO TO 50 292 40 CONTINUE 293 J = L 294 GO TO 70 295* 296 50 CONTINUE 297 DO 60 J = JP1, L 298 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 299 $ GO TO 80 300 60 CONTINUE 301 J = JP1 - 1 302* 303 70 CONTINUE 304 M = L 305 IFLOW = 1 306 GO TO 160 307 80 CONTINUE 308 GO TO 100 309* 310* Find column with one nonzero in rows K through N 311* 312 90 CONTINUE 313 K = K + 1 314* 315 100 CONTINUE 316 DO 150 J = K, L 317 DO 110 I = K, LM1 318 IP1 = I + 1 319 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 320 $ GO TO 120 321 110 CONTINUE 322 I = L 323 GO TO 140 324 120 CONTINUE 325 DO 130 I = IP1, L 326 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 327 $ GO TO 150 328 130 CONTINUE 329 I = IP1 - 1 330 140 CONTINUE 331 M = K 332 IFLOW = 2 333 GO TO 160 334 150 CONTINUE 335 GO TO 190 336* 337* Permute rows M and I 338* 339 160 CONTINUE 340 LSCALE( M ) = I 341 IF( I.EQ.M ) 342 $ GO TO 170 343 CALL SSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA ) 344 CALL SSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB ) 345* 346* Permute columns M and J 347* 348 170 CONTINUE 349 RSCALE( M ) = J 350 IF( J.EQ.M ) 351 $ GO TO 180 352 CALL SSWAP( L, A( 1, J ), 1, A( 1, M ), 1 ) 353 CALL SSWAP( L, B( 1, J ), 1, B( 1, M ), 1 ) 354* 355 180 CONTINUE 356 GO TO ( 20, 90 )IFLOW 357* 358 190 CONTINUE 359 ILO = K 360 IHI = L 361* 362 IF( LSAME( JOB, 'P' ) ) THEN 363 DO 195 I = ILO, IHI 364 LSCALE( I ) = ONE 365 RSCALE( I ) = ONE 366 195 CONTINUE 367 RETURN 368 END IF 369* 370 IF( ILO.EQ.IHI ) 371 $ RETURN 372* 373* Balance the submatrix in rows ILO to IHI. 374* 375 NR = IHI - ILO + 1 376 DO 200 I = ILO, IHI 377 RSCALE( I ) = ZERO 378 LSCALE( I ) = ZERO 379* 380 WORK( I ) = ZERO 381 WORK( I+N ) = ZERO 382 WORK( I+2*N ) = ZERO 383 WORK( I+3*N ) = ZERO 384 WORK( I+4*N ) = ZERO 385 WORK( I+5*N ) = ZERO 386 200 CONTINUE 387* 388* Compute right side vector in resulting linear equations 389* 390 BASL = LOG10( SCLFAC ) 391 DO 240 I = ILO, IHI 392 DO 230 J = ILO, IHI 393 TB = B( I, J ) 394 TA = A( I, J ) 395 IF( TA.EQ.ZERO ) 396 $ GO TO 210 397 TA = LOG10( ABS( TA ) ) / BASL 398 210 CONTINUE 399 IF( TB.EQ.ZERO ) 400 $ GO TO 220 401 TB = LOG10( ABS( TB ) ) / BASL 402 220 CONTINUE 403 WORK( I+4*N ) = WORK( I+4*N ) - TA - TB 404 WORK( J+5*N ) = WORK( J+5*N ) - TA - TB 405 230 CONTINUE 406 240 CONTINUE 407* 408 COEF = ONE / REAL( 2*NR ) 409 COEF2 = COEF*COEF 410 COEF5 = HALF*COEF2 411 NRP2 = NR + 2 412 BETA = ZERO 413 IT = 1 414* 415* Start generalized conjugate gradient iteration 416* 417 250 CONTINUE 418* 419 GAMMA = SDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) + 420 $ SDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 ) 421* 422 EW = ZERO 423 EWC = ZERO 424 DO 260 I = ILO, IHI 425 EW = EW + WORK( I+4*N ) 426 EWC = EWC + WORK( I+5*N ) 427 260 CONTINUE 428* 429 GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2 430 IF( GAMMA.EQ.ZERO ) 431 $ GO TO 350 432 IF( IT.NE.1 ) 433 $ BETA = GAMMA / PGAMMA 434 T = COEF5*( EWC-THREE*EW ) 435 TC = COEF5*( EW-THREE*EWC ) 436* 437 CALL SSCAL( NR, BETA, WORK( ILO ), 1 ) 438 CALL SSCAL( NR, BETA, WORK( ILO+N ), 1 ) 439* 440 CALL SAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 ) 441 CALL SAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 ) 442* 443 DO 270 I = ILO, IHI 444 WORK( I ) = WORK( I ) + TC 445 WORK( I+N ) = WORK( I+N ) + T 446 270 CONTINUE 447* 448* Apply matrix to vector 449* 450 DO 300 I = ILO, IHI 451 KOUNT = 0 452 SUM = ZERO 453 DO 290 J = ILO, IHI 454 IF( A( I, J ).EQ.ZERO ) 455 $ GO TO 280 456 KOUNT = KOUNT + 1 457 SUM = SUM + WORK( J ) 458 280 CONTINUE 459 IF( B( I, J ).EQ.ZERO ) 460 $ GO TO 290 461 KOUNT = KOUNT + 1 462 SUM = SUM + WORK( J ) 463 290 CONTINUE 464 WORK( I+2*N ) = REAL( KOUNT )*WORK( I+N ) + SUM 465 300 CONTINUE 466* 467 DO 330 J = ILO, IHI 468 KOUNT = 0 469 SUM = ZERO 470 DO 320 I = ILO, IHI 471 IF( A( I, J ).EQ.ZERO ) 472 $ GO TO 310 473 KOUNT = KOUNT + 1 474 SUM = SUM + WORK( I+N ) 475 310 CONTINUE 476 IF( B( I, J ).EQ.ZERO ) 477 $ GO TO 320 478 KOUNT = KOUNT + 1 479 SUM = SUM + WORK( I+N ) 480 320 CONTINUE 481 WORK( J+3*N ) = REAL( KOUNT )*WORK( J ) + SUM 482 330 CONTINUE 483* 484 SUM = SDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) + 485 $ SDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 ) 486 ALPHA = GAMMA / SUM 487* 488* Determine correction to current iteration 489* 490 CMAX = ZERO 491 DO 340 I = ILO, IHI 492 COR = ALPHA*WORK( I+N ) 493 IF( ABS( COR ).GT.CMAX ) 494 $ CMAX = ABS( COR ) 495 LSCALE( I ) = LSCALE( I ) + COR 496 COR = ALPHA*WORK( I ) 497 IF( ABS( COR ).GT.CMAX ) 498 $ CMAX = ABS( COR ) 499 RSCALE( I ) = RSCALE( I ) + COR 500 340 CONTINUE 501 IF( CMAX.LT.HALF ) 502 $ GO TO 350 503* 504 CALL SAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 ) 505 CALL SAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 ) 506* 507 PGAMMA = GAMMA 508 IT = IT + 1 509 IF( IT.LE.NRP2 ) 510 $ GO TO 250 511* 512* End generalized conjugate gradient iteration 513* 514 350 CONTINUE 515 SFMIN = SLAMCH( 'S' ) 516 SFMAX = ONE / SFMIN 517 LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE ) 518 LSFMAX = INT( LOG10( SFMAX ) / BASL ) 519 DO 360 I = ILO, IHI 520 IRAB = ISAMAX( N-ILO+1, A( I, ILO ), LDA ) 521 RAB = ABS( A( I, IRAB+ILO-1 ) ) 522 IRAB = ISAMAX( N-ILO+1, B( I, ILO ), LDB ) 523 RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) ) 524 LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE ) 525 IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) ) 526 IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB ) 527 LSCALE( I ) = SCLFAC**IR 528 ICAB = ISAMAX( IHI, A( 1, I ), 1 ) 529 CAB = ABS( A( ICAB, I ) ) 530 ICAB = ISAMAX( IHI, B( 1, I ), 1 ) 531 CAB = MAX( CAB, ABS( B( ICAB, I ) ) ) 532 LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE ) 533 JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) ) 534 JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB ) 535 RSCALE( I ) = SCLFAC**JC 536 360 CONTINUE 537* 538* Row scaling of matrices A and B 539* 540 DO 370 I = ILO, IHI 541 CALL SSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA ) 542 CALL SSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB ) 543 370 CONTINUE 544* 545* Column scaling of matrices A and B 546* 547 DO 380 J = ILO, IHI 548 CALL SSCAL( IHI, RSCALE( J ), A( 1, J ), 1 ) 549 CALL SSCAL( IHI, RSCALE( J ), B( 1, J ), 1 ) 550 380 CONTINUE 551* 552 RETURN 553* 554* End of SGGBAL 555* 556 END 557