Annotation of rpl/lapack/lapack/dgebal.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
! 2: *
! 3: * -- LAPACK routine (version 3.2) --
! 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 6: * November 2006
! 7: *
! 8: * .. Scalar Arguments ..
! 9: CHARACTER JOB
! 10: INTEGER IHI, ILO, INFO, LDA, N
! 11: * ..
! 12: * .. Array Arguments ..
! 13: DOUBLE PRECISION A( LDA, * ), SCALE( * )
! 14: * ..
! 15: *
! 16: * Purpose
! 17: * =======
! 18: *
! 19: * DGEBAL balances a general real matrix A. This involves, first,
! 20: * permuting A by a similarity transformation to isolate eigenvalues
! 21: * in the first 1 to ILO-1 and last IHI+1 to N elements on the
! 22: * diagonal; and second, applying a diagonal similarity transformation
! 23: * to rows and columns ILO to IHI to make the rows and columns as
! 24: * close in norm as possible. Both steps are optional.
! 25: *
! 26: * Balancing may reduce the 1-norm of the matrix, and improve the
! 27: * accuracy of the computed eigenvalues and/or eigenvectors.
! 28: *
! 29: * Arguments
! 30: * =========
! 31: *
! 32: * JOB (input) CHARACTER*1
! 33: * Specifies the operations to be performed on A:
! 34: * = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0
! 35: * for i = 1,...,N;
! 36: * = 'P': permute only;
! 37: * = 'S': scale only;
! 38: * = 'B': both permute and scale.
! 39: *
! 40: * N (input) INTEGER
! 41: * The order of the matrix A. N >= 0.
! 42: *
! 43: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
! 44: * On entry, the input matrix A.
! 45: * On exit, A is overwritten by the balanced matrix.
! 46: * If JOB = 'N', A is not referenced.
! 47: * See Further Details.
! 48: *
! 49: * LDA (input) INTEGER
! 50: * The leading dimension of the array A. LDA >= max(1,N).
! 51: *
! 52: * ILO (output) INTEGER
! 53: * IHI (output) INTEGER
! 54: * ILO and IHI are set to integers such that on exit
! 55: * A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
! 56: * If JOB = 'N' or 'S', ILO = 1 and IHI = N.
! 57: *
! 58: * SCALE (output) DOUBLE PRECISION array, dimension (N)
! 59: * Details of the permutations and scaling factors applied to
! 60: * A. If P(j) is the index of the row and column interchanged
! 61: * with row and column j and D(j) is the scaling factor
! 62: * applied to row and column j, then
! 63: * SCALE(j) = P(j) for j = 1,...,ILO-1
! 64: * = D(j) for j = ILO,...,IHI
! 65: * = P(j) for j = IHI+1,...,N.
! 66: * The order in which the interchanges are made is N to IHI+1,
! 67: * then 1 to ILO-1.
! 68: *
! 69: * INFO (output) INTEGER
! 70: * = 0: successful exit.
! 71: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 72: *
! 73: * Further Details
! 74: * ===============
! 75: *
! 76: * The permutations consist of row and column interchanges which put
! 77: * the matrix in the form
! 78: *
! 79: * ( T1 X Y )
! 80: * P A P = ( 0 B Z )
! 81: * ( 0 0 T2 )
! 82: *
! 83: * where T1 and T2 are upper triangular matrices whose eigenvalues lie
! 84: * along the diagonal. The column indices ILO and IHI mark the starting
! 85: * and ending columns of the submatrix B. Balancing consists of applying
! 86: * a diagonal similarity transformation inv(D) * B * D to make the
! 87: * 1-norms of each row of B and its corresponding column nearly equal.
! 88: * The output matrix is
! 89: *
! 90: * ( T1 X*D Y )
! 91: * ( 0 inv(D)*B*D inv(D)*Z ).
! 92: * ( 0 0 T2 )
! 93: *
! 94: * Information about the permutations P and the diagonal matrix D is
! 95: * returned in the vector SCALE.
! 96: *
! 97: * This subroutine is based on the EISPACK routine BALANC.
! 98: *
! 99: * Modified by Tzu-Yi Chen, Computer Science Division, University of
! 100: * California at Berkeley, USA
! 101: *
! 102: * =====================================================================
! 103: *
! 104: * .. Parameters ..
! 105: DOUBLE PRECISION ZERO, ONE
! 106: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 107: DOUBLE PRECISION SCLFAC
! 108: PARAMETER ( SCLFAC = 2.0D+0 )
! 109: DOUBLE PRECISION FACTOR
! 110: PARAMETER ( FACTOR = 0.95D+0 )
! 111: * ..
! 112: * .. Local Scalars ..
! 113: LOGICAL NOCONV
! 114: INTEGER I, ICA, IEXC, IRA, J, K, L, M
! 115: DOUBLE PRECISION C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
! 116: $ SFMIN2
! 117: * ..
! 118: * .. External Functions ..
! 119: LOGICAL LSAME
! 120: INTEGER IDAMAX
! 121: DOUBLE PRECISION DLAMCH
! 122: EXTERNAL LSAME, IDAMAX, DLAMCH
! 123: * ..
! 124: * .. External Subroutines ..
! 125: EXTERNAL DSCAL, DSWAP, XERBLA
! 126: * ..
! 127: * .. Intrinsic Functions ..
! 128: INTRINSIC ABS, MAX, MIN
! 129: * ..
! 130: * .. Executable Statements ..
! 131: *
! 132: * Test the input parameters
! 133: *
! 134: INFO = 0
! 135: IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
! 136: $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
! 137: INFO = -1
! 138: ELSE IF( N.LT.0 ) THEN
! 139: INFO = -2
! 140: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 141: INFO = -4
! 142: END IF
! 143: IF( INFO.NE.0 ) THEN
! 144: CALL XERBLA( 'DGEBAL', -INFO )
! 145: RETURN
! 146: END IF
! 147: *
! 148: K = 1
! 149: L = N
! 150: *
! 151: IF( N.EQ.0 )
! 152: $ GO TO 210
! 153: *
! 154: IF( LSAME( JOB, 'N' ) ) THEN
! 155: DO 10 I = 1, N
! 156: SCALE( I ) = ONE
! 157: 10 CONTINUE
! 158: GO TO 210
! 159: END IF
! 160: *
! 161: IF( LSAME( JOB, 'S' ) )
! 162: $ GO TO 120
! 163: *
! 164: * Permutation to isolate eigenvalues if possible
! 165: *
! 166: GO TO 50
! 167: *
! 168: * Row and column exchange.
! 169: *
! 170: 20 CONTINUE
! 171: SCALE( M ) = J
! 172: IF( J.EQ.M )
! 173: $ GO TO 30
! 174: *
! 175: CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
! 176: CALL DSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
! 177: *
! 178: 30 CONTINUE
! 179: GO TO ( 40, 80 )IEXC
! 180: *
! 181: * Search for rows isolating an eigenvalue and push them down.
! 182: *
! 183: 40 CONTINUE
! 184: IF( L.EQ.1 )
! 185: $ GO TO 210
! 186: L = L - 1
! 187: *
! 188: 50 CONTINUE
! 189: DO 70 J = L, 1, -1
! 190: *
! 191: DO 60 I = 1, L
! 192: IF( I.EQ.J )
! 193: $ GO TO 60
! 194: IF( A( J, I ).NE.ZERO )
! 195: $ GO TO 70
! 196: 60 CONTINUE
! 197: *
! 198: M = L
! 199: IEXC = 1
! 200: GO TO 20
! 201: 70 CONTINUE
! 202: *
! 203: GO TO 90
! 204: *
! 205: * Search for columns isolating an eigenvalue and push them left.
! 206: *
! 207: 80 CONTINUE
! 208: K = K + 1
! 209: *
! 210: 90 CONTINUE
! 211: DO 110 J = K, L
! 212: *
! 213: DO 100 I = K, L
! 214: IF( I.EQ.J )
! 215: $ GO TO 100
! 216: IF( A( I, J ).NE.ZERO )
! 217: $ GO TO 110
! 218: 100 CONTINUE
! 219: *
! 220: M = K
! 221: IEXC = 2
! 222: GO TO 20
! 223: 110 CONTINUE
! 224: *
! 225: 120 CONTINUE
! 226: DO 130 I = K, L
! 227: SCALE( I ) = ONE
! 228: 130 CONTINUE
! 229: *
! 230: IF( LSAME( JOB, 'P' ) )
! 231: $ GO TO 210
! 232: *
! 233: * Balance the submatrix in rows K to L.
! 234: *
! 235: * Iterative loop for norm reduction
! 236: *
! 237: SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
! 238: SFMAX1 = ONE / SFMIN1
! 239: SFMIN2 = SFMIN1*SCLFAC
! 240: SFMAX2 = ONE / SFMIN2
! 241: 140 CONTINUE
! 242: NOCONV = .FALSE.
! 243: *
! 244: DO 200 I = K, L
! 245: C = ZERO
! 246: R = ZERO
! 247: *
! 248: DO 150 J = K, L
! 249: IF( J.EQ.I )
! 250: $ GO TO 150
! 251: C = C + ABS( A( J, I ) )
! 252: R = R + ABS( A( I, J ) )
! 253: 150 CONTINUE
! 254: ICA = IDAMAX( L, A( 1, I ), 1 )
! 255: CA = ABS( A( ICA, I ) )
! 256: IRA = IDAMAX( N-K+1, A( I, K ), LDA )
! 257: RA = ABS( A( I, IRA+K-1 ) )
! 258: *
! 259: * Guard against zero C or R due to underflow.
! 260: *
! 261: IF( C.EQ.ZERO .OR. R.EQ.ZERO )
! 262: $ GO TO 200
! 263: G = R / SCLFAC
! 264: F = ONE
! 265: S = C + R
! 266: 160 CONTINUE
! 267: IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
! 268: $ MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
! 269: F = F*SCLFAC
! 270: C = C*SCLFAC
! 271: CA = CA*SCLFAC
! 272: R = R / SCLFAC
! 273: G = G / SCLFAC
! 274: RA = RA / SCLFAC
! 275: GO TO 160
! 276: *
! 277: 170 CONTINUE
! 278: G = C / SCLFAC
! 279: 180 CONTINUE
! 280: IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
! 281: $ MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
! 282: F = F / SCLFAC
! 283: C = C / SCLFAC
! 284: G = G / SCLFAC
! 285: CA = CA / SCLFAC
! 286: R = R*SCLFAC
! 287: RA = RA*SCLFAC
! 288: GO TO 180
! 289: *
! 290: * Now balance.
! 291: *
! 292: 190 CONTINUE
! 293: IF( ( C+R ).GE.FACTOR*S )
! 294: $ GO TO 200
! 295: IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
! 296: IF( F*SCALE( I ).LE.SFMIN1 )
! 297: $ GO TO 200
! 298: END IF
! 299: IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
! 300: IF( SCALE( I ).GE.SFMAX1 / F )
! 301: $ GO TO 200
! 302: END IF
! 303: G = ONE / F
! 304: SCALE( I ) = SCALE( I )*F
! 305: NOCONV = .TRUE.
! 306: *
! 307: CALL DSCAL( N-K+1, G, A( I, K ), LDA )
! 308: CALL DSCAL( L, F, A( 1, I ), 1 )
! 309: *
! 310: 200 CONTINUE
! 311: *
! 312: IF( NOCONV )
! 313: $ GO TO 140
! 314: *
! 315: 210 CONTINUE
! 316: ILO = K
! 317: IHI = L
! 318: *
! 319: RETURN
! 320: *
! 321: * End of DGEBAL
! 322: *
! 323: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgebal.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack