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
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