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