![[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|
Return to dgebal.f CVS log | Up to [local] / rpl / lapack / lapack |
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