|
2480
|
1 SUBROUTINE DGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, |
|
|
2 $ RSCALE, WORK, INFO ) |
|
|
3 * |
|
3333
|
4 * -- LAPACK routine (version 3.0) -- |
|
2480
|
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
6 * Courant Institute, Argonne National Lab, and Rice University |
|
|
7 * September 30, 1994 |
|
|
8 * |
|
|
9 * .. Scalar Arguments .. |
|
|
10 CHARACTER JOB |
|
|
11 INTEGER IHI, ILO, INFO, LDA, LDB, N |
|
|
12 * .. |
|
|
13 * .. Array Arguments .. |
|
|
14 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), LSCALE( * ), |
|
|
15 $ RSCALE( * ), WORK( * ) |
|
|
16 * .. |
|
|
17 * |
|
|
18 * Purpose |
|
|
19 * ======= |
|
|
20 * |
|
|
21 * DGGBAL balances a pair of general real matrices (A,B). This |
|
|
22 * involves, first, permuting A and B by similarity transformations to |
|
|
23 * isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N |
|
|
24 * elements on the diagonal; and second, applying a diagonal similarity |
|
|
25 * transformation to rows and columns ILO to IHI to make the rows |
|
|
26 * and columns as close in norm as possible. Both steps are optional. |
|
|
27 * |
|
|
28 * Balancing may reduce the 1-norm of the matrices, and improve the |
|
|
29 * accuracy of the computed eigenvalues and/or eigenvectors in the |
|
|
30 * generalized eigenvalue problem A*x = lambda*B*x. |
|
|
31 * |
|
|
32 * Arguments |
|
|
33 * ========= |
|
|
34 * |
|
|
35 * JOB (input) CHARACTER*1 |
|
|
36 * Specifies the operations to be performed on A and B: |
|
|
37 * = 'N': none: simply set ILO = 1, IHI = N, LSCALE(I) = 1.0 |
|
|
38 * and RSCALE(I) = 1.0 for i = 1,...,N. |
|
|
39 * = 'P': permute only; |
|
|
40 * = 'S': scale only; |
|
|
41 * = 'B': both permute and scale. |
|
|
42 * |
|
|
43 * N (input) INTEGER |
|
|
44 * The order of the matrices A and B. N >= 0. |
|
|
45 * |
|
|
46 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) |
|
|
47 * On entry, the input matrix A. |
|
|
48 * On exit, A is overwritten by the balanced matrix. |
|
|
49 * If JOB = 'N', A is not referenced. |
|
|
50 * |
|
|
51 * LDA (input) INTEGER |
|
|
52 * The leading dimension of the array A. LDA >= max(1,N). |
|
|
53 * |
|
|
54 * B (input/output) DOUBLE PRECISION array, dimension (LDB,N) |
|
|
55 * On entry, the input matrix B. |
|
|
56 * On exit, B is overwritten by the balanced matrix. |
|
|
57 * If JOB = 'N', B is not referenced. |
|
|
58 * |
|
|
59 * LDB (input) INTEGER |
|
|
60 * The leading dimension of the array B. LDB >= max(1,N). |
|
|
61 * |
|
|
62 * ILO (output) INTEGER |
|
|
63 * IHI (output) INTEGER |
|
|
64 * ILO and IHI are set to integers such that on exit |
|
|
65 * A(i,j) = 0 and B(i,j) = 0 if i > j and |
|
|
66 * j = 1,...,ILO-1 or i = IHI+1,...,N. |
|
|
67 * If JOB = 'N' or 'S', ILO = 1 and IHI = N. |
|
|
68 * |
|
|
69 * LSCALE (output) DOUBLE PRECISION array, dimension (N) |
|
|
70 * Details of the permutations and scaling factors applied |
|
|
71 * to the left side of A and B. If P(j) is the index of the |
|
|
72 * row interchanged with row j, and D(j) |
|
|
73 * is the scaling factor applied to row j, then |
|
|
74 * LSCALE(j) = P(j) for J = 1,...,ILO-1 |
|
|
75 * = D(j) for J = ILO,...,IHI |
|
|
76 * = P(j) for J = IHI+1,...,N. |
|
|
77 * The order in which the interchanges are made is N to IHI+1, |
|
|
78 * then 1 to ILO-1. |
|
|
79 * |
|
|
80 * RSCALE (output) DOUBLE PRECISION array, dimension (N) |
|
|
81 * Details of the permutations and scaling factors applied |
|
|
82 * to the right side of A and B. If P(j) is the index of the |
|
|
83 * column interchanged with column j, and D(j) |
|
|
84 * is the scaling factor applied to column j, then |
|
|
85 * LSCALE(j) = P(j) for J = 1,...,ILO-1 |
|
|
86 * = D(j) for J = ILO,...,IHI |
|
|
87 * = P(j) for J = IHI+1,...,N. |
|
|
88 * The order in which the interchanges are made is N to IHI+1, |
|
|
89 * then 1 to ILO-1. |
|
|
90 * |
|
|
91 * WORK (workspace) DOUBLE PRECISION array, dimension (6*N) |
|
|
92 * |
|
|
93 * INFO (output) INTEGER |
|
|
94 * = 0: successful exit |
|
|
95 * < 0: if INFO = -i, the i-th argument had an illegal value. |
|
|
96 * |
|
|
97 * Further Details |
|
|
98 * =============== |
|
|
99 * |
|
|
100 * See R.C. WARD, Balancing the generalized eigenvalue problem, |
|
|
101 * SIAM J. Sci. Stat. Comp. 2 (1981), 141-152. |
|
|
102 * |
|
|
103 * ===================================================================== |
|
|
104 * |
|
|
105 * .. Parameters .. |
|
|
106 DOUBLE PRECISION ZERO, HALF, ONE |
|
|
107 PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) |
|
|
108 DOUBLE PRECISION THREE, SCLFAC |
|
|
109 PARAMETER ( THREE = 3.0D+0, SCLFAC = 1.0D+1 ) |
|
|
110 * .. |
|
|
111 * .. Local Scalars .. |
|
|
112 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1, |
|
|
113 $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN, |
|
|
114 $ M, NR, NRP2 |
|
|
115 DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2, |
|
|
116 $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX, |
|
|
117 $ SFMIN, SUM, T, TA, TB, TC |
|
|
118 * .. |
|
|
119 * .. External Functions .. |
|
|
120 LOGICAL LSAME |
|
|
121 INTEGER IDAMAX |
|
|
122 DOUBLE PRECISION DDOT, DLAMCH |
|
|
123 EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH |
|
|
124 * .. |
|
|
125 * .. External Subroutines .. |
|
|
126 EXTERNAL DAXPY, DSCAL, DSWAP, XERBLA |
|
|
127 * .. |
|
|
128 * .. Intrinsic Functions .. |
|
|
129 INTRINSIC ABS, DBLE, INT, LOG10, MAX, MIN, SIGN |
|
|
130 * .. |
|
|
131 * .. Executable Statements .. |
|
|
132 * |
|
|
133 * Test the input parameters |
|
|
134 * |
|
|
135 INFO = 0 |
|
|
136 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. |
|
|
137 $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN |
|
|
138 INFO = -1 |
|
|
139 ELSE IF( N.LT.0 ) THEN |
|
|
140 INFO = -2 |
|
|
141 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN |
|
|
142 INFO = -4 |
|
|
143 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN |
|
|
144 INFO = -5 |
|
|
145 END IF |
|
|
146 IF( INFO.NE.0 ) THEN |
|
|
147 CALL XERBLA( 'DGGBAL', -INFO ) |
|
|
148 RETURN |
|
|
149 END IF |
|
|
150 * |
|
|
151 K = 1 |
|
|
152 L = N |
|
|
153 * |
|
|
154 * Quick return if possible |
|
|
155 * |
|
|
156 IF( N.EQ.0 ) |
|
|
157 $ RETURN |
|
|
158 * |
|
|
159 IF( LSAME( JOB, 'N' ) ) THEN |
|
|
160 ILO = 1 |
|
|
161 IHI = N |
|
|
162 DO 10 I = 1, N |
|
|
163 LSCALE( I ) = ONE |
|
|
164 RSCALE( I ) = ONE |
|
|
165 10 CONTINUE |
|
|
166 RETURN |
|
|
167 END IF |
|
|
168 * |
|
|
169 IF( K.EQ.L ) THEN |
|
|
170 ILO = 1 |
|
|
171 IHI = 1 |
|
|
172 LSCALE( 1 ) = ONE |
|
|
173 RSCALE( 1 ) = ONE |
|
|
174 RETURN |
|
|
175 END IF |
|
|
176 * |
|
|
177 IF( LSAME( JOB, 'S' ) ) |
|
|
178 $ GO TO 190 |
|
|
179 * |
|
|
180 GO TO 30 |
|
|
181 * |
|
|
182 * Permute the matrices A and B to isolate the eigenvalues. |
|
|
183 * |
|
|
184 * Find row with one nonzero in columns 1 through L |
|
|
185 * |
|
|
186 20 CONTINUE |
|
|
187 L = LM1 |
|
|
188 IF( L.NE.1 ) |
|
|
189 $ GO TO 30 |
|
|
190 * |
|
|
191 RSCALE( 1 ) = 1 |
|
|
192 LSCALE( 1 ) = 1 |
|
|
193 GO TO 190 |
|
|
194 * |
|
|
195 30 CONTINUE |
|
|
196 LM1 = L - 1 |
|
|
197 DO 80 I = L, 1, -1 |
|
|
198 DO 40 J = 1, LM1 |
|
|
199 JP1 = J + 1 |
|
|
200 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) |
|
|
201 $ GO TO 50 |
|
|
202 40 CONTINUE |
|
|
203 J = L |
|
|
204 GO TO 70 |
|
|
205 * |
|
|
206 50 CONTINUE |
|
|
207 DO 60 J = JP1, L |
|
|
208 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) |
|
|
209 $ GO TO 80 |
|
|
210 60 CONTINUE |
|
|
211 J = JP1 - 1 |
|
|
212 * |
|
|
213 70 CONTINUE |
|
|
214 M = L |
|
|
215 IFLOW = 1 |
|
|
216 GO TO 160 |
|
|
217 80 CONTINUE |
|
|
218 GO TO 100 |
|
|
219 * |
|
|
220 * Find column with one nonzero in rows K through N |
|
|
221 * |
|
|
222 90 CONTINUE |
|
|
223 K = K + 1 |
|
|
224 * |
|
|
225 100 CONTINUE |
|
|
226 DO 150 J = K, L |
|
|
227 DO 110 I = K, LM1 |
|
|
228 IP1 = I + 1 |
|
|
229 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) |
|
|
230 $ GO TO 120 |
|
|
231 110 CONTINUE |
|
|
232 I = L |
|
|
233 GO TO 140 |
|
|
234 120 CONTINUE |
|
|
235 DO 130 I = IP1, L |
|
|
236 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) |
|
|
237 $ GO TO 150 |
|
|
238 130 CONTINUE |
|
|
239 I = IP1 - 1 |
|
|
240 140 CONTINUE |
|
|
241 M = K |
|
|
242 IFLOW = 2 |
|
|
243 GO TO 160 |
|
|
244 150 CONTINUE |
|
|
245 GO TO 190 |
|
|
246 * |
|
|
247 * Permute rows M and I |
|
|
248 * |
|
|
249 160 CONTINUE |
|
|
250 LSCALE( M ) = I |
|
|
251 IF( I.EQ.M ) |
|
|
252 $ GO TO 170 |
|
|
253 CALL DSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA ) |
|
|
254 CALL DSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB ) |
|
|
255 * |
|
|
256 * Permute columns M and J |
|
|
257 * |
|
|
258 170 CONTINUE |
|
|
259 RSCALE( M ) = J |
|
|
260 IF( J.EQ.M ) |
|
|
261 $ GO TO 180 |
|
|
262 CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 ) |
|
|
263 CALL DSWAP( L, B( 1, J ), 1, B( 1, M ), 1 ) |
|
|
264 * |
|
|
265 180 CONTINUE |
|
|
266 GO TO ( 20, 90 )IFLOW |
|
|
267 * |
|
|
268 190 CONTINUE |
|
|
269 ILO = K |
|
|
270 IHI = L |
|
|
271 * |
|
|
272 IF( ILO.EQ.IHI ) |
|
|
273 $ RETURN |
|
|
274 * |
|
|
275 IF( LSAME( JOB, 'P' ) ) |
|
|
276 $ RETURN |
|
|
277 * |
|
|
278 * Balance the submatrix in rows ILO to IHI. |
|
|
279 * |
|
|
280 NR = IHI - ILO + 1 |
|
|
281 DO 200 I = ILO, IHI |
|
|
282 RSCALE( I ) = ZERO |
|
|
283 LSCALE( I ) = ZERO |
|
|
284 * |
|
|
285 WORK( I ) = ZERO |
|
|
286 WORK( I+N ) = ZERO |
|
|
287 WORK( I+2*N ) = ZERO |
|
|
288 WORK( I+3*N ) = ZERO |
|
|
289 WORK( I+4*N ) = ZERO |
|
|
290 WORK( I+5*N ) = ZERO |
|
|
291 200 CONTINUE |
|
|
292 * |
|
|
293 * Compute right side vector in resulting linear equations |
|
|
294 * |
|
|
295 BASL = LOG10( SCLFAC ) |
|
|
296 DO 240 I = ILO, IHI |
|
|
297 DO 230 J = ILO, IHI |
|
|
298 TB = B( I, J ) |
|
|
299 TA = A( I, J ) |
|
|
300 IF( TA.EQ.ZERO ) |
|
|
301 $ GO TO 210 |
|
|
302 TA = LOG10( ABS( TA ) ) / BASL |
|
|
303 210 CONTINUE |
|
|
304 IF( TB.EQ.ZERO ) |
|
|
305 $ GO TO 220 |
|
|
306 TB = LOG10( ABS( TB ) ) / BASL |
|
|
307 220 CONTINUE |
|
|
308 WORK( I+4*N ) = WORK( I+4*N ) - TA - TB |
|
|
309 WORK( J+5*N ) = WORK( J+5*N ) - TA - TB |
|
|
310 230 CONTINUE |
|
|
311 240 CONTINUE |
|
|
312 * |
|
|
313 COEF = ONE / DBLE( 2*NR ) |
|
|
314 COEF2 = COEF*COEF |
|
|
315 COEF5 = HALF*COEF2 |
|
|
316 NRP2 = NR + 2 |
|
|
317 BETA = ZERO |
|
|
318 IT = 1 |
|
|
319 * |
|
|
320 * Start generalized conjugate gradient iteration |
|
|
321 * |
|
|
322 250 CONTINUE |
|
|
323 * |
|
|
324 GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) + |
|
|
325 $ DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 ) |
|
|
326 * |
|
|
327 EW = ZERO |
|
|
328 EWC = ZERO |
|
|
329 DO 260 I = ILO, IHI |
|
|
330 EW = EW + WORK( I+4*N ) |
|
|
331 EWC = EWC + WORK( I+5*N ) |
|
|
332 260 CONTINUE |
|
|
333 * |
|
|
334 GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2 |
|
|
335 IF( GAMMA.EQ.ZERO ) |
|
|
336 $ GO TO 350 |
|
|
337 IF( IT.NE.1 ) |
|
|
338 $ BETA = GAMMA / PGAMMA |
|
|
339 T = COEF5*( EWC-THREE*EW ) |
|
|
340 TC = COEF5*( EW-THREE*EWC ) |
|
|
341 * |
|
|
342 CALL DSCAL( NR, BETA, WORK( ILO ), 1 ) |
|
|
343 CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 ) |
|
|
344 * |
|
|
345 CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 ) |
|
|
346 CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 ) |
|
|
347 * |
|
|
348 DO 270 I = ILO, IHI |
|
|
349 WORK( I ) = WORK( I ) + TC |
|
|
350 WORK( I+N ) = WORK( I+N ) + T |
|
|
351 270 CONTINUE |
|
|
352 * |
|
|
353 * Apply matrix to vector |
|
|
354 * |
|
|
355 DO 300 I = ILO, IHI |
|
|
356 KOUNT = 0 |
|
|
357 SUM = ZERO |
|
|
358 DO 290 J = ILO, IHI |
|
|
359 IF( A( I, J ).EQ.ZERO ) |
|
|
360 $ GO TO 280 |
|
|
361 KOUNT = KOUNT + 1 |
|
|
362 SUM = SUM + WORK( J ) |
|
|
363 280 CONTINUE |
|
|
364 IF( B( I, J ).EQ.ZERO ) |
|
|
365 $ GO TO 290 |
|
|
366 KOUNT = KOUNT + 1 |
|
|
367 SUM = SUM + WORK( J ) |
|
|
368 290 CONTINUE |
|
|
369 WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM |
|
|
370 300 CONTINUE |
|
|
371 * |
|
|
372 DO 330 J = ILO, IHI |
|
|
373 KOUNT = 0 |
|
|
374 SUM = ZERO |
|
|
375 DO 320 I = ILO, IHI |
|
|
376 IF( A( I, J ).EQ.ZERO ) |
|
|
377 $ GO TO 310 |
|
|
378 KOUNT = KOUNT + 1 |
|
|
379 SUM = SUM + WORK( I+N ) |
|
|
380 310 CONTINUE |
|
|
381 IF( B( I, J ).EQ.ZERO ) |
|
|
382 $ GO TO 320 |
|
|
383 KOUNT = KOUNT + 1 |
|
|
384 SUM = SUM + WORK( I+N ) |
|
|
385 320 CONTINUE |
|
|
386 WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM |
|
|
387 330 CONTINUE |
|
|
388 * |
|
|
389 SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) + |
|
|
390 $ DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 ) |
|
|
391 ALPHA = GAMMA / SUM |
|
|
392 * |
|
|
393 * Determine correction to current iteration |
|
|
394 * |
|
|
395 CMAX = ZERO |
|
|
396 DO 340 I = ILO, IHI |
|
|
397 COR = ALPHA*WORK( I+N ) |
|
|
398 IF( ABS( COR ).GT.CMAX ) |
|
|
399 $ CMAX = ABS( COR ) |
|
|
400 LSCALE( I ) = LSCALE( I ) + COR |
|
|
401 COR = ALPHA*WORK( I ) |
|
|
402 IF( ABS( COR ).GT.CMAX ) |
|
|
403 $ CMAX = ABS( COR ) |
|
|
404 RSCALE( I ) = RSCALE( I ) + COR |
|
|
405 340 CONTINUE |
|
|
406 IF( CMAX.LT.HALF ) |
|
|
407 $ GO TO 350 |
|
|
408 * |
|
|
409 CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 ) |
|
|
410 CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 ) |
|
|
411 * |
|
|
412 PGAMMA = GAMMA |
|
|
413 IT = IT + 1 |
|
|
414 IF( IT.LE.NRP2 ) |
|
|
415 $ GO TO 250 |
|
|
416 * |
|
|
417 * End generalized conjugate gradient iteration |
|
|
418 * |
|
|
419 350 CONTINUE |
|
|
420 SFMIN = DLAMCH( 'S' ) |
|
|
421 SFMAX = ONE / SFMIN |
|
|
422 LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE ) |
|
|
423 LSFMAX = INT( LOG10( SFMAX ) / BASL ) |
|
|
424 DO 360 I = ILO, IHI |
|
|
425 IRAB = IDAMAX( N-ILO+1, A( I, ILO ), LDA ) |
|
|
426 RAB = ABS( A( I, IRAB+ILO-1 ) ) |
|
|
427 IRAB = IDAMAX( N-ILO+1, B( I, ILO ), LDA ) |
|
|
428 RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) ) |
|
|
429 LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE ) |
|
|
430 IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) ) |
|
|
431 IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB ) |
|
|
432 LSCALE( I ) = SCLFAC**IR |
|
|
433 ICAB = IDAMAX( IHI, A( 1, I ), 1 ) |
|
|
434 CAB = ABS( A( ICAB, I ) ) |
|
|
435 ICAB = IDAMAX( IHI, B( 1, I ), 1 ) |
|
|
436 CAB = MAX( CAB, ABS( B( ICAB, I ) ) ) |
|
|
437 LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE ) |
|
|
438 JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) ) |
|
|
439 JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB ) |
|
|
440 RSCALE( I ) = SCLFAC**JC |
|
|
441 360 CONTINUE |
|
|
442 * |
|
|
443 * Row scaling of matrices A and B |
|
|
444 * |
|
|
445 DO 370 I = ILO, IHI |
|
|
446 CALL DSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA ) |
|
|
447 CALL DSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB ) |
|
|
448 370 CONTINUE |
|
|
449 * |
|
|
450 * Column scaling of matrices A and B |
|
|
451 * |
|
|
452 DO 380 J = ILO, IHI |
|
|
453 CALL DSCAL( IHI, RSCALE( J ), A( 1, J ), 1 ) |
|
|
454 CALL DSCAL( IHI, RSCALE( J ), B( 1, J ), 1 ) |
|
|
455 380 CONTINUE |
|
|
456 * |
|
|
457 RETURN |
|
|
458 * |
|
|
459 * End of DGGBAL |
|
|
460 * |
|
|
461 END |
