|
2329
|
1 SUBROUTINE ZGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO ) |
|
|
2 * |
|
7034
|
3 * -- LAPACK routine (version 3.1) -- |
|
|
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
5 * November 2006 |
|
2329
|
6 * |
|
|
7 * .. Scalar Arguments .. |
|
|
8 CHARACTER JOB |
|
|
9 INTEGER IHI, ILO, INFO, LDA, N |
|
|
10 * .. |
|
|
11 * .. Array Arguments .. |
|
|
12 DOUBLE PRECISION SCALE( * ) |
|
|
13 COMPLEX*16 A( LDA, * ) |
|
|
14 * .. |
|
|
15 * |
|
|
16 * Purpose |
|
|
17 * ======= |
|
|
18 * |
|
|
19 * ZGEBAL balances a general complex 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) COMPLEX*16 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 CBAL. |
|
|
98 * |
|
3333
|
99 * Modified by Tzu-Yi Chen, Computer Science Division, University of |
|
|
100 * California at Berkeley, USA |
|
|
101 * |
|
2329
|
102 * ===================================================================== |
|
|
103 * |
|
|
104 * .. Parameters .. |
|
|
105 DOUBLE PRECISION ZERO, ONE |
|
|
106 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) |
|
|
107 DOUBLE PRECISION SCLFAC |
|
7034
|
108 PARAMETER ( SCLFAC = 2.0D+0 ) |
|
2329
|
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 COMPLEX*16 CDUM |
|
|
118 * .. |
|
|
119 * .. External Functions .. |
|
|
120 LOGICAL LSAME |
|
|
121 INTEGER IZAMAX |
|
|
122 DOUBLE PRECISION DLAMCH |
|
|
123 EXTERNAL LSAME, IZAMAX, DLAMCH |
|
|
124 * .. |
|
|
125 * .. External Subroutines .. |
|
|
126 EXTERNAL XERBLA, ZDSCAL, ZSWAP |
|
|
127 * .. |
|
|
128 * .. Intrinsic Functions .. |
|
|
129 INTRINSIC ABS, DBLE, DIMAG, MAX, MIN |
|
|
130 * .. |
|
|
131 * .. Statement Functions .. |
|
|
132 DOUBLE PRECISION CABS1 |
|
|
133 * .. |
|
|
134 * .. Statement Function definitions .. |
|
|
135 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) |
|
|
136 * .. |
|
|
137 * .. Executable Statements .. |
|
|
138 * |
|
|
139 * Test the input parameters |
|
|
140 * |
|
|
141 INFO = 0 |
|
|
142 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. |
|
|
143 $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN |
|
|
144 INFO = -1 |
|
|
145 ELSE IF( N.LT.0 ) THEN |
|
|
146 INFO = -2 |
|
|
147 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN |
|
|
148 INFO = -4 |
|
|
149 END IF |
|
|
150 IF( INFO.NE.0 ) THEN |
|
|
151 CALL XERBLA( 'ZGEBAL', -INFO ) |
|
|
152 RETURN |
|
|
153 END IF |
|
|
154 * |
|
|
155 K = 1 |
|
|
156 L = N |
|
|
157 * |
|
|
158 IF( N.EQ.0 ) |
|
|
159 $ GO TO 210 |
|
|
160 * |
|
|
161 IF( LSAME( JOB, 'N' ) ) THEN |
|
|
162 DO 10 I = 1, N |
|
|
163 SCALE( I ) = ONE |
|
|
164 10 CONTINUE |
|
|
165 GO TO 210 |
|
|
166 END IF |
|
|
167 * |
|
|
168 IF( LSAME( JOB, 'S' ) ) |
|
|
169 $ GO TO 120 |
|
|
170 * |
|
|
171 * Permutation to isolate eigenvalues if possible |
|
|
172 * |
|
|
173 GO TO 50 |
|
|
174 * |
|
|
175 * Row and column exchange. |
|
|
176 * |
|
|
177 20 CONTINUE |
|
|
178 SCALE( M ) = J |
|
|
179 IF( J.EQ.M ) |
|
|
180 $ GO TO 30 |
|
|
181 * |
|
|
182 CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 ) |
|
|
183 CALL ZSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA ) |
|
|
184 * |
|
|
185 30 CONTINUE |
|
|
186 GO TO ( 40, 80 )IEXC |
|
|
187 * |
|
|
188 * Search for rows isolating an eigenvalue and push them down. |
|
|
189 * |
|
|
190 40 CONTINUE |
|
|
191 IF( L.EQ.1 ) |
|
|
192 $ GO TO 210 |
|
|
193 L = L - 1 |
|
|
194 * |
|
|
195 50 CONTINUE |
|
|
196 DO 70 J = L, 1, -1 |
|
|
197 * |
|
|
198 DO 60 I = 1, L |
|
|
199 IF( I.EQ.J ) |
|
|
200 $ GO TO 60 |
|
|
201 IF( DBLE( A( J, I ) ).NE.ZERO .OR. DIMAG( A( J, I ) ).NE. |
|
|
202 $ ZERO )GO TO 70 |
|
|
203 60 CONTINUE |
|
|
204 * |
|
|
205 M = L |
|
|
206 IEXC = 1 |
|
|
207 GO TO 20 |
|
|
208 70 CONTINUE |
|
|
209 * |
|
|
210 GO TO 90 |
|
|
211 * |
|
|
212 * Search for columns isolating an eigenvalue and push them left. |
|
|
213 * |
|
|
214 80 CONTINUE |
|
|
215 K = K + 1 |
|
|
216 * |
|
|
217 90 CONTINUE |
|
|
218 DO 110 J = K, L |
|
|
219 * |
|
|
220 DO 100 I = K, L |
|
|
221 IF( I.EQ.J ) |
|
|
222 $ GO TO 100 |
|
|
223 IF( DBLE( A( I, J ) ).NE.ZERO .OR. DIMAG( A( I, J ) ).NE. |
|
|
224 $ ZERO )GO TO 110 |
|
|
225 100 CONTINUE |
|
|
226 * |
|
|
227 M = K |
|
|
228 IEXC = 2 |
|
|
229 GO TO 20 |
|
|
230 110 CONTINUE |
|
|
231 * |
|
|
232 120 CONTINUE |
|
|
233 DO 130 I = K, L |
|
|
234 SCALE( I ) = ONE |
|
|
235 130 CONTINUE |
|
|
236 * |
|
|
237 IF( LSAME( JOB, 'P' ) ) |
|
|
238 $ GO TO 210 |
|
|
239 * |
|
|
240 * Balance the submatrix in rows K to L. |
|
|
241 * |
|
|
242 * Iterative loop for norm reduction |
|
|
243 * |
|
|
244 SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' ) |
|
|
245 SFMAX1 = ONE / SFMIN1 |
|
|
246 SFMIN2 = SFMIN1*SCLFAC |
|
|
247 SFMAX2 = ONE / SFMIN2 |
|
|
248 140 CONTINUE |
|
|
249 NOCONV = .FALSE. |
|
|
250 * |
|
|
251 DO 200 I = K, L |
|
|
252 C = ZERO |
|
|
253 R = ZERO |
|
|
254 * |
|
|
255 DO 150 J = K, L |
|
|
256 IF( J.EQ.I ) |
|
|
257 $ GO TO 150 |
|
|
258 C = C + CABS1( A( J, I ) ) |
|
|
259 R = R + CABS1( A( I, J ) ) |
|
|
260 150 CONTINUE |
|
|
261 ICA = IZAMAX( L, A( 1, I ), 1 ) |
|
|
262 CA = ABS( A( ICA, I ) ) |
|
|
263 IRA = IZAMAX( N-K+1, A( I, K ), LDA ) |
|
|
264 RA = ABS( A( I, IRA+K-1 ) ) |
|
|
265 * |
|
|
266 * Guard against zero C or R due to underflow. |
|
|
267 * |
|
|
268 IF( C.EQ.ZERO .OR. R.EQ.ZERO ) |
|
|
269 $ GO TO 200 |
|
|
270 G = R / SCLFAC |
|
|
271 F = ONE |
|
|
272 S = C + R |
|
|
273 160 CONTINUE |
|
|
274 IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR. |
|
|
275 $ MIN( R, G, RA ).LE.SFMIN2 )GO TO 170 |
|
|
276 F = F*SCLFAC |
|
|
277 C = C*SCLFAC |
|
|
278 CA = CA*SCLFAC |
|
|
279 R = R / SCLFAC |
|
|
280 G = G / SCLFAC |
|
|
281 RA = RA / SCLFAC |
|
|
282 GO TO 160 |
|
|
283 * |
|
|
284 170 CONTINUE |
|
|
285 G = C / SCLFAC |
|
|
286 180 CONTINUE |
|
|
287 IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR. |
|
|
288 $ MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190 |
|
|
289 F = F / SCLFAC |
|
|
290 C = C / SCLFAC |
|
|
291 G = G / SCLFAC |
|
|
292 CA = CA / SCLFAC |
|
|
293 R = R*SCLFAC |
|
|
294 RA = RA*SCLFAC |
|
|
295 GO TO 180 |
|
|
296 * |
|
|
297 * Now balance. |
|
|
298 * |
|
|
299 190 CONTINUE |
|
|
300 IF( ( C+R ).GE.FACTOR*S ) |
|
|
301 $ GO TO 200 |
|
|
302 IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN |
|
|
303 IF( F*SCALE( I ).LE.SFMIN1 ) |
|
|
304 $ GO TO 200 |
|
|
305 END IF |
|
|
306 IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN |
|
|
307 IF( SCALE( I ).GE.SFMAX1 / F ) |
|
|
308 $ GO TO 200 |
|
|
309 END IF |
|
|
310 G = ONE / F |
|
|
311 SCALE( I ) = SCALE( I )*F |
|
|
312 NOCONV = .TRUE. |
|
|
313 * |
|
|
314 CALL ZDSCAL( N-K+1, G, A( I, K ), LDA ) |
|
|
315 CALL ZDSCAL( L, F, A( 1, I ), 1 ) |
|
|
316 * |
|
|
317 200 CONTINUE |
|
|
318 * |
|
|
319 IF( NOCONV ) |
|
|
320 $ GO TO 140 |
|
|
321 * |
|
|
322 210 CONTINUE |
|
|
323 ILO = K |
|
|
324 IHI = L |
|
|
325 * |
|
|
326 RETURN |
|
|
327 * |
|
|
328 * End of ZGEBAL |
|
|
329 * |
|
|
330 END |
