Annotation of rpl/lapack/lapack/zggbal.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
! 2: $ RSCALE, WORK, INFO )
! 3: *
! 4: * -- LAPACK routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 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) REAL array, dimension (lwork)
! 92: * lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
! 93: * at least 1 when JOB = 'N' or 'P'.
! 94: *
! 95: * INFO (output) INTEGER
! 96: * = 0: successful exit
! 97: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 98: *
! 99: * Further Details
! 100: * ===============
! 101: *
! 102: * See R.C. WARD, Balancing the generalized eigenvalue problem,
! 103: * SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
! 104: *
! 105: * =====================================================================
! 106: *
! 107: * .. Parameters ..
! 108: DOUBLE PRECISION ZERO, HALF, ONE
! 109: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
! 110: DOUBLE PRECISION THREE, SCLFAC
! 111: PARAMETER ( THREE = 3.0D+0, SCLFAC = 1.0D+1 )
! 112: COMPLEX*16 CZERO
! 113: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
! 114: * ..
! 115: * .. Local Scalars ..
! 116: INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
! 117: $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
! 118: $ M, NR, NRP2
! 119: DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
! 120: $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
! 121: $ SFMIN, SUM, T, TA, TB, TC
! 122: COMPLEX*16 CDUM
! 123: * ..
! 124: * .. External Functions ..
! 125: LOGICAL LSAME
! 126: INTEGER IZAMAX
! 127: DOUBLE PRECISION DDOT, DLAMCH
! 128: EXTERNAL LSAME, IZAMAX, DDOT, DLAMCH
! 129: * ..
! 130: * .. External Subroutines ..
! 131: EXTERNAL DAXPY, DSCAL, XERBLA, ZDSCAL, ZSWAP
! 132: * ..
! 133: * .. Intrinsic Functions ..
! 134: INTRINSIC ABS, DBLE, DIMAG, INT, LOG10, MAX, MIN, SIGN
! 135: * ..
! 136: * .. Statement Functions ..
! 137: DOUBLE PRECISION CABS1
! 138: * ..
! 139: * .. Statement Function definitions ..
! 140: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
! 141: * ..
! 142: * .. Executable Statements ..
! 143: *
! 144: * Test the input parameters
! 145: *
! 146: INFO = 0
! 147: IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
! 148: $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
! 149: INFO = -1
! 150: ELSE IF( N.LT.0 ) THEN
! 151: INFO = -2
! 152: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 153: INFO = -4
! 154: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 155: INFO = -6
! 156: END IF
! 157: IF( INFO.NE.0 ) THEN
! 158: CALL XERBLA( 'ZGGBAL', -INFO )
! 159: RETURN
! 160: END IF
! 161: *
! 162: * Quick return if possible
! 163: *
! 164: IF( N.EQ.0 ) THEN
! 165: ILO = 1
! 166: IHI = N
! 167: RETURN
! 168: END IF
! 169: *
! 170: IF( N.EQ.1 ) THEN
! 171: ILO = 1
! 172: IHI = N
! 173: LSCALE( 1 ) = ONE
! 174: RSCALE( 1 ) = ONE
! 175: RETURN
! 176: END IF
! 177: *
! 178: IF( LSAME( JOB, 'N' ) ) THEN
! 179: ILO = 1
! 180: IHI = N
! 181: DO 10 I = 1, N
! 182: LSCALE( I ) = ONE
! 183: RSCALE( I ) = ONE
! 184: 10 CONTINUE
! 185: RETURN
! 186: END IF
! 187: *
! 188: K = 1
! 189: L = N
! 190: IF( LSAME( JOB, 'S' ) )
! 191: $ GO TO 190
! 192: *
! 193: GO TO 30
! 194: *
! 195: * Permute the matrices A and B to isolate the eigenvalues.
! 196: *
! 197: * Find row with one nonzero in columns 1 through L
! 198: *
! 199: 20 CONTINUE
! 200: L = LM1
! 201: IF( L.NE.1 )
! 202: $ GO TO 30
! 203: *
! 204: RSCALE( 1 ) = 1
! 205: LSCALE( 1 ) = 1
! 206: GO TO 190
! 207: *
! 208: 30 CONTINUE
! 209: LM1 = L - 1
! 210: DO 80 I = L, 1, -1
! 211: DO 40 J = 1, LM1
! 212: JP1 = J + 1
! 213: IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
! 214: $ GO TO 50
! 215: 40 CONTINUE
! 216: J = L
! 217: GO TO 70
! 218: *
! 219: 50 CONTINUE
! 220: DO 60 J = JP1, L
! 221: IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
! 222: $ GO TO 80
! 223: 60 CONTINUE
! 224: J = JP1 - 1
! 225: *
! 226: 70 CONTINUE
! 227: M = L
! 228: IFLOW = 1
! 229: GO TO 160
! 230: 80 CONTINUE
! 231: GO TO 100
! 232: *
! 233: * Find column with one nonzero in rows K through N
! 234: *
! 235: 90 CONTINUE
! 236: K = K + 1
! 237: *
! 238: 100 CONTINUE
! 239: DO 150 J = K, L
! 240: DO 110 I = K, LM1
! 241: IP1 = I + 1
! 242: IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
! 243: $ GO TO 120
! 244: 110 CONTINUE
! 245: I = L
! 246: GO TO 140
! 247: 120 CONTINUE
! 248: DO 130 I = IP1, L
! 249: IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
! 250: $ GO TO 150
! 251: 130 CONTINUE
! 252: I = IP1 - 1
! 253: 140 CONTINUE
! 254: M = K
! 255: IFLOW = 2
! 256: GO TO 160
! 257: 150 CONTINUE
! 258: GO TO 190
! 259: *
! 260: * Permute rows M and I
! 261: *
! 262: 160 CONTINUE
! 263: LSCALE( M ) = I
! 264: IF( I.EQ.M )
! 265: $ GO TO 170
! 266: CALL ZSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
! 267: CALL ZSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
! 268: *
! 269: * Permute columns M and J
! 270: *
! 271: 170 CONTINUE
! 272: RSCALE( M ) = J
! 273: IF( J.EQ.M )
! 274: $ GO TO 180
! 275: CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
! 276: CALL ZSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
! 277: *
! 278: 180 CONTINUE
! 279: GO TO ( 20, 90 )IFLOW
! 280: *
! 281: 190 CONTINUE
! 282: ILO = K
! 283: IHI = L
! 284: *
! 285: IF( LSAME( JOB, 'P' ) ) THEN
! 286: DO 195 I = ILO, IHI
! 287: LSCALE( I ) = ONE
! 288: RSCALE( I ) = ONE
! 289: 195 CONTINUE
! 290: RETURN
! 291: END IF
! 292: *
! 293: IF( ILO.EQ.IHI )
! 294: $ RETURN
! 295: *
! 296: * Balance the submatrix in rows ILO to IHI.
! 297: *
! 298: NR = IHI - ILO + 1
! 299: DO 200 I = ILO, IHI
! 300: RSCALE( I ) = ZERO
! 301: LSCALE( I ) = ZERO
! 302: *
! 303: WORK( I ) = ZERO
! 304: WORK( I+N ) = ZERO
! 305: WORK( I+2*N ) = ZERO
! 306: WORK( I+3*N ) = ZERO
! 307: WORK( I+4*N ) = ZERO
! 308: WORK( I+5*N ) = ZERO
! 309: 200 CONTINUE
! 310: *
! 311: * Compute right side vector in resulting linear equations
! 312: *
! 313: BASL = LOG10( SCLFAC )
! 314: DO 240 I = ILO, IHI
! 315: DO 230 J = ILO, IHI
! 316: IF( A( I, J ).EQ.CZERO ) THEN
! 317: TA = ZERO
! 318: GO TO 210
! 319: END IF
! 320: TA = LOG10( CABS1( A( I, J ) ) ) / BASL
! 321: *
! 322: 210 CONTINUE
! 323: IF( B( I, J ).EQ.CZERO ) THEN
! 324: TB = ZERO
! 325: GO TO 220
! 326: END IF
! 327: TB = LOG10( CABS1( B( I, J ) ) ) / BASL
! 328: *
! 329: 220 CONTINUE
! 330: WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
! 331: WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
! 332: 230 CONTINUE
! 333: 240 CONTINUE
! 334: *
! 335: COEF = ONE / DBLE( 2*NR )
! 336: COEF2 = COEF*COEF
! 337: COEF5 = HALF*COEF2
! 338: NRP2 = NR + 2
! 339: BETA = ZERO
! 340: IT = 1
! 341: *
! 342: * Start generalized conjugate gradient iteration
! 343: *
! 344: 250 CONTINUE
! 345: *
! 346: GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
! 347: $ DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
! 348: *
! 349: EW = ZERO
! 350: EWC = ZERO
! 351: DO 260 I = ILO, IHI
! 352: EW = EW + WORK( I+4*N )
! 353: EWC = EWC + WORK( I+5*N )
! 354: 260 CONTINUE
! 355: *
! 356: GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
! 357: IF( GAMMA.EQ.ZERO )
! 358: $ GO TO 350
! 359: IF( IT.NE.1 )
! 360: $ BETA = GAMMA / PGAMMA
! 361: T = COEF5*( EWC-THREE*EW )
! 362: TC = COEF5*( EW-THREE*EWC )
! 363: *
! 364: CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
! 365: CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
! 366: *
! 367: CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
! 368: CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
! 369: *
! 370: DO 270 I = ILO, IHI
! 371: WORK( I ) = WORK( I ) + TC
! 372: WORK( I+N ) = WORK( I+N ) + T
! 373: 270 CONTINUE
! 374: *
! 375: * Apply matrix to vector
! 376: *
! 377: DO 300 I = ILO, IHI
! 378: KOUNT = 0
! 379: SUM = ZERO
! 380: DO 290 J = ILO, IHI
! 381: IF( A( I, J ).EQ.CZERO )
! 382: $ GO TO 280
! 383: KOUNT = KOUNT + 1
! 384: SUM = SUM + WORK( J )
! 385: 280 CONTINUE
! 386: IF( B( I, J ).EQ.CZERO )
! 387: $ GO TO 290
! 388: KOUNT = KOUNT + 1
! 389: SUM = SUM + WORK( J )
! 390: 290 CONTINUE
! 391: WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
! 392: 300 CONTINUE
! 393: *
! 394: DO 330 J = ILO, IHI
! 395: KOUNT = 0
! 396: SUM = ZERO
! 397: DO 320 I = ILO, IHI
! 398: IF( A( I, J ).EQ.CZERO )
! 399: $ GO TO 310
! 400: KOUNT = KOUNT + 1
! 401: SUM = SUM + WORK( I+N )
! 402: 310 CONTINUE
! 403: IF( B( I, J ).EQ.CZERO )
! 404: $ GO TO 320
! 405: KOUNT = KOUNT + 1
! 406: SUM = SUM + WORK( I+N )
! 407: 320 CONTINUE
! 408: WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
! 409: 330 CONTINUE
! 410: *
! 411: SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
! 412: $ DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
! 413: ALPHA = GAMMA / SUM
! 414: *
! 415: * Determine correction to current iteration
! 416: *
! 417: CMAX = ZERO
! 418: DO 340 I = ILO, IHI
! 419: COR = ALPHA*WORK( I+N )
! 420: IF( ABS( COR ).GT.CMAX )
! 421: $ CMAX = ABS( COR )
! 422: LSCALE( I ) = LSCALE( I ) + COR
! 423: COR = ALPHA*WORK( I )
! 424: IF( ABS( COR ).GT.CMAX )
! 425: $ CMAX = ABS( COR )
! 426: RSCALE( I ) = RSCALE( I ) + COR
! 427: 340 CONTINUE
! 428: IF( CMAX.LT.HALF )
! 429: $ GO TO 350
! 430: *
! 431: CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
! 432: CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
! 433: *
! 434: PGAMMA = GAMMA
! 435: IT = IT + 1
! 436: IF( IT.LE.NRP2 )
! 437: $ GO TO 250
! 438: *
! 439: * End generalized conjugate gradient iteration
! 440: *
! 441: 350 CONTINUE
! 442: SFMIN = DLAMCH( 'S' )
! 443: SFMAX = ONE / SFMIN
! 444: LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
! 445: LSFMAX = INT( LOG10( SFMAX ) / BASL )
! 446: DO 360 I = ILO, IHI
! 447: IRAB = IZAMAX( N-ILO+1, A( I, ILO ), LDA )
! 448: RAB = ABS( A( I, IRAB+ILO-1 ) )
! 449: IRAB = IZAMAX( N-ILO+1, B( I, ILO ), LDB )
! 450: RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
! 451: LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
! 452: IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )
! 453: IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
! 454: LSCALE( I ) = SCLFAC**IR
! 455: ICAB = IZAMAX( IHI, A( 1, I ), 1 )
! 456: CAB = ABS( A( ICAB, I ) )
! 457: ICAB = IZAMAX( IHI, B( 1, I ), 1 )
! 458: CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
! 459: LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
! 460: JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) )
! 461: JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
! 462: RSCALE( I ) = SCLFAC**JC
! 463: 360 CONTINUE
! 464: *
! 465: * Row scaling of matrices A and B
! 466: *
! 467: DO 370 I = ILO, IHI
! 468: CALL ZDSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
! 469: CALL ZDSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
! 470: 370 CONTINUE
! 471: *
! 472: * Column scaling of matrices A and B
! 473: *
! 474: DO 380 J = ILO, IHI
! 475: CALL ZDSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
! 476: CALL ZDSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
! 477: 380 CONTINUE
! 478: *
! 479: RETURN
! 480: *
! 481: * End of ZGGBAL
! 482: *
! 483: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zggbal.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack