![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.2.2.
1: SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO ) 2: * 3: * -- LAPACK routine (version 3.2.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: * June 2010 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 DISNAN, LSAME 120: INTEGER IDAMAX 121: DOUBLE PRECISION DLAMCH 122: EXTERNAL DISNAN, 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: IF( DISNAN( C+F+CA+R+G+RA ) ) THEN 270: * 271: * Exit if NaN to avoid infinite loop 272: * 273: INFO = -3 274: CALL XERBLA( 'DGEBAL', -INFO ) 275: RETURN 276: END IF 277: F = F*SCLFAC 278: C = C*SCLFAC 279: CA = CA*SCLFAC 280: R = R / SCLFAC 281: G = G / SCLFAC 282: RA = RA / SCLFAC 283: GO TO 160 284: * 285: 170 CONTINUE 286: G = C / SCLFAC 287: 180 CONTINUE 288: IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR. 289: $ MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190 290: F = F / SCLFAC 291: C = C / SCLFAC 292: G = G / SCLFAC 293: CA = CA / SCLFAC 294: R = R*SCLFAC 295: RA = RA*SCLFAC 296: GO TO 180 297: * 298: * Now balance. 299: * 300: 190 CONTINUE 301: IF( ( C+R ).GE.FACTOR*S ) 302: $ GO TO 200 303: IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN 304: IF( F*SCALE( I ).LE.SFMIN1 ) 305: $ GO TO 200 306: END IF 307: IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN 308: IF( SCALE( I ).GE.SFMAX1 / F ) 309: $ GO TO 200 310: END IF 311: G = ONE / F 312: SCALE( I ) = SCALE( I )*F 313: NOCONV = .TRUE. 314: * 315: CALL DSCAL( N-K+1, G, A( I, K ), LDA ) 316: CALL DSCAL( L, F, A( 1, I ), 1 ) 317: * 318: 200 CONTINUE 319: * 320: IF( NOCONV ) 321: $ GO TO 140 322: * 323: 210 CONTINUE 324: ILO = K 325: IHI = L 326: * 327: RETURN 328: * 329: * End of DGEBAL 330: * 331: END